Production Rate

After fitting a fluorescence model, the column density \(\log_{10} N\) (\(N\) in molecules cm^-2) can be converted to a molecular production rate \(Q\) (molecules s^-1) using one of two approaches:

• Haser model — models a spherically expanding coma with parent and daughter photodissociation scale lengths.

• ρ^-1 profile — assumes a pure radial-outflow column-density profile \(N(\rho) \propto \rho^{-1}\) and uses a pre-computed aperture integral.

The examples below build a short MCMC run first so that posterior samples are available. To skip the sampling step, pass use_samples=False to use the median logN directly (no error estimate is returned in that case); if you use this approach, there is no need to fit.

Sections

• Haser model

• ρ-1 profile via aperture integral

• Multi-isotopologue production rates

• Seeing-corrected production rates

• References

---

Let’s prepare a fitted model from which we will compute production rates.

from cometspec.helper import open_kurucz_irradiance, open_hall_anderson_irradiance
from cometspec import FluorescenceModel
from matplotlib import pyplot as plt
import numpy as np
import pandas as pd

irradiance = open_kurucz_irradiance()

model = FluorescenceModel(
    window=(3870, 3885),
    pumping=irradiance,
    isotopologues="12C14N",
    systems=['BX', 'AX_dv1'],
    sigma=0.05,
    logQ=-1,
    logN=14.0,
    T=200,
    v_kms=0,
    lsf_method="Gauss",
    wave_col="WAVE",
    flux_col="FLUX",
    error_col="ERR",
    continuum_col="CONTINUUM",
)

# add noise
data_flux = model.best_model + np.random.normal(0, np.max(model.best_model)*0.05,
                                                    size=model.model_wave.shape)
data = pd.DataFrame({'WAVE': model.model_wave,
                    'FLUX': data_flux,
                    'ERR':np.max(model.best_model)*0.05,
                    'CONTINUUM': np.zeros_like(model.model_wave)})

results = model.fit_mcmc(
                data=data,
                nwalkers=10,
                nsteps=500,
                priors={
                    "logN":  (12.0, 16.0),
                    "T":     (50.0, 300.0),
                    "v_kms": (-2.0, 2.0),
                },
                make_plots=False,
                )

Haser model

compute_production_rate() uses the sbpy Haser implementation. You must supply the geocentric distance, aperture geometry, CN scale lengths, and the outflow velocity.

Typical CN scale lengths at 1 AU are: parent ~ 1.3 × 10⁴ km, daughter ~ 2.1 × 10⁵ km [1]. Both scale as \(r_h^2\). The outflow velocity can be estimated as \(v = 0.85 r_h^{-0.5}\) km s^-1 [2], where \(r_h\) is the heliocentric distance in AU.

rh = 1.0   # your heliocentric distance [AU]

logQ, logQ_err = model.compute_production_rate(
    delta_au=1.0, # your geocentric distance [AU]
    aperture={"type": "circular", "radius_arcsec": 5.0},
    parent_length_km=1.3e4 * rh**2,
    daughter_length_km=2.1e5 * rh**2,
    v_outflow_km_s=0.85 * rh**(-0.5),
)
print(f"log Q(CN) = {logQ:.2f} ± {logQ_err:.2f}  [s^-1]")

Rectangular (slit) apertures are also supported:

logQ_slit, logQ_slit_err = model.compute_production_rate(
    delta_au=1.0,
    aperture={"type": "rectangular", "width_arcsec": 2.0, "length_arcsec": 10.0},
    parent_length_km=1.3e4 * rh**2,
    daughter_length_km=2.1e5 * rh**2,
    v_outflow_km_s=0.85 * rh**(-0.5),
)
print(f"log Q(CN, slit) = {logQ_slit:.2f} ± {logQ_slit_err:.2f}")

The production rate and its error are stored in the attributes model.q and model.q_err after calling compute_production_rate.

---

ρ^-1 profile via aperture integral

For a \(N(\rho) \propto \rho^{-1}\) profile where \(\rho\) is the radial distance from the nucleus, first compute the geometric aperture integral G with compute_aperture_integral(), then pass it to compute_production_rate_from_profile().

G = model.compute_aperture_integral(
    aperture={"type": "circular", "radius_arcsec": 5.0},
    delta_au=1.0,
)
print("Aperture integral G:", G)

logQ_prof, logQ_prof_err = model.compute_production_rate_from_profile(
    G=G,
    delta_au=1.0,
    aperture={"type": "circular", "radius_arcsec": 5.0},
    v_outflow_km_s=0.85 * rh**(-0.5),
)
print(f"log Q (rho^-1 profile) = {logQ_prof:.2f} ± {logQ_prof_err:.2f}")

