# poshist.py: GBM Position History class
#
# Authors: William Cleveland (USRA),
# Adam Goldstein (USRA) and
# Daniel Kocevski (NASA)
#
# Portions of the code are Copyright 2020 William Cleveland and
# Adam Goldstein, Universities Space Research Association
# All rights reserved.
#
# Written for the Fermi Gamma-ray Burst Monitor (Fermi-GBM)
#
# This program is free software: you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with this program. If not, see <https://www.gnu.org/licenses/>.
#
import astropy.io.fits as fits
import numpy as np
from numpy.lib.recfunctions import stack_arrays
from collections import OrderedDict
from scipy.interpolate import interp1d
from warnings import warn
from gbm import coords
from .data import DataFile
from gbm.detectors import Detector
[docs]class PosHist(DataFile):
"""Class for Fermi Position History
Attributes:
datatype (str): The datatype of the file
detector (str): The GBM detector the file is associated with
directory (str): The directory the file is located in
filename (str): The filename
full_path (str): The full path+filename
gti ([(float, float), ...]): A list of time intervals where Fermi is
outside the SAA
headers (dict): The headers for each extension of the poshist file
id (str): The GBM file ID
is_gbm_file (bool): True if the file is a valid GBM standard file,
False if it is not.
is_trigger (bool): True if the file is a GBM trigger file, False if not
sun_occulted ([(float, float), ...]): A list of time intervals where the
sun is occulted by the Earth
time_range (float, float): The time range of the poshist
"""
def __init__(self):
super(PosHist, self).__init__()
self._headers = OrderedDict()
self._data = None
self._times = None
self._eic_interp = None
self._quat_interp = None
self._lat_interp = None
self._lon_interp = None
self._alt_interp = None
self._vel_interp = None
self._angvel_interp = None
self._sun_interp = None
self._saa_interp = None
self._gti = None
self._sun_occulted = None
self._earth_radius_interp = None
self._geocenter_interp = None
self._detectors = [det.short_name for det in Detector]
@property
def headers(self):
return self._headers
@property
def time_range(self):
return (self._times[0], self._times[-1])
@property
def gti(self):
return self._gti
@property
def sun_occulted(self):
return self._sun_occulted
[docs] def get_eic(self, times):
"""Retrieve the position of Fermi in Earth inertial coordinates
Args:
times (float or np.array): Time(s) in MET
Returns:
np.array: A (3, `n`) array of coordinates at `n` times
"""
return self._eic_interp(times)
[docs] def get_quaternions(self, times):
"""Retrieve the Fermi attitude quaternions
Args:
times (float or np.array): Time(s) in MET
Returns:
np.array: A (4, `n`) array of quaternions at `n` times
"""
return self._quat_interp(times)
[docs] def get_latitude(self, times):
"""Retrieve the position of Fermi in Earth latitude
Args:
times (float or np.array): Time(s) in MET
Returns:
np.array: An array of latitudes at the requested times
"""
return self._lat_interp(times)
[docs] def get_longitude(self, times):
"""Retrieve the position of Fermi in Earth East longitude
Args:
times (float or np.array): Time(s) in MET
np.array: An array of East longitudes at the requested times
"""
return self._lon_interp(times)
[docs] def get_altitude(self, times):
"""Retrieve the altitude of Fermi in orbit
Args:
times (float or np.array): Time(s) in MET
Returns:
np.array: An array of altitudes at the requested times
"""
return self._alt_interp(times)
[docs] def get_velocity(self, times):
"""Retrieve the velocity of Fermi in Earth inertial coordinates
Args:
times (float or np.array): Time(s) in MET
Returns:
np.array: A (3, `n`) array of velocity components at `n` times
"""
return self._vel_interp(times)
[docs] def get_angular_velocity(self, times):
"""Retrieve the angular velocity of Fermi
Args:
times (float or np.array): Time(s) in MET
Returns
np.array: A (3, `n`) array of angular velocity components at `n` times
"""
return self._angvel_interp(times)
[docs] def get_sun_visibility(self, times):
"""Determine if the sun is visible (not occulted)
Args:
times (float or np.array): Time(s) in MET
Returns:
np.array(dtype=bool): A boolean array, where True indicates the \
sun is visible
"""
return self._sun_interp(times).astype(bool)
[docs] def get_saa_passage(self, times):
"""Determine if Fermi is in an SAA passage
Args:
times (float or np.array): Time(s) in MET
Returns:
np.array(dtype=bool): A boolean array, where True indicates Fermi \
is in SAA
"""
return self._saa_interp(times).astype(bool)
[docs] def get_earth_radius(self, times):
"""Retrieve the angular radius of the Earth as observed by Fermi
Args:
times (float or np.array): Time(s) in MET
Returns:
np.array: An array of the apparent angular radius of the Earth
"""
return self._earth_radius_interp(times)
[docs] def get_geocenter_radec(self, times):
"""Retrieve the Equatorial position of the Geocenter as observed by Fermi
Args:
times (float or np.array): Time(s) in MET
Returns:
(np.array, np.array): The RA and Dec of the geocenter
"""
eic = self.get_eic(times)
return coords.geocenter_in_radec(eic)
[docs] def get_mcilwain_l(self, times):
"""Retrieve the approximate McIlwain L value as determined by the
orbital position of Fermi
Args:
times (float or np.array): Time(s) in MET
Returns:
np.array: An array of McIlwain L values
"""
lat = self.get_latitude(times)
lon = self.get_longitude(times)
ml = coords.calc_mcilwain_l(lat, lon)
return ml
[docs] def to_fermi_frame(self, ra, dec, times):
"""Convert an equatorial position to a position in spacecraft coordinates
Args:
ra (float): The RA of a sky position
dec (float): The Dec of a sky position
times (float or np.array): Time(s) in MET
Returns:
(np.array, np.array): The Fermi azimuth, zenith
"""
quats = self.get_quaternions(times)
az, zen = coords.radec_to_spacecraft(ra, dec, quats)
return az, zen
[docs] def to_equatorial(self, fermi_az, fermi_zen, times):
"""Convert a position in spacecraft coordinates to equatorial coordinates
Args:
fermi_az (float): The Fermi azimuth of a position
fermi_zen (float): The Fermi zenith of a position
times (float or np.array): Time(s) in MET
Returns:
(np.array, np.array): The RA, Dec of the position
"""
quats = self.get_quaternions(times)
ra, dec = coords.spacecraft_to_radec(fermi_az, fermi_zen, quats)
return ra, dec
[docs] def location_visible(self, ra, dec, times):
"""Determine if a sky location is visible or occulted by the Earth
Args:
ra (float): The RA of a sky position
dec (float): The Dec of a sky position
times (float or np.array): Time(s) in MET
Returns:
np.array(dtype=bool): A boolean array where True indicates the \
position is visible.
"""
geo = np.array(self.get_geocenter_radec(times)).reshape(2, -1)
angle = coords.haversine(geo[0, :], geo[1, :], ra, dec)
visible = (angle > self.get_earth_radius(times))
return visible
[docs] def detector_pointing(self, det, times):
"""Retrieve the pointing of a GBM detector in equatorial coordinates
Args:
det (str): The GBM detector
times (float or np.array): Time(s) in MET
Returns:
(np.array, np.array): The RA, Dec of the detector pointing
"""
det = det.lower()
if det not in self._detectors:
raise ValueError('Illegal detector name')
d = Detector.from_str(det)
az, zen = d.azimuth, d.zenith
ra, dec = self.to_equatorial(az, zen, times)
return ra, dec
[docs] def detector_angle(self, ra, dec, det, times):
"""Determine the angle between a sky location and a GBM detector
Args:
ra (float): The RA of a sky position
dec (float): The Dec of a sky position
det (str): The GBM detector
times (float or np.array): Time(s) in MET
Returns:
np.array: The angle, in degrees, between the position and \
detector pointing
"""
det_ra, det_dec = self.detector_pointing(det, times)
angle = coords.haversine(det_ra, det_dec, ra, dec)
return angle
[docs] @classmethod
def open(cls, filename):
"""Open and read a position history file
Args:
filename (str): The filename of the FITS file
Returns:
:class:`PosHist`: The PosHist object
"""
obj = cls()
obj._file_properties(filename)
# open FITS file
with fits.open(filename) as hdulist:
for hdu in hdulist:
obj._headers.update({hdu.name: hdu.header})
data = hdulist['GLAST POS HIST'].data
times = data['SCLK_UTC']
obj._times = times
obj._data = data
# set the interpolators
obj._set_interpolators()
return obj
[docs] @classmethod
def merge(cls, poshists):
"""Merge multiple PosHists into a single PosHist object
Args:
poshists (list of :class:`PosHist`): List of PosHist objects
Returns:
:class:`PosHist`: The merged PosHist object
"""
# sort by tstart, and remove overlapping times
idx = np.argsort([poshist.time_range[0] for poshist in poshists])
times = [poshists[idx[0]]._times]
data = [poshists[idx[0]]._data]
for i in idx[1:]:
mask = (poshists[i]._times > times[-1][-1])
times.append(poshists[i]._times[mask])
data.append(poshists[i]._data[mask])
# create object with merged data
obj = cls()
obj._file_properties(poshists[idx[0]].full_path)
obj._times = np.concatenate(times)
obj._data = stack_arrays(data, asrecarray=True, usemask=False)
obj._headers = poshists[idx[0]].headers
obj._set_interpolators()
return obj
def _set_interpolators(self):
times = self._times
data = self._data
# Earth inertial coordinates interpolator
eic = np.array((data['POS_X'], data['POS_Y'], data['POS_Z']))
self._eic_interp = interp1d(times, eic)
# quaternions interpolator
quat = np.array(
(data['QSJ_1'], data['QSJ_2'], data['QSJ_3'], data['QSJ_4']))
self._quat_interp = interp1d(times, quat)
# Orbital position interpolator
# mark TODO: poshist uses the "simple" version of lat/lon calc
if 'SC_LAT' in data.dtype.names:
self._lat_interp = interp1d(times, data['SC_LAT'])
self._lon_interp = interp1d(times, data['SC_LON'])
alt = self._altitude_from_scpos(eic)
else:
lat, lon, alt = self._lat_lon_from_scpos(times, eic, complex=False)
self._lat_interp = interp1d(times, lat)
self._lon_interp = interp1d(times, lon)
self._alt_interp = interp1d(times, alt)
# Earth radius and geocenter interpolators
self._earth_radius_interp = interp1d(times, self._geo_half_angle(alt))
self._geocenter_interp = interp1d(times,
coords.geocenter_in_radec(eic))
# Velocity interpolator
vel = np.array((data['VEL_X'], data['VEL_Y'], data['VEL_Z']))
self._vel_interp = interp1d(times, vel)
# Angular velocity interpolator
angvel = np.array((data['WSJ_1'], data['WSJ_2'], data['WSJ_3']))
self._angvel_interp = interp1d(times, angvel)
# mark FIXME: the FLAGS=0 and FLAGS=2 never appear in daily poshist
# Interpolators for sun visibility and SAA passage
sun_visible = ((data['FLAGS'] == 1) | (data['FLAGS'] == 3))
in_saa = ((data['FLAGS']) == 2 | (data['FLAGS'] == 3))
self._sun_interp = interp1d(times, sun_visible)
self._saa_interp = interp1d(times, in_saa)
# set GTI based on SAA passages
# also times where sun is occulted
self._gti = self._split_bool_mask(in_saa, times)
self._sun_occulted = self._split_bool_mask(sun_visible, times)
def _from_trigdat(self, times, quats, scpos):
"""Initialize the PosHist object with info from a Trigdat object
Args:
times (np.array): An array of Trigdat bin start times
quats (np.array): An array of quaternions from the Trigdat
scpos (np.array): An array of EICs from the Trigdat
"""
if quats.shape[0] == times.shape[0]:
quats = quats.T
if scpos.shape[0] == times.shape[0]:
scpos = scpos.T
# EIC is in km in Trigdat. We need it in m.
scpos *= 1000.0
# interpolators for EIC and quaternions
self._eic_interp = interp1d(times, scpos, fill_value='extrapolate')
self._quat_interp = interp1d(times, quats, fill_value='extrapolate')
# orbital position interpolators
lat, lon, alt = self._lat_lon_from_scpos(times, scpos, complex=True)
self._lat_interp = interp1d(times, lat, fill_value='extrapolate')
self._lon_interp = interp1d(times, lon, fill_value='extrapolate')
self._alt_interp = interp1d(times, alt, fill_value='extrapolate')
# georadius and geocenter interpolators
self._earth_radius_interp = interp1d(times, self._geo_half_angle(alt),
fill_value='extrapolate')
self._geocenter_interp = interp1d(times,
coords.geocenter_in_radec(scpos),
fill_value='extrapolate')
# velocity and angular velocity interpolators
vel = self._velocity_from_scpos(times, scpos)
self._vel_interp = interp1d(times, vel, fill_value='extrapolate')
angvel = self._angular_velocity_from_quaternion(times, quats)
self._angvel_interp = interp1d(times, angvel, fill_value='extrapolate')
# sun visibility
sun_visible = self._sun_visible_from_times(times)
self._sun_interp = interp1d(times, sun_visible,
fill_value='extrapolate')
self._sun_occulted = self._split_bool_mask(sun_visible, times)
def _split_bool_mask(self, mask, times):
"""Split a boolean mask into segments of contiguous values and apply
to array of times
Args:
mask (np.array(dtype=bool)): The boolean mask
times (np.array): The array of times. Must be the same size as mask
Returns
[(float, float),...]: A list of time ranges
"""
# split a boolean mask array into segments based on True/False
indices = np.nonzero(mask[1:] != mask[:-1])[0] + 1
time_segs = np.split(times, indices)
mask_segs = np.split(mask, indices)
# retrieve the start and stop times for the "off" intervals
segs = []
numsegs = len(indices) + 1
for i in range(numsegs):
if mask_segs[i][0] == 0:
segs.append((time_segs[i][0], time_segs[i][-1]))
# mark FIXME: null if mask is all True or all False
# currently assuming that it must be all True
if len(segs) == 0:
segs = [self.time_range]
return segs
def _altitude_from_scpos(self, scpos):
"""Calculate altitude from the EIC.
Will attempt to do the more complex calculation taking into account
the shape of the Earth and the variable rotational velocity of the
Earth. This requires having up-to-date IERS tables. If local tables
are not up-to-date and astropy cannot query for a new table, will
default to a simpler approximation (spherical Earth, constant rotational
velocity).
Args:
scpos (np.array): The Earth inertial coordinates
Returns:
np.array: The altitudes
"""
try:
_, alt = coords.latitude_from_geocentric_coords_complex(scpos)
except:
warn('Using simple spheroidal Earth approximation')
_, alt = coords.latitude_from_geocentric_coords_simple(scpos)
return alt
def _angular_velocity_from_quaternion(self, times, quats):
"""Calculate angular velocity from changes in quaternion over time.
This is accurate for small rotations over very short times.
Args:
times (np.array): Times in MET
quats (np.array): Quaternions; one for each time
Returns:
np.array: The angular velocity
"""
dt = times[1:] - times[0:-1]
dquat = quats[:, 1:] - quats[:, :-1]
prod_quat = dquat
for i in range(len(dt)):
conj_quat = coords.quaternion_conj(quats[:, i])
prod_quat[:, i] = coords.quaternion_prod(conj_quat, dquat[:, i])
ang_vel = 2.0 / dt[np.newaxis, :] * prod_quat[0:3, :]
ang_vel = np.append(ang_vel, ang_vel[:, -1:], axis=1)
return ang_vel
def _velocity_from_scpos(self, times, scpos):
"""Calculate velocity from changes in position over time.
Args:
times (np.array): Times in MET
scpos (np.array): Positions of the spacecraft in Earth inertial
coordinates
Returns:
np.array: The velocity
"""
dt = times[1:] - times[0:-1]
dpos = scpos[:, 1:] - scpos[:, 0:-1]
velocity = dpos / dt
# copy last entry
velocity = np.append(velocity, velocity[:, -1:], axis=1)
return velocity
def _lat_lon_from_scpos(self, times, scpos, complex=False):
"""Convert coordinates of the spacecraft in Earth inertial coordinates
to latitude, East longitude, and altitude.
Args:
times (np.array): Times in MET
scpos (np.array): Positions of the spacecraft in Earth inertial
coordinates
complex (bool, optional):
If True, then uses the non-spherical Earth model and accurate
Earth rotation model. Default is False; Earth is a sphere and
the Earth rotation model is less accurate.
GBM FSW uses the latter option.
Returns:
tuple - 3-tuple containing:
- *np.array*: latitude
- *np.array*: longitude
- *np.array*: altitude
"""
if complex is False:
lat, alt = coords.latitude_from_geocentric_coords_simple(scpos)
else:
lat, alt = coords.latitude_from_geocentric_coords_complex(scpos)
try:
lon = coords.longitude_from_geocentric_coords(scpos, times,
ut1=complex)
except:
lon = coords.longitude_from_geocentric_coords(scpos, times,
ut1=False)
return (lat, lon, alt)
def _geo_half_angle(self, altitudes):
"""Calculate the apparent angular radius of the Earth
Args:
altitudes (np.array): The altitude of Fermi
Returns:
np.array: The angular radius, in degrees
"""
r = 6371.0 * 1000.0
half_angle = np.rad2deg(np.arcsin(r / (r + altitudes)))
return half_angle
def _sun_visible_from_times(self, times):
"""Calculate the sun visibility mask from the MET
Args:
times (np.array): The METs
Returns:
np.array(dtype=bool): A boolean mask, where True indicates the sun \
is visible.
"""
sun = np.array(coords.get_sun_loc(times))
return self.location_visible(sun[0, :], sun[1, :], times)