GUIBRUSHR.core.functions.spectroscopy module

Spectroscopy utility functions for GUIBRUSHR atmospheric retrieval.

This module contains functions for spectroscopic data processing extracted from general_functions.py:

  • kernel_solid_body_rotation: Compute a rotational broadening kernel.

  • convolve_solid_body_rotation: Convolve a spectrum with the rotational kernel.

  • convolve_resolution: Degrade a spectrum to instrumental resolution with Doppler shift.

  • calculate_mass_molecule: Calculate molecular mass from elemental composition.

  • generate_chemcat: Generate a chemical catalog for atmospheric modeling.

  • skewed_gaussian: Evaluate a skewed Gaussian (helper for estimate_continuum).

  • estimate_continuum: Estimate and remove the continuum from a spectrum.

  • compute_adjusted_abundance: Compute abundance accounting for dissociation effects.

  • spectral_convolution_ptr: Convolve a spectrum using petitRADTRANS.

GUIBRUSHR.core.functions.spectroscopy.kernel_solid_body_rotation(omegad: float, f_rot_arr: ndarray, Teq: float, mu: float, gs: float, Rpj: float, rad_mode: str, index_nights_total: int, rv_sampling: float = 0.1) → ndarray[source]

Compute a rotational broadening kernel for spectroscopic observations.

A planet rotating as a rigid body with angular frequency omega broadens spectral lines because different surface elements of the visible disk move at different line-of-sight velocities: an element at projected distance x from the rotation axis contributes flux at a radial velocity offset v = omega * x. The rotational kernel K(v) is the disk-integrated weight of all surface elements emitting at velocity v, and convolving a model spectrum with K(v) reproduces the line broadening induced by planetary rotation. The overall width of the kernel is set by the equatorial rotational velocity v_eq = omega * Rp.

The rotation rate omega is a free parameter of the retrieval and is not constrained to any particular value. For a tidally locked planet it equals 2*pi / P_orb, but the kernel is valid for any rotation rate.

The kernel geometry differs between emission and transmission spectroscopy:

  • In emission spectroscopy, only the planetary disk itself contributes. Each surface strip at projected distance x has a chord length 2*sqrt(Rp^2 - x^2), so K(v) is proportional to sqrt(1 - (v/v_eq)^2) for |v| < v_eq, and zero outside. This is the classical solid-body rotational kernel for a uniformly emitting disk: it peaks at v=0 (sub-observer point, longest chord) and drops to zero at v=±v_eq (planetary limb, chord=0).

  • In transmission spectroscopy, the relevant region is the atmospheric annulus between Rp and Rp + z, where z = 5*H is the effective atmospheric thickness and H = k_B * T_eq / (mu * m_H * g) is the scale height. Defining v_eq_atm = omega * (Rp + z), two velocity regions contribute differently:

    • v_eq <= |v| < v_eq_atm (annulus only): the line of sight passes only through the annulus, with weight sqrt((Rp+z)^2 - (v/omega)^2) / z.

    • |v| < v_eq (annulus above opaque disk): the line of sight passes through the full annulus minus the opaque planetary disk, with weight (sqrt((Rp+z)^2 - (v/omega)^2) - sqrt(Rp^2 - (v/omega)^2)) / z.

    An asymmetry parameter f_rot (stored and passed in linear units in f_rot_arr, converted from log10 by to_linear at extraction time) is applied to the receding hemisphere (v > 0) to allow for day/night asymmetries in the transmission signal. It is a free parameter of the retrieval.

In both modes the kernel is normalized to unit sum before being returned, so that convolution preserves the total flux of the spectrum.

