Pixel spectrum

Spectrum

A class representing a spectrum, with wavelength and flux arrays, and methods to manipulate them.

Parameters:

Name	Type	Description	Default
wave	array - like	Wavelength array of the spectrum [Angstrom].	required
flux	array - like	Flux array of the spectrum, same shape as wave.	required
step	float	Step size for the CCF, in km/s. Default is 0.5 km/s.	0.5

Source code in SOAP/classes.py

class Spectrum:
    """ 
    A class representing a spectrum, with wavelength and flux arrays, and methods to manipulate them.

    Args:
        wave (array-like):
            Wavelength array of the spectrum [Angstrom].
        flux (array-like):
            Flux array of the spectrum, same shape as wave.
        step (float, optional):
            Step size for the CCF, in km/s. Default is 0.5 km/s.

    """
    # is this a CCF?
    _ccf: bool = False

    def __init__(self, wave: np.ndarray, flux: np.ndarray, step=0.5) -> None:
        self.wave = wave
        self.flux = flux
        self._rv_units = kms
        self.step = step
        self._vrot = 0.0

    @property
    def n(self):
        return self.wave.size

    @property
    def vrot(self):
        """Rotation velocity of the star"""
        return self._vrot

    @property
    def n_v(self):
        """
        Total number of RV bins in the CCF after considering a stellar rotation
        velocity equal to self.vrot.
        """
        # The CCF is assumed to have a rotational velocity equal to 0 because it
        # is taken in the stellar disk center. To consider rotation when
        # simulating the emerging spectrum on the limb of the star (velocity
        # different from 0), we have to extrapolate outside of the range
        # [-self.width, self.width]. This is done by calculating self.n_v which
        # is the number of additional RV bins required to consider a stellar
        # rotation of star.vrot

        # self.n_v must by an odd integer so that we have as many values on the
        # positive side of the CCF as on the negative one (because 0 is present)
        r = np.round(((1.1 * self.vrot) / self.step) * 2)
        if r % 2 == 0:
            return int(self.n + r)
        else:
            return int(self.n + r + 1)

    @property
    def v_interval(self):
        # self.v_interval gives the difference in radial velocity
        # between the minimum and maximum values of the CCF,
        # once the extrapolation is done
        # self.n_v gives the number of points of the extrapolated CCF.
        # (self.n_v-1) gives the number of intervals,
        # which is then multiplied by the sampling of the CCF
        # self.step to give self.v_interval
        return self.step * (self.n_v - 1) / 2.0

    def __call__(self, *args, **kwargs):
        return self.flux

    def plot(self, ax=None, thin=100, **kwargs):
        import matplotlib.pyplot as plt

        if ax is None:
            fig, ax = plt.subplots()
        else:
            fig = ax.figure

        if self.flux.ndim == 1:
            kwargs.setdefault("color", "k")

        if self.flux.shape[0] > 50000:
            ax.plot(self.wave[::thin], self.flux[::thin], **kwargs)
        else:
            ax.plot(self.wave, self.flux, **kwargs)
        ax.set(xlabel=r"wavelength [$\AA$]", ylabel="flux")
        return fig, ax

    def doppler_shift(self, rv: float, inplace=False):
        shifted = doppler_shift(self.wave, self.flux, rv)
        if inplace:
            self.flux = shifted
        else:
            return shifted

    def interpolate_to(self, wave, inplace=False, **kwargs):
        # interpolate this spectrum to wavelengths `wave`
        interpolated = np.interp(wave, self.wave, self.flux)
        if inplace:
            self.flux = interpolated
            self.wave = wave
        else:
            return interpolated

    def convolve_to(self, R, inplace=False):
        if self.flux.ndim == 2:
            flux = np.zeros_like(self.flux)
            for i in range(self.flux.shape[1]):
                _, f = ip_convolution(self.wave, self.flux[:, i], R)
                flux[:, i] = f
        else:
            _, flux = ip_convolution(self.wave, self.flux, R)

        if inplace:
            self.flux = flux
        else:
            return flux

    def resample(self, every, inplace=False):
        wave = self.wave[::every]
        flux = self.flux[::every]
        if inplace:
            self.wave = wave
            self.flux = flux
        else:
            return wave, flux

    def to_numba(self):
        return SpectrumNumba(self.wave.astype(float), self.flux.astype(float))

    @classmethod
    def random(cls, min_wave: float = 6000.0, max_wave: float = 6010.0):
        # define high-level parameters, especially including spectrograph parameters
        R = 135000  # resolution
        SNR = 100.0  # s/n ratio in the continuum
        continuum_ivar = SNR**2  # inverse variance of the noise in the continuum
        sigma_x = 1.0 / R  # LSF sigma in x units
        Delta_x = 1.0 / (3.0 * R)  # pixel spacing
        x_min = np.log(min_wave)  # minimum ln wavelength
        x_max = np.log(max_wave)  # maximum ln wavelength
        lines_per_x = (
            2.0e4  # mean density (Poisson rate) of lines per unit ln wavelength
        )
        ew_max_x = 3.0e-5  # maximum equivalent width in x units
        ew_power = 5.0  # power parameter in EW maker
        # set up the true line list for the true spectral model
        x_margin = 1.0e6 / c.value  # hoping no velocities are bigger than 1000 km/s
        x_range = (
            x_max - x_min + 2.0 * x_margin
        )  # make lines in a bigger x range than the data range
        nlines = np.random.poisson(
            x_range * lines_per_x
        )  # set the total number of lines
        line_xs = (x_min - x_margin) + x_range * np.random.uniform(size=nlines)
        # give those lines equivalent widths from a power-law distribution
        line_ews = ew_max_x * np.random.uniform(size=nlines) ** ew_power

        xs = np.arange(x_min, x_max, Delta_x)
        true_doppler = 0.0
        ys, y_ivars = Spectrum._noisy_true_spectrum(
            xs, true_doppler, continuum_ivar, line_xs, line_ews, sigma_x
        )
        # y_ivars_empirical = Spectrum._ivar(ys, continuum_ivar)
        return cls(np.exp(xs), ys)

    @staticmethod
    def _oned_gaussian(dxs, sigma):
        return np.exp(-0.5 * dxs**2 / sigma**2) / (sqrt2pi * sigma)

    @staticmethod
    def _true_spectrum(xs, doppler, lxs, ews, sigma):
        g = Spectrum._oned_gaussian(xs[:, None] - doppler - lxs[None, :], sigma)
        return np.exp(-1.0 * np.sum(ews[None, :] * g, axis=1))

    @staticmethod
    def _ivar(ys, continuum_ivar):
        return continuum_ivar / ys

    @staticmethod
    def _noisy_true_spectrum(xs, doppler, continuum_ivar, lxs, ews, sigma):
        ys_true = Spectrum._true_spectrum(xs, doppler, lxs, ews, sigma)
        y_ivars = Spectrum._ivar(ys_true, continuum_ivar)
        ys = ys_true + np.random.normal(size=xs.shape) / np.sqrt(y_ivars)
        return ys, y_ivars

