diff --git a/Ska/Sun/data/pitch_roll.fits.gz b/Ska/Sun/data/pitch_roll.fits.gz new file mode 100644 index 0000000..f74bc93 Binary files /dev/null and b/Ska/Sun/data/pitch_roll.fits.gz differ diff --git a/Ska/Sun/sun.py b/Ska/Sun/sun.py index 468f05f..9d23efc 100755 --- a/Ska/Sun/sun.py +++ b/Ska/Sun/sun.py @@ -2,15 +2,46 @@ """ Utility for calculating sun position and pitch angle. """ +from pathlib import Path +import numpy as np +from math import cos, sin, acos, atan2, asin, pi, radians, degrees +from astropy.table import Table +from ska_helpers.utils import LazyVal from Quaternion import Quat from chandra_aca.transform import radec_to_eci from Chandra.Time import DateTime -from math import cos, sin, acos, atan2, asin, pi, radians, degrees -import numpy as np import Ska.quatutil +def load_roll_table(): + dat = Table.read(Path(__file__).parent / 'data' / 'pitch_roll.fits.gz') + + # Add a terminating row to the data such that for pitch at or greater + # than 180 the allowed roll deviation is defined as 0. + dat.add_row({'pitch': 180, 'rolldev': 0}) + return dat + + +ROLL_TABLE = LazyVal(load_roll_table) + + +def allowed_rolldev(pitch): + """Get allowed roll deviation (off-nominal roll) for the given ``pitch``. + :param pitch: Sun pitch angle (deg) + :returns: Roll deviation (deg) + """ + idx1 = np.searchsorted(ROLL_TABLE.val['pitch'], pitch, side='right') + idx0 = idx1 - 1 + idx_max = len(ROLL_TABLE.val) - 1 + idx0 = np.clip(idx0, 0, idx_max) + idx1 = np.clip(idx1, 0, idx_max) + val0 = ROLL_TABLE.val['rolldev'][idx0] + val1 = ROLL_TABLE.val['rolldev'][idx1] + out = np.minimum(val0, val1) # works even for a vector input for `pitch` + return out + + # The position() method is a modification of # http://idlastro.gsfc.nasa.gov/ftp/pro/astro/sunpos.pro # diff --git a/Ska/Sun/tests/__init__.py b/Ska/Sun/tests/__init__.py new file mode 100644 index 0000000..e69de29 diff --git a/Ska/Sun/tests/test_sun.py b/Ska/Sun/tests/test_sun.py index d11c0b7..4c420fe 100644 --- a/Ska/Sun/tests/test_sun.py +++ b/Ska/Sun/tests/test_sun.py @@ -3,7 +3,31 @@ import numpy as np from Quaternion import Quat from ..sun import (apply_sun_pitch_yaw, get_sun_pitch_yaw, nominal_roll, - off_nominal_roll, position, pitch as sun_pitch) + off_nominal_roll, position, allowed_rolldev, + pitch as sun_pitch) + + +def test_allowed_rolldev(): + + # Test array of pitchs and allowed roll dev + testarr = [[135, 13.979], + [138, 14.516], + [0, 0], + [40, 0], + [179.9, 18.748772], + [179.997, 0], + [180, 0], + [181, 0], + [85.49229, 13.677669], + [85.52, 18.756727], + [124.99, 18.748772], + [125, 17.0]] + for pitch, rolldev in testarr: + assert np.isclose(allowed_rolldev(pitch), rolldev) + + # Also test with pitch as vector + assert np.allclose(allowed_rolldev(np.array(testarr)[:, 0]), + np.array(testarr)[:, 1]) def test_position(): diff --git a/make_rolldev_file.py b/make_rolldev_file.py new file mode 100644 index 0000000..f523ef9 --- /dev/null +++ b/make_rolldev_file.py @@ -0,0 +1,20 @@ +# For the pitch/roll table. I used google sheets to save the xlsx file as +# a CSV (after adding headers: pitch, negrolldev, rolldev). +# For that CSV file, this little script checked for matching neg and +# positive roll values and then saved out a compressed file. + +# For future versions of the file, we should convert directly +# with something like +# df = pandas.read_excel('pitch_roll_2022.xlsx' +# ...: , names=['pitch', 'rolldev_minus', 'rolldev_plus']) +# dat = Table.from_pandas(df) + +import numpy as np +from astropy.table import Table +dat = Table.read('pitch_roll_2022.csv') + +assert np.allclose(np.abs(dat['negrolldev']), dat['rolldev'], rtol=0, atol=1e-08) + +dat['pitch'] = dat['pitch'].astype(np.float32) +dat['rolldev'] = dat['rolldev'].astype(np.float32) +dat['pitch', 'rolldev'].write('pitch_roll.fits.gz', overwrite=True) diff --git a/setup.py b/setup.py index 58a1095..d274b18 100755 --- a/setup.py +++ b/setup.py @@ -14,7 +14,7 @@ setup_requires=['setuptools_scm', 'setuptools_scm_git_archive'], zip_safe=False, packages=['Ska', 'Ska.Sun', 'Ska.Sun.tests'], - package_data={}, + package_data={'Ska.Sun': ['data/*fits.gz']}, tests_require=['pytest'], cmdclass=cmdclass, )