Skip to content

WIP: Reflectance concepts #329

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
jianyangli wants to merge 7 commits into
base: main
Choose a base branch
from
Open

WIP: Reflectance concepts #329

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
7 commits
Select commit Hold shift + click to select a range
f35e481
add albedo unit
jianyangli May 1, 2022
7a6cf74
expand equivalencies for disk-integrated albedo
jianyangli May 2, 2022
2bf739e
update tests for `sbpy.unit`
jianyangli May 2, 2022
92e902f
fix `photometry` module
jianyangli May 3, 2022
301092e
PEP8 fixes
jianyangli May 3, 2022
666ce60
update documents
jianyangli May 3, 2022
ebf192b
Merge branch 'main' into update_reflectance_size_conversion
jianyangli Jun 28, 2022
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
4 changes: 2 additions & 2 deletions docs/sbpy/photometry.rst
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -65,9 +65,9 @@ Calculate geometric albedo, Bond albedo, and phase integral:
>>> solar_fluxd.set({'V': -26.77 * u.mag})
<ScienceState solar_fluxd: {'V': <Quantity -26.77 mag>}>
>>> print(m.geomalb) # doctest: +FLOAT_CMP
0.09557298727307795
0.09557298727307795 albedo
>>> print(m.bondalb) # doctest: +FLOAT_CMP
0.03482207291799989
0.03482207291799989 albedo
>>> print(m.phaseint) # doctest: +FLOAT_CMP
0.3643505755292945

Expand Down
26 changes: 14 additions & 12 deletions docs/sbpy/units.rst
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -78,7 +78,8 @@ phase angle can be calculated:
>>> import numpy as np
>>> from astropy import units as u
>>> from sbpy.calib import solar_fluxd, vega_fluxd
>>> from sbpy.units import reflectance, VEGAmag, spectral_density_vega
>>> from sbpy.units import (dimensionless_albedo, albedo_unit, VEGAmag,
... spectral_density_vega)
>>>
>>> solar_fluxd.set({'V': -26.775 * VEGAmag})
<ScienceState solar_fluxd: {'V': <Magnitude -26.775 mag(VEGA)>}>
Expand All @@ -87,39 +88,40 @@ phase angle can be calculated:
>>> mag = 3.4 * VEGAmag
>>> radius = 460 * u.km
>>> cross_sec = np.pi * (radius)**2
>>> ref = mag.to('1/sr', reflectance('V', cross_section=cross_sec))
>>> ref = mag.to(albedo_unit, dimensionless_albedo('V',
... cross_section=cross_sec))
>>> print('{0:.4f}'.format(ref))
0.0287 1 / sr
0.0900 albedo

`~sbpy.units.reflectance` works with `sbpy`'s spectral calibration system (see :ref:`sbpy-calib`):

>>> from sbpy.photometry import bandpass
>>> V = bandpass('Johnson V')
>>> ref = 0.0287 / u.sr
>>> cross_sec = mag.to('km2', reflectance(V, reflectance=ref))
>>> ref = 0.09 * albedo_unit
>>> cross_sec = mag.to('km2', dimensionless_albedo(V, albedo=ref))
>>> radius = np.sqrt(cross_sec / np.pi)
>>> print('{0:.2f}'.format(radius))
459.69 km
460.11 km

`~sbpy.units.reflectance` also supports conversion between a flux spectrum and a reflectance spectrum:

>>> wave = [1.046, 1.179, 1.384, 1.739, 2.416] * u.um
>>> flux = [1.636e-18, 1.157e-18, 8.523e-19, 5.262e-19, 1.9645e-19] \
... * u.Unit('W/(m2 um)')
>>> xsec = 0.0251 * u.km**2
>>> ref = flux.to('1/sr', reflectance(wave, cross_section=xsec))
>>> ref = flux.to(albedo_unit, dimensionless_albedo(wave, cross_section=xsec))
>>> print(ref) # doctest: +FLOAT_CMP
[0.0021763 0.00201223 0.0022041 0.00269637 0.00292785] 1 / sr
[0.00683705 0.00632161 0.00692438 0.0084709 0.00919811] albedo