Parameters:

  • omegad (float) – Rotational angular frequency of the planet, in inverse days (d^-1). The equatorial rotational velocity is v_eq = omega * Rp, where omega is omegad converted to s^-1. This is a free parameter of the retrieval and is not assumed to equal any particular value.

  • f_rot_arr (np.ndarray) – Array of linear rotational asymmetry scaling factors (converted from log10 by to_linear at extraction time), one per observing night. Used only in transmission mode: f_rot_arr[index_nights_total] is multiplied onto the kernel for v > 0 to model day/night hemisphere asymmetries. It is a free parameter of the retrieval.

  • Teq (float) – Equilibrium temperature of the planet, in Kelvin. Used to compute the atmospheric scale height H = k_B * T_eq / (mu * m_H * g), which defines the vertical extent z = 5*H of the transmission annulus. Only used in transmission mode.

  • mu (float) – Mean molecular weight of the atmosphere, in units of the hydrogen mass m_H. Used to compute the scale height. Only used in transmission mode.

  • gs (float) – Surface gravity of the planet, in cgs units (cm s^-2). Used to compute the scale height. Only used in transmission mode.

  • Rpj (float) – Planetary radius, in centimeters. Defines v_eq = omega * Rp (emission) and the inner boundary of the transmission annulus (transmission).

  • rad_mode (str) – Observing mode. Must be one of ConstantVariables.LIST_RAD_MODE: - LIST_RAD_MODE[0]: transmission spectroscopy (annulus kernel). - any other value: emission spectroscopy (solid disk kernel).

  • index_nights_total (int) – Index into f_rot_arr selecting the asymmetry factor for the current observing night. Only used in transmission mode.

  • rv_sampling (float, optional) – Velocity step of the radial velocity grid on which the kernel is computed, in km/s. Default is 0.1 km/s, fine enough to well-sample the kernel profile for any realistic v_eq.

Returns:

Normalized rotational broadening kernel, sampled on a velocity grid spanning -30 to +30 km/s with step rv_sampling. The kernel sums to unity by construction.

Return type:

np.ndarray

GUIBRUSHR.core.functions.spectroscopy.convolve_solid_body_rotation(wl: ndarray, model: ndarray, ker_rot: ndarray, rv_sampling: float = 0.1, n_el: int = 321) → Tuple[ndarray, ndarray][source]

Convolve a model spectrum with a rotational broadening kernel.

This function is the second convolution step in the model preparation pipeline, applied after convolve_resolution. It takes the rotational kernel produced by kernel_solid_body_rotation (defined on a fine velocity grid with step rv_sampling) and convolves it with the model spectrum.

The procedure mirrors that of convolve_resolution:

  1. The pixel velocity step of the model wavelength grid is computed as dv = c * dlambda / lambda, constant by construction for a log-uniform grid (constant resolving power R_model).

  2. The kernel is resampled from its fine velocity grid (step rv_sampling, default 0.1 km/s) onto a coarser grid whose spacing matches the model pixel velocity step, using spline interpolation. This step is mandatory because discrete convolution requires kernel and signal to share the same grid spacing.

  3. The resampled kernel is convolved with the model in ‘valid’ mode (scipy.signal.convolve), discarding the n_el//2 edge pixels on each side where the kernel does not fully overlap the signal, and the result is normalized by the kernel sum to ensure flux conservation.

Note: unlike convolve_resolution, this function does not apply a radial velocity shift. The shift is handled entirely by convolve_resolution via kernel translation. This function is responsible only for the rotational line broadening.

Parameters:

  • wl (np.ndarray) – 1D array of wavelength values for the model spectrum. Must be log-uniform (constant resolving power R_model), so that the pixel velocity step dv = c / R_model is constant across the grid.

  • model (np.ndarray) – 1D array of spectral flux values corresponding to wl. This is the spectrum already convolved with the instrumental LSF by convolve_resolution.

  • ker_rot (np.ndarray) – Rotational broadening kernel computed by kernel_solid_body_rotation, sampled on a velocity grid from -30 to +30 km/s with step rv_sampling. Must be normalized to unit sum.

  • rv_sampling (float, optional) – Velocity step of the input kernel grid, in km/s. Must match the rv_sampling used in kernel_solid_body_rotation. Default is 0.1 km/s.

  • n_el (int, optional) – Number of elements in the resampled kernel. Must be odd to ensure symmetry around zero and avoid introducing a spurious velocity shift. The half-extent of the resampled grid is (n_el - 1) / 2 pixels. Default is 321, giving a half-extent of 160 pixels.