n_v property

Total number of RV bins in the CCF after considering a stellar rotation velocity equal to self.vrot.

vrot property

Rotation velocity of the star

solarFTS

Bases: Spectrum

Solar spectrum obtained with the FTS spectrograph, for the quiet Sun (Wallace+1998) or a sunspot umbra(Wallace+1995). Source: https://nso.edu/data/historical-archive/

Parameters:

Name	Type	Description	Default
spot	bool	If True, get the spectrum for a sunspot. If False, get the spectrum for the quiet Sun. Default is False.	False
wave_range	tuple of float	Wavelength range to get the spectrum, in Angstrom. Default is (3500, 7500).	(3500, 7500)
resolution	float	If not None, convolve the spectrum to this resolution. Default is None (no convolution).	None
air_wave	bool	If True, convert the wavelengths from vacuum to air. Default is True.	True

Source code in SOAP/classes.py

class solarFTS(Spectrum):
    """
    Solar spectrum obtained with the FTS spectrograph, for the quiet Sun (Wallace+1998) or a sunspot umbra(Wallace+1995).
    Source: https://nso.edu/data/historical-archive/

    Args:
        spot (bool, optional):
            If True, get the spectrum for a sunspot. If False, get the spectrum for the quiet Sun. Default is False.
        wave_range (tuple of float, optional):
            Wavelength range to get the spectrum, in Angstrom. Default is (3500, 7500).
        resolution (float, optional):
            If not None, convolve the spectrum to this resolution. Default is None (no convolution).
        air_wave (bool, optional):
            If True, convert the wavelengths from vacuum to air. Default is True.
    """

    def __init__(
        self, spot=False, wave_range=(3500, 7500), resolution=None, air_wave=True
    ) -> None:
        here = os.path.dirname(__file__)
        if spot:
            file = os.path.join(here, "../data/sunspot_FTS.npy")
        else:
            file = os.path.join(here, "../data/quiet_sun_FTS.npy")

        wave, flux = np.load(file)

        if air_wave:
            wave = vacuum_to_air(wave)

        mask = (wave > wave_range[0]) & (wave < wave_range[1])
        wave = wave[mask]
        flux = flux[mask]

        if resolution is not None:
            from scipy.ndimage.filters import gaussian_filter1d

            # resolution R = λ / Δλ => Δλ = λ / R
            inst_profile_sig = np.median(wave) / resolution
            # inst_profile_sig = inst_profile_FWHM / (2 * np.sqrt(2 * np.log(2)))
            flux = gaussian_filter1d(flux, inst_profile_sig)

        super().__init__(wave, flux)
        self.wave_range = wave_range