---

Multi-isotopologue production rates

When multiple isotopologues are fitted, compute_production_rate() returns a dict keyed by isotopologue name. You can also use the attribute model.q to access the production rates and errors for each isotopologue after calling compute_production_rate.

# Two-iso fit
model = FluorescenceModel(
    pumping=irradiance,
    isotopologues=["12C14N", "13C14N"], systems="BX",
    logN_by_iso={"12C14N": 14, "13C14N": 13.0},
    T=100.0, sigma=0.02, logQ=-1, lsf_method="Gauss",
    wave_col="WAVE", flux_col="FLUX", error_col="ERR",
    continuum_col="CONTINUUM",
)
data_flux = model.best_model + np.random.normal(0, np.max(model.best_model)*0.05,
                                                    size=model.model_wave.shape)
data = pd.DataFrame({'WAVE': model.model_wave,
                    'FLUX': data_flux,
                    'ERR':np.max(model.best_model)*0.05,
                    'CONTINUUM': np.zeros_like(model.model_wave)})


model.fit_mcmc(
    data=data,
    nwalkers=20, nsteps=300,
    priors={
        "logN_12C14N": (12.0, 16.0),
        "logN_13C14N": (11.0, 14.0),
        "T":           (50.0, 300.0),
    },
    make_plots=False, progress=False,
)

rates = model.compute_production_rate(
    delta_au=1.0,
    aperture={"type": "circular", "radius_arcsec": 5.0},
    parent_length_km=1.3e4 * rh**2,
    daughter_length_km=2.1e5 * rh**2,
    v_outflow_km_s=0.85 * rh**(-0.5),
)
for iso, (logQ_i, logQ_i_err) in rates.items():
    print(f"log Q({iso}) = {logQ_i:.2f} ± {logQ_i_err:.2f}")

Note

If you have multiple isotopologues and want to use use_samples=False, you need to have both present in logN_by_iso. On the other hand, if multiple isotopologues are present and you want to use use_samples=True, you will need to fit both isotopologues; otherwise the code raises an error. However, you can work around this by adding the missing isotopologue later to logN_by_iso, appending the missing logN_<isotopologue> at the end of param_keys, and adding an array of the repeated value to samples_pruned. An example of this can be found in the Jupyter notebook WorkFlow.ipynb of the GitHub repository (link).

---

Seeing-corrected production rates

add_slit_loss_error() inflates the column-density uncertainty to account for seeing-induced slit losses. Call it on the fitted model after compute_production_rate. This will also compute the slit-loss error of logN.

print('Before adding slit-loss error:')
for iso in model.q:
    print(f"log Q({iso}) = {model.q[iso]:.3f} ± {model.q_err[iso]:.3f}")
for iso in model.logN_by_iso:
    print(f"log N({iso}) = {model.logN_by_iso[iso]:.2f} ± {model.logN_err_by_iso[iso][0]:.2f}") # for N you have up and low error, we just show 1

model.add_slit_loss_error(
lambda_nm=388.3, # the wavelength of the fitted band or line [nm]
aperture={"type": "circular", "radius_arcsec": .5},
            # small aperture to illustrate the effect
eps_min_arcsec_500=0.7,   # best-case seeing at 500 nm [arcsec]
eps_max_arcsec_500=1.5,   # worst-case seeing at 500 nm [arcsec]
zmin_deg=20, # minimum zenith angle of observations [degrees]
zmax_deg=60, # maximum zenith angle of observations [degrees]
)

print('After adding slit-loss error:')
for iso in model.q:
    print(f"log Q({iso}) = {model.q[iso]:.3f} ± {model.q_err[iso]:.3f}")
for iso in model.logN_by_iso:
    print(f"log N({iso}) = {model.logN_by_iso[iso]:.2f} ± {model.logN_err_by_iso[iso][0]:.2f}") # for N you have up and low error, we just show 1

If only one isotopologue is fitted, you have logN and logN_err. We encourage you to read the documentation of add_slit_loss_error() for more details on the method and its assumptions. Since this method adds the initial error in quadrature with the slit-loss error, the slit-loss error is not recalculated if add_slit_loss_error is called multiple times, in order to avoid accumulating it. If you want to recalculate it, you will need to either 1) fit again, or 2) reset the attributes q_seeing_corrected = False and logN_seeing_corrected = False and run the production rate calculation again (if the latter is not applied, the error will start accumulating repeatedly).

References

[1]

A’Hearn, M. F., Millis, R. C., Schleicher, D. O., Osip, D. J., & Birch, P. V. 1995, Icarus, 118, 223 (link)

[2]

Cochran, A. L. & Schleicher, D. G. 1993, Icarus, 105, 235 (link)