`~sbpy.units.reflectance` also supports dimensionless logarithmic unit for
solar flux, which can be specified with `~sbpy.calib.solar_fluxd.set`:
`~sbpy.units.dimensionless_albedo` also supports dimensionless logarithmic unit
for solar flux, which can be specified with `~sbpy.calib.solar_fluxd.set`:

>>> with solar_fluxd.set({'V': -26.77 * u.mag}):
... mag = 3.4 * u.mag
... ref = mag.to('1/sr', reflectance('V', cross_section =
... ref = mag.to(albedo_unit, dimensionless_albedo('V', cross_section =
... np.pi*(460*u.km)**2))
>>> print(ref) # doctest: +FLOAT_CMP
0.028786262941247264 1 / sr
0.09043471218052651 albedo

Projected Sizes
---------------
Expand Down
7 changes: 4 additions & 3 deletions sbpy/photometry/bandpass.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -4,13 +4,14 @@

"""

import os
from astropy.utils.data import get_pkg_data_filename


__all__ = [
'bandpass'
]

import os
from astropy.utils.data import get_pkg_data_filename


def bandpass(name):
"""Retrieve bandpass transmission spectrum from sbpy.
Expand Down
50 changes: 24 additions & 26 deletions sbpy/photometry/core.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -6,9 +6,6 @@

"""

__all__ = ['DiskIntegratedPhaseFunc', 'LinearPhaseFunc', 'HG', 'HG12BaseClass',
'HG12', 'HG1G2', 'HG12_Pen16', 'NonmonotonicPhaseFunctionWarning']

from collections import OrderedDict
import warnings
import numpy as np
Expand All @@ -20,10 +17,14 @@
from ..data import (Phys, Obs, Ephem, dataclass_input,
quantity_to_dataclass)
from ..bib import cite
from ..units import reflectance
from ..units import dimensionless_albedo, albedo_unit
from ..exceptions import SbpyWarning


__all__ = ['DiskIntegratedPhaseFunc', 'LinearPhaseFunc', 'HG', 'HG12BaseClass',
'HG12', 'HG1G2', 'HG12_Pen16', 'NonmonotonicPhaseFunctionWarning']


class _spline(object):

"""Cubic spline class
Expand Down Expand Up @@ -131,9 +132,9 @@ class DiskIntegratedPhaseFunc(Fittable1DModel):
... geomalb = linear_phasefunc.geomalb
... phaseint = linear_phasefunc.phaseint
... bondalb = linear_phasefunc.bondalb
>>> print('Geometric albedo is {0:.3}'.format(geomalb))
>>> print('Geometric albedo is {0:.3}'.format(geomalb.value))
Geometric albedo is 0.0487
>>> print('Bond albedo is {0:.3}'.format(bondalb))
>>> print('Bond albedo is {0:.3}'.format(bondalb.value))
Bond albedo is 0.0179
>>> print('Phase integral is {0:.3}'.format(phaseint))
Phase integral is 0.367
Expand Down Expand Up @@ -235,10 +236,7 @@ def _check_unit(self):
@property
def geomalb(self):
"""Geometric albedo"""
alb = np.pi*self.to_ref(0.*u.rad)
if hasattr(alb, 'unit') and (alb.unit == 1/u.sr):
alb = alb*u.sr
return alb
return self.to_ref(0. * u.rad)

@property
def bondalb(self):
Expand Down Expand Up @@ -375,9 +373,9 @@ def to_phys(self):
<class 'sbpy.data.phys.Phys'>
>>> p.table.pprint(max_width=-1) # doctest: +SKIP
targetname diameter H G pv A
km mag
----------------- -------- ---- ---- ------------------- --------------------
1 Ceres (A801 AA) 939.4 3.31 0.12 0.07624470768627523 0.027779803126557152
km mag albedo albedo
----------------- -------- ---- ---- ------------------- -------------------
1 Ceres (A801 AA) 939.4 3.31 0.12 0.09423445077857852 0.03433437637586201
""" # noqa: E501
cols = {}
if (self.meta is not None) and ('targetname' in self.meta.keys()):
Expand Down Expand Up @@ -611,10 +609,9 @@ def to_mag(self, eph, unit=None, append_results=False, **kwargs):
' size of object is unknown.')
if self.wfb is None:
raise ValueError('Wavelength/Frequency/Band is unknown.')
out = out.to(
unit,
reflectance(self.wfb, cross_section=np.pi * self.radius**2)
)
xsec = np.pi * self.radius**2
out = out.to(unit,
dimensionless_albedo(self.wfb, cross_section=xsec))
dist_corr = self._distance_module(eph)
dist_corr = u.Quantity(dist_corr).to(u.mag, u.logarithmic())
out = out - dist_corr
Expand Down Expand Up @@ -704,10 +701,11 @@ def to_ref(self, eph, normalized=None, append_results=False, **kwargs):
' phase function can be calculated.')
if self.wfb is None:
raise ValueError('Wavelength/Frequency/Band is unknown.')
xsec = np.pi * self.radius**2
out = out.to(
'1/sr',
reflectance(self.wfb, cross_section=np.pi*self.radius**2)
)
albedo_unit,
dimensionless_albedo(self.wfb, cross_section=xsec)
)
else:
out = out - norm
out = out.to('', u.logarithmic())
Expand Down Expand Up @@ -773,9 +771,9 @@ class LinearPhaseFunc(DiskIntegratedPhaseFunc):
... geomalb = linear_phasefunc.geomalb
... phaseint = linear_phasefunc.phaseint
... bondalb = linear_phasefunc.bondalb
>>> print('Geometric albedo is {0:.3}'.format(geomalb))
>>> print('Geometric albedo is {0:.3}'.format(geomalb.value))
Geometric albedo is 0.0487
>>> print('Bond albedo is {0:.3}'.format(bondalb))
>>> print('Bond albedo is {0:.3}'.format(bondalb.value))
Bond albedo is 0.0179
>>> print('Phase integral is {0:.3}'.format(phaseint))
Phase integral is 0.367
Expand Down Expand Up @@ -817,7 +815,7 @@ class HG(DiskIntegratedPhaseFunc):
>>> from sbpy.photometry import HG
>>> ceres = HG(3.34 * u.mag, 0.12, radius = 480 * u.km, wfb = 'V')
>>> with solar_fluxd.set({'V': -26.77 * u.mag}):
... print('geometric albedo = {0:.4f}'.format(ceres.geomalb))
... print('geometric albedo = {0:.4f}'.format(ceres.geomalb.value))
... print('phase integral = {0:.4f}'.format(ceres.phaseint))
geometric albedo = 0.0878
phase integral = 0.3644
Expand Down Expand Up @@ -988,7 +986,7 @@ class HG1G2(HG12BaseClass):
>>> 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('geometric albedo = {0:.4f}'.format(themis.geomalb.value))
... print('phase integral = {0:.4f}'.format(themis.phaseint))
geometric albedo = 0.0656
phase integral = 0.3742
Expand Down Expand Up @@ -1082,7 +1080,7 @@ class HG12(HG12BaseClass):
>>> from sbpy.photometry import HG12
>>> themis = HG12(7.121 * u.mag, 0.68, radius = 100 * u.km, wfb = 'V')
>>> with solar_fluxd.set({'V': -26.77 * u.mag}):
... print('geometric albedo = {0:.4f}'.format(themis.geomalb))
... print('geometric albedo = {0:.4f}'.format(themis.geomalb.value))
... print('phase integral = {0:.4f}'.format(themis.phaseint))
geometric albedo = 0.0622
phase integral = 0.3949
Expand Down Expand Up @@ -1180,7 +1178,7 @@ class HG12_Pen16(HG12):
>>> themis = HG12_Pen16(7.121 * u.mag, 0.68, radius = 100 * u.km,
... wfb = 'V')
>>> with solar_fluxd.set({'V': -26.77 * u.mag}):
... print('geometric albedo = {0:.4f}'.format(themis.geomalb))
... print('geometric albedo = {0:.4f}'.format(themis.geomalb.value))
... print('phase integral = {0:.4f}'.format(themis.phaseint))
geometric albedo = 0.0622
phase integral = 0.3804
Expand Down
Loading