Spec_mu

Bases: Spectrum

A class representing local input spectra defined on a grid of μ values, with wavelength and flux arrays, and methods to manipulate them. The grid of spectra will be linearly interpolated to get the spectrum at any μ value between 0 and 1. Near the limb (μ close to 0), the spectra are set to match the smallest μ value spectrum. Near the disk center (μ close to 1), the spectra can be extrapolated by taking the spectrum at the largest μ value. The spectra must be defined on the same wavelength grid for all μ values.

Parameters:

Name	Type	Description	Default
mu_array	array - like	Array of μ values where the spectrum is defined, between 0 and 1.	required
wavelength	array - like	Wavelength array of the spectrum [Angstrom]. Must be the same for all μ values.	required
spectra	array - like	Array of spectra for each μ value.	required
wave_range	tuple of float	Wavelength range to get the spectrum, in Angstrom. Default is (4198, 7998).	required
air_wave	bool	If True, convert the wavelengths from vacuum to air. Default is True.	True

Attributes:

Name	Type	Description
μ	array - like	Array of μ values where the spectrum is defined, between 0 and 1.
wave	array - like	Wavelength array of the spectrum [Angstrom].
flux	array - like	Flux array of the spectrum, same shape as wave.

Outputs: A Spectrum object that can be called with a μ value to get the corresponding spectrum, and that can be interpolated to other wavelengths or convolved to other resolutions.

Source code in SOAP/classes.py

