PyMolFit

PyMolFit is an independent pure-Python package for modelling and correcting telluric absorption in astronomical spectra. It is based on the concepts, workflow, and radiative-transfer approach of the original ESO Molecfit, with a Python API for both spectrum files and wavelength-flux arrays.

PyMolFit does not import or run Molecfit and is not affiliated with or endorsed by ESO.

Installation

python -m pip install pymolfit

Install plotting support for the tutorial notebooks:

python -m pip install "pymolfit[interactive]"

PyMolFit supports Python 3.10 and newer.

Correct A Spectrum File

The main API is correct():

from pymolfit import correct

result = correct(
    input_path="spectrum.fits",
    wavelength_medium="air",
)

wavelength = result.corrected.wavelength
corrected_flux = result.corrected.flux

The same function accepts a previously loaded spectrum:

from pymolfit import correct, load_spectrum

spectrum = load_spectrum("spectrum.fits", wavelength_medium="air")
result = correct(spectrum=spectrum)

Use wavelength_medium="vacuum" for vacuum wavelengths. The argument may be omitted only when the FITS metadata declares the wavelength medium unambiguously.

For table-based echelle products, PyMolFit uses common order/detector columns to preserve real physical groups and common quality columns such as QUAL or DQ to exclude flagged pixels. Wide orders may be divided into smaller radiative-transfer chunks for memory control, but those chunks continue to share one smooth wavelength-alignment model for their physical order.

Results remain in memory unless an output is requested:

from pymolfit import load_fit_product, save_corrected_txt, save_fit_product

save_corrected_txt(result, "corrected_spectrum.txt")
save_fit_product(result, "fit_product.ecsv")
reloaded = load_fit_product("fit_product.ecsv")

The text file contains a compact corrected spectrum. The ECSV product also contains the fitted transmission, model, masks, metadata, and provenance. The printed fit summary reports the resolved atmosphere and line data, numerical chunks versus physical groups, fitted parameters, masked-pixel counts, and residual line-alignment diagnostics.

Correct Wavelength And Flux Arrays

Arrays do not have a FITS header, so observing metadata is supplied with an Observation:

from pymolfit import Observation, correct

observation = Observation(
    time="2025-09-28T07:44:51",
    latitude_deg=-24.627,
    longitude_deg=-70.404,
    altitude_m=2635.0,
    airmass=1.18,
    resolving_power=140_000,
    wavelength_frame="observatory",
    instrument="ESPRESSO",
)

result = correct(
    wavelength=wavelength,
    flux=flux,
    wavelength_unit="angstrom",
    wavelength_medium="air",
    observation=observation,
)

wavelength_medium states whether the numbers are air or vacuum wavelengths. wavelength_frame separately states their velocity reference frame. Supported frames are observatory, barycentric, and heliocentric.

Select Fitting Regions Interactively

PyMolFit can save fit and exclusion windows without manually copying wavelength values from a plot:

from pymolfit import TheoreticalSpectrum, load_spectrum, select_telluric_regions

stellar_template = TheoreticalSpectrum(
    path="stellar_model.dat",
    radial_velocity_kms=21.3,
    vsini_kms=124.0,
)

spectrum = load_spectrum(
    "spectrum.fits",
    wavelength_medium="air",
)

selector = select_telluric_regions(
    spectrum,
    theoretical_spectrum=stellar_template,
    output_path="telluric_regions.ecsv",
)

Candidate telluric transitions from the AER catalogue are marked automatically and colored by molecule. To create an initial selection automatically, enter a line count such as 100 and press Automatic. PyMolFit proposes fit windows around the transitions with the strongest expected atmospheric absorption, using both AER line intensity and the representative atmospheric column of each molecule. Automatic windows extend 12 sampled detector pixels to either side of each line by default, retaining the line core, wings, and nearby continuum needed to constrain alignment and instrumental broadening. It skips lines in detector/order gaps and merges overlapping windows. These proposals can be edited like any manual selection.

For manual editing, first zoom or pan to the desired area. Choose Fit or Exclude, enable the Draw regions checkbox, and drag rectangles around the lines; only their horizontal wavelength limits are stored. Drawing stays active until you clear the checkbox. Every stored region is numbered on the plot and listed in the side panel. Edit the filename field if needed, then press Save All once to write the complete collection. In a Jupyter notebook, run %matplotlib widget before opening the selector.

When a theoretical spectrum is supplied, automatically identified stellar features appear immediately as red exclusion regions. To edit the physical template parameters in the selector, pass enable_theoretical_controls=True. The optional panel then exposes radial velocity, v sin(i), resolving power, mask depth, mask padding, limb darkening, continuum window, and residual velocity search.

