Calculate Simple Cometary Magnitude

import numpy as np
import pandas as pd
import astropy.units as u
from sorcha_addons.activity.lsst_comet.model import Comet
from sbpy.activity import Afrho
import synphot

from sorcha.modules.PPCalculateSimpleCometaryMagnitude import PPCalculateSimpleCometaryMagnitude

/home/docs/checkouts/readthedocs.org/user_builds/sorcha-addons/envs/stable/lib/python3.10/site-packages/rebound/__init__.py:58: UserWarning: pkg_resources is deprecated as an API. See https://setuptools.pypa.io/en/latest/pkg_resources.html. The pkg_resources package is slated for removal as early as 2025-11-30. Refrain from using this package or pin to Setuptools<81.
  import pkg_resources

The lsstcomet code used by sorcha validates its results by comparing them to the coma magnitude calculated by sbpy. We will do the same.

First, calculating using sbpy:

g = {'rh': 2.0 * u.au, 'delta': 1.0 * u.au, 'phase': 0 * u.deg}
afrho = Afrho(100 * 2**-2, 'cm')
tab = np.loadtxt('lsst-total-r.dat').T
r = synphot.SpectralElement(synphot.Empirical1D, points=tab[0] * u.nm,
                            lookup_table=tab[1])
rap = 1 * u.arcsec
m0 = afrho.to_fluxd(r, rap, g, unit=u.ABmag).value

comet = Comet(R=1, afrho1=100, k=-2)
m_sbpy = comet.mag(g, 'r', rap=rap.value)

Now a test dataset must be created using the same values. TrailedSourceMag here is a placeholder value: the function calculates the total apparent magnitude of coma and nucleus, and thus needs the “nucleus” apparent magnitude.

test_dict = {'MJD': [2459215.5],
              'H_r': [7.3],
              'afrho1': [100],
              'k':[-2],
              'optFilter':'r',
              'seeingFwhmEff': [1],
              'trailedSourceMagTrue': 18}

test_data = pd.DataFrame(test_dict)

rho = 2.
delta = 1.
alpha = 0

Calculating coma apparent magnitude using the SSPP function and comparing:

from sorcha_addons.activity.lsst_comet.lsst_comet_activity import LSSTCometActivity
from sorcha.activity.activity_registration import update_activity_subclasses

update_activity_subclasses()

test_data = PPCalculateSimpleCometaryMagnitude(test_data.copy(), ['r'], rho, delta, alpha, 'lsst_comet')

test_data

	MJD	H_r	afrho1	k	optFilter	seeingFwhmEff	trailedSourceMagTrue	coma_magnitude
0	2459215.5	7.3	100	-2	r	1	17.578331	18.809232

m_sspp = test_data['coma_magnitude']

assert np.isclose(m_sspp, m_sbpy, atol=0.05)