Returns:

  • wavelength (np.ndarray) – Wavelength array of the convolved spectrum, shorter than the input by n_el // 2 = 160 points on each side due to ‘valid’ mode convolution.

  • rconv (np.ndarray) – Rotationally broadened flux array, flux-conserving by construction.

GUIBRUSHR.core.functions.spectroscopy.convolve_resolution(wavelength_model: ndarray, flux_model: ndarray, lsf_hwhm: float, lsf_grid_step: float = 0.1, kernel_size: int = 321) → Tuple[ndarray, ndarray][source]

Degrade a high-resolution model spectrum to instrumental resolution.

A real spectrograph cannot reproduce an infinitely narrow spectral line: even a perfectly monochromatic source produces a broadened profile on the detector, due to diffraction, optical aberrations, slit width, and pixel sampling. This broadening is described by the Line Spread Function (LSF). To compare a synthetic model with observed data, the model must be degraded to the same instrumental resolution, i.e. convolved with the LSF.

The LSF is approximated as a Gaussian whose width is expressed in velocity units. Working in velocity space is convenient because the LSF width is approximately constant in velocity across the full spectral range, whereas in wavelength units it varies with lambda. The Half Width at Half Maximum is HWHM = c / (2R), where R is the nominal resolving power of the instrument, and the corresponding Gaussian standard deviation is sigma = HWHM / sqrt(2*ln(2)). The factor sqrt(2*ln(2)) ~ 1.177 comes from setting the Gaussian equal to half its peak value and solving for v: G(v_1/2) = exp(-v_1/2^2 / 2*sigma^2) = 1/2 gives v_1/2 = sigma * sqrt(2*ln(2)), which is by definition the HWHM.

This function performs ONLY the resolution degradation. It does not apply any radial velocity shift: the kernel is built symmetric around zero velocity, so convolution preserves the spectral positions. All Doppler shifts (orbital, systemic, barycentric and the retrieval rv via rv_dynamic) are applied as a single velocity term on the data wavelength grid by the caller, through calculate_rv_planet_and_star (see cross_correlation.py). This keeps the radial velocity computed in one place rather than split between the kernel and vtot.

Parameters:

  • wavelength_model (np.ndarray) – 1D array of wavelength values for the model spectrum. Must be computed at constant resolving power R_model, so that the ratio dlambda/lambda = 1/R_model is constant by construction (log-uniform grid), yielding a constant velocity step dv = c / R_model per pixel.

  • flux_model (np.ndarray) – 1D array of flux values corresponding to wavelength_model.

  • lsf_hwhm (float) – Half Width at Half Maximum of the instrumental LSF, in km/s. Related to spectral resolution by: HWHM = c / (2R). Example: IGRINS (R ~ 45000) has HWHM ~ 3.3 km/s.

  • lsf_grid_step (float, optional) – Velocity step for the initial fine LSF kernel grid, in km/s. Must be small enough to well-sample the Gaussian peak. Default is 0.1 km/s.

  • kernel_size (int, optional) – Number of elements in the resampled convolution kernel. Must be odd to ensure the grid is symmetric around zero, which avoids introducing a spurious velocity shift during convolution. The half-extent of the resampled grid is (kernel_size - 1) / 2 pixels. Default is 321, giving a half-extent of 160 pixels.

Returns:

  • wavelength_out (np.ndarray) – Wavelength array of the convolved spectrum. Shorter than the input by kernel_size // 2 = 160 points on each side, because scipy.signal.convolve in ‘valid’ mode discards the edge pixels where the kernel does not fully overlap the signal (see Step 5).

  • flux_convolved (np.ndarray) – Convolved flux array, flux-conserving by construction.