On later runs, the same output_path is detected and loaded automatically, so the selector window is skipped. Pass reuse_existing=False to reopen the saved regions for editing.

The ECSV file records both region types and their wavelength unit and air/vacuum medium. Apply it directly without transcribing any endpoints:

from pymolfit import correct

result = correct(
    input_path="spectrum.fits",
    wavelength_medium="air",
    region_file="telluric_regions.ecsv",
)

region_file cannot be combined with explicit fit_ranges or exclude_ranges.

Protect Stellar Features With A Theoretical Spectrum

An optional theoretical stellar spectrum can identify astrophysical lines that must not bias the telluric fit. The recommended workflow is to pass it to the region selector as shown above, producing one file containing telluric fit intervals and stellar exclusions. It can also be supplied directly to correct() when an interactive selector is not needed:

from pymolfit import TheoreticalSpectrum, correct

stellar_template = TheoreticalSpectrum(
    path="stellar_model.dat",
    radial_velocity_kms=21.3,
    vsini_kms=124.0,
    mask_padding_kms="auto",
)

result = correct(
    input_path="spectrum.fits",
    theoretical_spectrum=stellar_template,
    stellar_mask_path="stellar_exclusions.ecsv",
)

The input is a two-column ASCII wavelength/flux table. Files from the SVO Theoretical Spectra server are supported directly. PyMolFit reads their coordinate metadata, normalizes physical flux, applies the supplied radial velocity and projected rotation, and detects stellar features automatically. Those features become hard exclusions from atmospheric parameter estimation by default. Set confidence_weighted_masking=True on TheoreticalSpectrum to opt into continuous confidence weighting instead. PyMolFit uses observation metadata for frame and instrumental broadening, and checks a small residual alignment. The generated stellar model affects parameter estimation only. It does not replace the observed flux, and atmospheric transmission is still evaluated inside downweighted or excluded stellar regions. Automatic mask padding scales with the rotational and instrumental broadening so the fitted mask also protects broad stellar-line wings. A numeric mask_padding_kms can override that width.

stellar_mask_path is optional. When supplied, it saves the generated exclusions in the input spectrum's original wavelength unit and air/vacuum medium for inspection.

The selector displays the automatically detected exclusions when a template is supplied. Its editable theoretical-parameter panel is disabled by default; pass enable_theoretical_controls=True to the selector to show it.

For blended stellar and telluric lines, correct() also provides an opt-in joint forward model:

result = correct(
    input_path="spectrum.fits",
    theoretical_spectrum=stellar_template,
    joint_stellar_model=True,
)

This evaluates continuum * LSF(stellar * atmosphere) rather than using the stellar template only as fit weights. The correction divides by the stellar-aware atmospheric effect, preserving the convolved stellar spectrum. The mode is disabled by default and is available only through correct().

Tutorials

The runnable notebooks in tutorials/ cover:

1. Correcting part of a spectrum
2. Correcting a full spectrum
3. Expert parameters
4. Correcting wavelength-flux arrays
5. Selecting telluric fit regions interactively

Start with Tutorial 1 for a short example or Tutorial 2 for a complete one-dimensional echelle spectrum.

Documentation

• Software guide and public API
• Troubleshooting
• Physics parity audit
• Changelog
• Contributing

Molecular And Atmospheric Data

The standard workflow automatically obtains the versioned AER molecular line catalogue, verifies its checksum, and caches it under ~/.cache/pymolfit/aer. It is not bundled in the Python wheel, and the normal workflow does not require a HITRAN API key.

Users may instead provide their own HITRAN .par files or other supported line and continuum data.

Scientific Basis And References

PyMolFit builds on the methods and scientific data developed by several projects. Publications using PyMolFit should cite the relevant upstream resources:

• ESO Molecfit: Smette et al. 2015, A&A 576, A77 and Kausch et al. 2015, A&A 576, A78
• AER LBLRTM line-by-line radiative-transfer methods and data
• AER MT_CKD atmospheric continuum model
• HITRAN molecular spectroscopic data
• MIPAS atmospheric reference profiles
• NOAA GDAS meteorological profiles

PyMolFit uses NumPy, SciPy, and Astropy.

Development

git clone https://github.com/nikolasvalsamidis97/PyMolFit.git
cd PyMolFit
python -m pip install -e ".[dev,plot]"
python -m pytest

Report problems through the GitHub issue tracker.

Licensing And Third-Party Data

PyMolFit's original Python code is available under the BSD 3-Clause License. Scientific datasets and derived coefficient tables retain their upstream terms and attribution requirements. See THIRD_PARTY_NOTICES.md for the included and automatically downloaded data.