class Spec_mu(Spectrum):
    """
    A class representing local input spectra defined on a grid of μ values, with wavelength and flux arrays, and methods to manipulate them. 
    The grid of spectra will be linearly interpolated to get the spectrum at any μ value between 0 and 1. 
    Near the limb (μ close to 0), the spectra are set to match the smallest μ value spectrum. 
    Near the disk center (μ close to 1), the spectra can be extrapolated by taking the spectrum at the largest μ value. 
    The spectra must be defined on the same wavelength grid for all μ values.

    Args:
        mu_array (array-like):
            Array of μ values where the spectrum is defined, between 0 and 1.
        wavelength (array-like):
            Wavelength array of the spectrum [Angstrom]. Must be the same for all μ values.
        spectra (array-like):
            Array of spectra for each μ value.
        wave_range (tuple of float, optional):
            Wavelength range to get the spectrum, in Angstrom. Default is (4198, 7998).
        air_wave (bool, optional):
            If True, convert the wavelengths from vacuum to air. Default is True.

    Attributes:
        μ (array-like):
            Array of μ values where the spectrum is defined, between 0 and 1.
        wave (array-like):
            Wavelength array of the spectrum [Angstrom].
        flux (array-like):
            Flux array of the spectrum, same shape as wave.
    Outputs:
        A Spectrum object that can be called with a μ value to get the corresponding spectrum, and that can be interpolated to other wavelengths or convolved to other resolutions.
    """
    def __init__(self,
                 mu_array,
                 wavelength,
                 spectra,
                 wave_range,
                 air_wave=True
                 ):
        self.μ = mu_array
        wave=wavelength
        flux=spectra

        # wave cuts
        mask = (wave > wave_range[0]) & (wave < wave_range[1])
        wave = wave[mask]
        flux = flux[mask]

        # remove NaNs
        nan_mask = np.zeros_like(wave, dtype=bool)
        for i,f in enumerate(flux.T):
            flux.T[i]=f/np.max(flux.T[i])

        for f in flux.T:
            nan_mask = np.logical_or(nan_mask, np.isnan(f))
        wave = wave[~nan_mask]
        flux = flux[~nan_mask, :]

        # conversions
        wave = wave.astype(float)
        flux = flux.astype(float)

        super().__init__(wave, flux)
        self.wave_range = wave_range
        self.interpolate()

    def interpolate(self):
        self.interp = RegularGridInterpolator((self.μ, self.wave), self.flux.T)

    def __call__(self, mu: float):
        xi = np.c_[np.full(self.wave.size, mu), self.wave]
        return self.interp(xi)

    def to_numba(self):
        return SpectrumNumbaInterpolated(
            self.wave.astype(float), self.flux.astype(float),self.μ.astype(float)
        )

solarIAGatlas

Bases: Spectrum

A class representing the Goettingen Solar Atlas, with wavelength and flux arrays, and methods to manipulate them. The grid of spectra is linearly interpolated to get the spectrum at any μ value between 0 and 1. Near the limb (μ close to 0), the spectra are set to match the smallest μ value spectrum. Near the disk center (μ close to 1), the spectra can be extrapolated by taking the spectrum at the largest μ value. Default: μ = [0.2, 0.3, 0.35, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 0.95, 0.97, 0.98, 0.99, 1] Source: http://www.astro.physik.uni-goettingen.de/research/solar-lib/

Parameters:

Name	Type	Description	Default
wave_range	tuple of float	Wavelength range to get the spectrum, in Angstrom. Default is (4198, 7998).	(4198, 7998)
air_wave	bool	If True, convert the wavelengths from vacuum to air. Default is True.	True

Attributes:

Name	Type	Description
μ	array - like	Array of μ values where the spectrum is defined, between 0 and 1.
wave	array - like	Wavelength array of the spectrum [Angstrom].
flux	array - like	Flux array of the spectrum, same shape as wave.

Outputs

A Spectrum object that can be called with a μ value to get the corresponding spectrum, and that can be interpolated to other wavelengths or convolved to other resolutions.

Source code in SOAP/classes.py

