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

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)