feature request: add sHG1G2 model accounting for aspect angle variation #426
Description
This is a request for
- a new feature
- an enhancement to existing sbpy functionality
- somethings else: [explain here]
The requested changes will be implemented by
- [50%] me
- [50%] the sbpy developers
High-level concept
Phase functions in sbpy do not account for changing aspect angle, but a simple analytical model can account for it (Carry+2024).
Explain the relevance to sbpy
This is 100% related to small bodies :-)
We already use this model in the LSST broker FINK, and we have released a stand-alone python package that handles it: PHUNK.
Advantage of adding it to sbpy:
- it will be in sbpy :-) so useful for the community
- it can then be used in SORCHA (LSST survey simulator)
Proposal details
it should be a new class of phase function, just like the other HG HG12 HG1G2 etc, in
sbpy/photometry/iau.py
main difference though: the new model does not only depend on the phase angle but also on the equatorlal coordinates (RA,Dec) of the target. So the input observable are 3 (phase,RA,Dec) and not longer one (phase) /
Example (pseudo-)code
I have tried to code it in sbpy/photometry/iau.py, but I have issues with the fact the new model uses 3 inputs (phase,RA,Dec) and not one (phase).
I coded the new class in iau.py, I'm copying the "code" here, but it generates errors too deep in the dependencies of sbpy/astropy for me to be able to sort it out, sorry
(in particular in astropy/modeling/core.py:
Too many input arguments - expected 1, got 3
)
Here comes the proposed additional class:
Thanks!
class sHG1G2(HG12BaseClass):
"""sHG1G2 photometric phase model (Carry et al. 2024)
Examples
--------
>>> # Define the phase function for Themis with
>>> # H = 7.063, G1 = 0.62, G2 = 0.14
>>>
>>> import astropy.units as u
>>> from sbpy.calib import solar_fluxd
>>> from sbpy.photometry import HG1G2
>>> themis = HG1G2(7.063 * u.mag, 0.62, 0.14, radius = 100 * u.km,
... wfb = 'V')
>>> with solar_fluxd.set({'V': -26.77 * u.mag}):
... print('geometric albedo = {0:.4f}'.format(themis.geomalb))
... print('phase integral = {0:.4f}'.format(themis.phaseint))
geometric albedo = 0.0656
phase integral = 0.3742
"""
H = Parameter(description='H parameter', default=8)
G1 = Parameter(description='G1 parameter', default=0.2,
bounds=_hg1g2_g1_bounds)
G2 = Parameter(description='G2 parameter', default=0.2,
bounds=_hg1g2_g2_bounds)
RA0 = Parameter(description='Spin right ascension', bounds=_shg1g2_ra0_bounds)
Dec0 = Parameter(description='Spin declination', bounds=_shg1g2_dec0_bounds)
R = Parameter(description='Oblateness parameter', bounds=_shg1g2_R_bounds)
def __init__(self, *args, **kwargs):
print(*args)
ineqcons = kwargs.pop('ineqcons', [])
ineqcons.append(lambda x, *args: 1 - x[1] - x[2])
super().__init__(*args, ineqcons=ineqcons, **kwargs)
@G1.validator
def G1(self, value):
_check_bounds(value, _hg1g2_g1_bounds, InvalidPhaseFunctionWarning,
'G1 results in an invalid phase function.')
if (value + self.G2 > 1).any():
warnings.warn('G1, G2 combination results in an invalid '
'phase function.')
@G2.validator
def G2(self, value):
_check_bounds(value, _hg1g2_g2_bounds, InvalidPhaseFunctionWarning,
'G2 results in an invalid phase function.')
if (value + self.G1 > 1).any():
warnings.warn('G1, G2 combination results in an invalid '
'phase function.')
@RA0.validator
def RA0(self, value):
_check_bounds(value, _shg1g2_ra0_bounds, InvalidPhaseFunctionWarning,
'Out of bound spin right ascension (RA0).')
@Dec0.validator
def Dec0(self, value):
_check_bounds(value, _shg1g2_dec0_bounds, InvalidPhaseFunctionWarning,
'Out of bound spin declination (Dec0).')
@R.validator
def R(self, value):
_check_bounds(value, _shg1g2_R_bounds, InvalidPhaseFunctionWarning,
'Non-physical oblateness (R).')
@property
def _G1(self):
return self.G1.value
@property
def _G2(self):
return self.G2.value
@property
def _RA0(self):
return self.RA0.value
@property
def _Dec0(self):
return self.Dec0.value
@property
def _R(self):
return self.R.value
def cos_aspect_angle(ra, dec, ra0, dec0):
"""Compute the cosine of the aspect angle
This angle is computed from the coordinates of the target and
the coordinates of its pole.
See Eq 12.4 "Introduction to Ephemerides and Astronomical Phenomena", IMCCE
ISBN-13: 978-2759824144
Parameters
----------
ra: float
Right ascension of the target in radians.
dec: float
Declination of the target in radians.
ra0: float
Right ascension of the pole in radians.
dec0: float
Declination of the pole in radians.
Returns
-------
float: The cosine of the aspect angle
"""
ra = u.Quantity(ra, 'rad').to_value('rad')
dec = u.Quantity(dec, 'rad').to_value('rad')
ra0 = u.Quantity(ra0, 'rad').to_value('rad')
dec0 = u.Quantity(dec0, 'rad').to_value('rad')
return np.sin(dec) * np.sin(dec0) + np.cos(dec) * np.cos(dec0) * np.cos(ra - ra0)
@staticmethod
def evaluate(ph, ra, dec, h, g1, g2, RA0, Dec0, R):
ph = u.Quantity(ph, 'rad').to_value('rad')
ra = u.Quantity(ra, 'rad').to_value('rad')
dec = u.Quantity(dec, 'rad').to_value('rad')
# Phase part
func = g1*HG1G2._phi1(ph)+g2*HG1G2._phi2(ph)+(1-g1-g2)*HG1G2._phi3(ph)
if isinstance(func, u.Quantity):
func = func.value
func = -2.5 * np.log10(func)
if isinstance(h, u.Quantity):
func = func * h.unit
# Geometry part
geo = cos_aspect_angle(ra, dec, RA0, Dec0)
func2 = 1 - (1 - R) * np.abs(geo)
func2 = 2.5 * np.log10(func2)
if isinstance(h, u.Quantity):
func2 = func2 * h.unit
return h + func + func2
@staticmethod
def fit_deriv(ph, h, g1, g2):
if hasattr(ph, '__iter__'):
ddh = np.ones_like(ph)
else:
ddh = 1.
phi1 = HG1G2._phi1(ph)
phi2 = HG1G2._phi2(ph)
phi3 = HG1G2._phi3(ph)
dom = (g1*phi1+g2*phi2+(1-g1-g2)*phi3)
ddg1 = 1.085736205*(phi3-phi1)/dom
ddg2 = 1.085736205*(phi3-phi2)/dom
return [ddh, ddg1, ddg2]
def _parameter_units_for_data_units(self, inputs_unit, outputs_unit):
return OrderedDict([('H', outputs_unit['y']),
('G1', u.dimensionless_unscaled),
('G2', u.dimensionless_unscaled),
('RA0', u.rad),
('Dec0', u.rad),
('R', u.dimensionless_unscaled)])