class solarIAGatlas(Spectrum):
    """
    A class representing the Goettingen Solar Atlas, with wavelength and flux arrays, and methods to manipulate them.
    The grid of spectra is linearly interpolated to get the spectrum at any μ value between 0 and 1. 
    Near the limb (μ close to 0), the spectra are set to match the smallest μ value spectrum. 
    Near the disk center (μ close to 1), the spectra can be extrapolated by taking the spectrum at the largest μ value.
    Default: μ = [0.2, 0.3, 0.35, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 0.95, 0.97, 0.98, 0.99, 1]
    Source: http://www.astro.physik.uni-goettingen.de/research/solar-lib/

    Args:
        wave_range (tuple of float, optional):
            Wavelength range to get the spectrum, in Angstrom. Default is (4198, 7998).
        air_wave (bool, optional):
            If True, convert the wavelengths from vacuum to air. Default is True.

    Attributes:
        μ (array-like):
            Array of μ values where the spectrum is defined, between 0 and 1.
        wave (array-like):
            Wavelength array of the spectrum [Angstrom].
        flux (array-like):
            Flux array of the spectrum, same shape as wave.

    Outputs:
        A Spectrum object that can be called with a μ value to get the corresponding spectrum, and that can be interpolated to other wavelengths or convolved to other resolutions.     
    """
    def __init__(self, wave_range=(4198, 7998), air_wave=True):
        here = os.path.dirname(__file__)
        files = sorted(glob(os.path.join(here, "../data/IAGatlas/*.fits")))

        # need to download the atlas
        if len(files) < 14:
            from .utils import download_goettingen_solar_atlas

            download_goettingen_solar_atlas()

        self.μ = np.array([float(f.split("mu")[1][:-5]) for f in files])
        self.atlas_files = files

        wave_file = os.path.join(here, "../data/IAGatlas/wave.npy")
        if os.path.exists(wave_file):
            wave = np.load(wave_file)
        else:
            wave = fits.getdata(files[0])[0]
            np.save(wave_file, wave)

        flux_file = os.path.join(here, "../data/IAGatlas/flux.npy")
        if os.path.exists(flux_file):
            flux = np.load(flux_file)
        else:
            flux = np.array([fits.getdata(f)[1] for f in files]).T
            np.save(flux_file, flux)

        if air_wave:
            wave = vacuum_to_air(wave)

        # wave cuts
        mask = (wave > wave_range[0]) & (wave < wave_range[1])
        wave = wave[mask]
        flux = flux[mask]

        # remove NaNs
        nan_mask = np.zeros_like(wave, dtype=bool)
        for f in flux.T:
            nan_mask = np.logical_or(nan_mask, np.isnan(f))
        wave = wave[~nan_mask]
        flux = flux[~nan_mask, :]

        # conversions
        wave = wave.astype(float)
        flux = flux.astype(float)

        super().__init__(wave, flux)
        self.wave_range = wave_range
        self.interpolate()

    def interpolate(self):
        self.interp = RegularGridInterpolator((self.μ, self.wave), self.flux.T)

    def __call__(self, mu: float):
        xi = np.c_[np.full(self.wave.size, mu), self.wave]
        return self.interp(xi)

    def to_numba(self):
        return SpectrumNumbaInterpolated(
            self.wave.astype(float), self.flux.astype(float),self.μ.astype(float)
        )

PHOENIX

Bases: Spectrum

A class representing a PHOENIX spectrum, with wavelength and flux arrays, and methods to manipulate them.

Parameters:

Name	Type	Description	Default
ar	bool	If True, get the spectrum for a spot (with a contrast in temperature). If False, get the spectrum for the quiet photosphere. Default is False.	False
wave_range	tuple of float	Wavelength range to get the spectrum, in Angstrom. Default is (3500, 7500).	(3500, 7500)
resolution	float	If not None, convolve the spectrum to this resolution. Default is None (no convolution).	None
teff	float	Effective temperature of the star in Kelvin. Default is 5780 K (solar).	5780
logg	float	Surface gravity of the star in cgs. Default is 4.4 (solar).	4.4
Z	float	Metallicity of the star. Default is 0.012 (solar).	0.012
St_alpha	float	Alpha element abundance of the star. Default is 0 (solar).	0
contrast	float	Temperature contrast of the spot in Kelvin. Default is 0 K (no contrast).	0
normalize	bool	If True, normalize the spectrum by fitting a linear continuum to the 10% highest flux points in the spectrum. Default is False.	False
cache	bool	If True, save the spectrum to a file and load it if the same parameters are used again. Default is True.	True

Attributes:

Name	Type	Description
wave	array - like	Wavelength array of the spectrum [Angstrom].
flux	array - like	Flux array of the spectrum, same shape as wave.

Source code in SOAP/classes.py