GUIBRUSHR.core.functions.spectroscopy.calculate_mass_molecule(elements, periodic_table)[source]

Calculate the total molecular mass of a molecule from its elemental composition.

This function parses element strings that may contain mass numbers, atom counts, and ionization states, then calculates the total molecular mass using atomic masses from the periodic table.

Parameters:

  • elements – List of element strings with format like ‘H2’, ‘16O’, ‘C1+’, etc.

  • periodic_table – DataFrame containing periodic table data with columns ‘Symbol’, ‘NumberofProtons’, ‘NumberofNeutrons’

Returns:

Total molecular mass in atomic mass units (u)

Return type:

float

GUIBRUSHR.core.functions.spectroscopy.generate_chemcat(min_pressure: float, max_pressure: float, nlayersf: int, path_default: str | Path) → Tuple[ndarray, chemcat.Network, List[str]][source]

Generate a chemical catalog and retrieve molecular masses for atmospheric modeling.

This function creates a chemical network with specified pressure layers and temperature, and prepares molecular data for atmospheric modeling.

Parameters:

  • min_pressure (float) – Logarithm of minimum pressure (bar)

  • max_pressure (float) – Logarithm of maximum pressure (bar)

  • nlayersf (int) – Number of pressure layers

  • path_default (Union[str, Path]) – Base path to default files directory

Returns:

  • np.ndarray: Array of molecular masses for selected molecules

  • chemcat.Network: Chemical network object

  • List[str]: List of species names compatible with petitRADTRANS

Return type:

Tuple containing

GUIBRUSHR.core.functions.spectroscopy.skewed_gaussian(x, amp, center, width, skew)[source]

Calculate a skewed Gaussian distribution.

Parameters:

  • x – Input values

  • amp – Amplitude of the distribution

  • center – Center position

  • width – Width parameter

  • skew – Skewness parameter

Returns:

Values of the skewed Gaussian distribution

GUIBRUSHR.core.functions.spectroscopy.estimate_continuum(spectrum_new, num_windows=40, top_window=10, check_gaussian=False)[source]

Estimate and remove the continuum from a spectrum using windowed median fitting.

Parameters:

  • spectrum_new (np.ndarray) – Input spectrum array.

  • num_windows (int) – Number of windows to divide the spectrum into.

  • top_window (int) – Number of top values per window used for continuum estimation.

  • check_gaussian (bool) – If True, fit a skewed Gaussian instead of a polynomial.

Returns:

(continuum, model_after_processing) where continuum is the estimated

continuum and model_after_processing is the spectrum normalized by it.

Return type:

tuple

GUIBRUSHR.core.functions.spectroscopy.compute_adjusted_abundance(pressure: ndarray, temperature: ndarray, coeff, base_vmr) → ndarray[source]

Compute adjusted abundance accounting for dissociation effects.

This function calculates the adjusted abundance of a species by considering both the base volume mixing ratio and dissociation effects. The dissociation is modeled using a pressure and temperature dependent function.

Parameters:

  • pressure (np.ndarray) – Pressure values in bar

  • temperature (np.ndarray) – Temperature values in Kelvin

  • coeff – Coefficients [a, b, c] for dissociation calculation

  • base_vmr – Base volume mixing ratio from reference layer

Returns:

Adjusted abundance values

Return type:

np.ndarray

The calculation follows these steps: 1. Compute dissociation factor (Ad) using pressure and temperature 2. Convert to inverse dissociation factor (1/Ad) 3. Combine with inverse base abundance (1/A0) in quadrature 4. Return final adjusted abundance (A)

GUIBRUSHR.core.functions.spectroscopy.spectral_convolution_ptr(wavelength, spectrum, pixel_dv)[source]

Convolve a spectrum to instrumental resolution using petitRADTRANS.