class PHOENIX(Spectrum):
    """
    A class representing a PHOENIX spectrum, with wavelength and flux arrays, and methods to manipulate them.

    Args:
        ar (bool, optional):
            If True, get the spectrum for a spot (with a contrast in temperature). If False, get the spectrum for the quiet photosphere. Default is False.
        wave_range (tuple of float, optional):
            Wavelength range to get the spectrum, in Angstrom. Default is (3500, 7500).
        resolution (float, optional):
            If not None, convolve the spectrum to this resolution. Default is None (no convolution).
        teff (float, optional):
            Effective temperature of the star in Kelvin. Default is 5780 K (solar).
        logg (float, optional):
            Surface gravity of the star in cgs. Default is 4.4 (solar).
        Z (float, optional):
            Metallicity of the star. Default is 0.012 (solar).
        St_alpha (float, optional):
            Alpha element abundance of the star. Default is 0 (solar).
        contrast (float, optional):
            Temperature contrast of the spot in Kelvin. Default is 0 K (no contrast).
        normalize (bool, optional):
            If True, normalize the spectrum by fitting a linear continuum to the 10% highest flux points in the spectrum. Default is False.
        cache (bool, optional):
            If True, save the spectrum to a file and load it if the same parameters are used again. Default is True.

    Attributes:
        wave (array-like):
            Wavelength array of the spectrum [Angstrom].
        flux (array-like):
            Flux array of the spectrum, same shape as wave.
    """
    def __init__(
        self,
        ar=False,
        wave_range=(3500, 7500),
        resolution=None,
        teff=5780,
        logg=4.4,
        Z=0.012,
        St_alpha=0,
        contrast=0,
        normalize=False,
        cache=True,
    ) -> None:
        if cache == True:
            read_dict = None
            directory = Path().resolve()
            matching_files = [
                f.name
                for f in directory.iterdir()
                if f.is_file()
                and f.name.startswith("spectrum")
                and f.name.endswith("npz")
            ]
            prop_dict = {
                "ar": ar,
                "wave_range": wave_range,
                "resolution": resolution,
                "teff": teff,
                "logg": logg,
                "Z": Z,
                "St_alpha": St_alpha,
                "contrast": contrast,
                "normalize": normalize,
                "cache": cache,
            }
            if len(matching_files) == 0:
                pass
            else:
                for saved_spec in matching_files:
                    data = np.load(saved_spec, allow_pickle=True)
                    flux_read = data["flux"]
                    wave_read = data["wavelength"]
                    read_dict = data["prop_dict"]
                    if prop_dict == read_dict:
                        flux = flux_read
                        wave = wave_read
                        super().__init__(wave, flux)
                        return None
                    else:
                        pass
        else:
            pass
        if ar:
            spectrum = get_spectrum(
                T_eff=teff + contrast, log_g=logg, Z=Z, alpha=St_alpha
            )
        else:
            spectrum = get_spectrum(T_eff=teff, log_g=logg, Z=Z, alpha=St_alpha)
        # Read the spectrum
        wave = vacuum_to_air((spectrum.wavelength.value).astype(float))
        flux = (spectrum.flux.value).astype(float)

        mask_0 = (wave > wave_range[0]) & (wave < wave_range[1])

        flux = flux[mask_0]
        wave = wave[mask_0]
        if normalize:
            n = len(flux)
            # Compute 10% of the points in the spectrum
            n_10percent = int(n * 0.1)

            # From the 10% of the points in the spectrum compute the 99 percentile in flux
            mask_thresh_1 = flux[:n_10percent] >= np.percentile(flux[:n_10percent], 99)
            mask_thresh_2 = flux[-n_10percent:] >= np.percentile(
                flux[-n_10percent:], 99
            )

            # Compute the respective points in flux and wavelength
            percent_continuum_flux = np.concatenate(
                (flux[:n_10percent][mask_thresh_1], flux[-n_10percent:][mask_thresh_2])
            )
            percent_continuum_wave = np.concatenate(
                (wave[:n_10percent][mask_thresh_1], wave[-n_10percent:][mask_thresh_2])
            )

            # Fit the points with a first degree polynomial
            slope_coefs = np.polyfit(percent_continuum_wave, percent_continuum_flux, 1)
            slope_spectrum = slope_coefs[0] * wave + slope_coefs[1]

            # Correct the continuum
            flux = flux / slope_spectrum
        else:
            None

        if resolution is not None:
            from scipy.ndimage.filters import gaussian_filter1d

            # resolution R = λ / Δλ => Δλ = λ / R
            inst_profile_sig = np.median(wave) / resolution
            # inst_profile_sig = inst_profile_FWHM / (2 * np.sqrt(2 * np.log(2)))
            flux = gaussian_filter1d(flux, inst_profile_sig)
        if cache == True:
            np.savez(
                "spectrum" + str(uuid.uuid4()) + ".npz",
                wavelength=wave,
                flux=flux,
                prop_dict=prop_dict,
            )
        else:
            pass
        super().__init__(wave, flux)