Source code for sunkit_instruments.fermi.fermi

"""
This module provides processing routines for Fermi Gamma-ray Space Telescope
(FGST), formerly called the Gamma-ray Large Area Space Telescope (GLAST).
"""
import os
import copy
import urllib
import tempfile
from collections import OrderedDict

import matplotlib.pyplot as plt
import numpy as np

import astropy.units as u
from astropy.coordinates import Latitude, Longitude
from astropy.time import TimeDelta
from sunpy.coordinates import sun
from sunpy.io.fits import fits
from sunpy.time import TimeRange, parse_time
from sunpy.time.time import _variables_for_parse_time_docstring
from sunpy.util.decorators import add_common_docstring

__all__ = [
    "download_weekly_pointing_file",
    "get_detector_sun_angles_for_time",
    "get_detector_sun_angles_for_date",
    "plot_detector_sun_angles",
    "met_to_utc",
]


[docs]@add_common_docstring(**_variables_for_parse_time_docstring()) def download_weekly_pointing_file(date): """ Downloads the FERMI/LAT weekly pointing file corresponding to the specified date. This file contains 1 minute cadence data on the spacecraft pointing, useful for calculating detector angles. Parameters ---------- date : {parse_time_types} A date specified as a parse_time-compatible time string, number, or a datetime object. Returns ------- `str`: The filepath to the downloaded file. """ date = parse_time(date) # use a temp directory to hold the file tmp_dir = tempfile.mkdtemp() # use Fermi data server to access weekly LAT pointing file. base_url = "https://heasarc.gsfc.nasa.gov/FTP/fermi/data/lat/weekly/spacecraft/" fbasename = "lat_spacecraft_weekly_w" # find out which file to get based on date # earliest full file in the FERMI server is for mission week 10, # beginning 2008 August 7. weekly_file_start = parse_time("2008-08-07") base_week = 10 # find out which mission week corresponds to date time_diff = date - weekly_file_start weekdiff = time_diff.to(u.day).value // 7 week = weekdiff + base_week # weekstr = ('%03.0f' % week) weekstr = f"{week:03.0f}" # construct the full url for the weekly pointing file full_fname = fbasename + weekstr + "_p202_v001.fits" pointing_file_url = base_url + full_fname # try to download the file try: # Use a context manager to avoid leaving a connection open with urllib.request.urlopen(pointing_file_url) as _: exists = True except urllib.error.HTTPError: exists = False # if no matches at all were found, then the pointing file doesn't exist if not exists: raise ValueError("No Fermi pointing files found for given date!") # download the file destination = os.path.join(tmp_dir, full_fname) urllib.request.urlretrieve(pointing_file_url, destination) # return the location of the downloaded file return destination
[docs]@add_common_docstring(**_variables_for_parse_time_docstring()) def get_detector_sun_angles_for_time(time, file): """ Get the GBM detector angles vs the Sun for a single time. Parameters ---------- time : {parse_time_types} A time specified as a parse_time-compatible time string, number, or a datetime object. file : `str` A filepath to a Fermi/LAT weekly pointing file (e.g. as obtained by the download_weekly_pointing_file function). Returns ------- `tuple`: A tuple of all the detector angles. """ time = parse_time(time) scx, scz, tt = get_scx_scz_at_time(time, file) # retrieve the detector angle information in spacecraft coordinates detectors = nai_detector_angles() # get the detector pointings in RA/DEC given the input spacecraft x and z # axes detector_radecs = nai_detector_radecs(detectors, scx, scz, tt) # this gets the sun position with RA in hours in decimal format (e.g. 4.3). # DEC is already in degrees sunpos_ra_not_in_deg = [ sun.apparent_rightascension(time), sun.apparent_declination(time), ] # now Sun position with RA in degrees sun_pos = [sunpos_ra_not_in_deg[0].to("deg"), sunpos_ra_not_in_deg[1]] # sun_pos = [(sunpos_ra_not_in_deg[0] / 24) * 360., sunpos_ra_not_in_deg[1]] # now get the angle between each detector and the Sun detector_to_sun_angles = get_detector_separation_angles(detector_radecs, sun_pos) return detector_to_sun_angles
[docs]@add_common_docstring(**_variables_for_parse_time_docstring()) def get_detector_sun_angles_for_date(date, file): """ Get the GBM detector angles vs the Sun as a function of time for a given date. Parameters ---------- date : {parse_time_types} A date specified as a parse_time-compatible time string, number, or a datetime object. file : `str` A filepath to a Fermi/LAT weekly pointing file (e.g. as obtained by the download_weekly_pointing_file function). Returns ------- `tuple`: A tuple of all the detector angles. """ date = parse_time(date) tran = TimeRange(date, date + TimeDelta(1 * u.day)) scx, scz, times = get_scx_scz_in_timerange(tran, file) # retrieve the detector angle information in spacecraft coordinates detectors = nai_detector_angles() detector_to_sun_angles = [] # get the detector vs Sun angles for each t and store in a list of # dictionaries. for i in range(len(scx)): detector_radecs = nai_detector_radecs(detectors, scx[i], scz[i], times[i]) # this gets the sun position with RA in hours in decimal format # (e.g. 4.3). DEC is already in degrees sunpos_ra_not_in_deg = [ sun.apparent_rightascension(times[i]), sun.apparent_declination(times[i]), ] # now Sun position with RA in degrees sun_pos = [sunpos_ra_not_in_deg[0].to("deg"), sunpos_ra_not_in_deg[1]] # now get the angle between each detector and the Sun detector_to_sun_angles.append( get_detector_separation_angles(detector_radecs, sun_pos) ) # slice the list of dictionaries to get the angles for each detector in a # list form angles = OrderedDict() key_list = [ "n0", "n1", "n2", "n3", "n4", "n5", "n6", "n7", "n8", "n9", "n10", "n11", "time", ] for i in range(13): if not key_list[i] == "time": angles[key_list[i]] = [ item[key_list[i]].value for item in detector_to_sun_angles ] * u.deg else: angles[key_list[i]] = [item[key_list[i]] for item in detector_to_sun_angles] return angles
[docs]def plot_detector_sun_angles(angles): """ Plots the Fermi/GBM detector angles as a function of time. Parameters ---------- angles : `dict` A dictionary containing the Fermi/GBM detector angle information as a function of time. Obtained from the `~sunkit_instruments.fermi.get_detector_separation_angles` function. """ # make a plot showing the angles vs time figure = plt.figure(1) for n in angles.keys(): if not n == "time": plt.plot( angles["time"], angles[n].value, label="{lab} ({val})".format( lab=n, val=str(np.mean(angles[n].value))[0:5] ), ) plt.ylim(180, 0) plt.ylabel("angle (degrees)") plt.xlabel("Start time: " + angles["time"][0].isoformat()) plt.title("Detector pointing angle from Sun") plt.legend(fontsize=10) figure.autofmt_xdate() plt.show()
@add_common_docstring(**_variables_for_parse_time_docstring()) def get_scx_scz_at_time(time, file): """ Read a downloaded FERMI weekly pointing file and extract "scx", "scz" for a single time. Parameters ---------- time : {parse_time_types} A time specified as a parse_time-compatible time string, number, or a datetime object. file : `str` A filepath to a Fermi/LAT weekly pointing file (e.g. as obtained by the `~sunkit_instruments.fermi.download_weekly_pointing_file` function). Returns ------- `tuple`, `tuple`, `list`: The pointing coordinates as a `~astropy.coordinates.Longitude` in a `tuple` and it's time. """ time = parse_time(time) hdulist = fits.open(file) timesinutc = [] for tim in hdulist[1].data["START"]: timesinutc.append(met_to_utc(tim)) ind = np.searchsorted(timesinutc, time) scx_radec = ( Longitude(hdulist[1].data["RA_SCX"][ind] * u.deg), Latitude(hdulist[1].data["DEC_SCX"][ind] * u.deg), ) scz_radec = ( Longitude(hdulist[1].data["RA_SCZ"][ind] * u.deg), Latitude(hdulist[1].data["DEC_SCZ"][ind] * u.deg), ) return scx_radec, scz_radec, timesinutc[ind] def get_scx_scz_in_timerange(timerange, file): """ Read a downloaded FERMI weekly pointing file and extract scx, scz for a timerange. Parameters ---------- timerange : `sunpy.time.TimeRange` A SunPy `~sunpy.time.TimeRange`. file : `str` A filepath to a Fermi/LAT weekly pointing file (e.g. as obtained by the `~sunkit_instruments.fermi.download_weekly_pointing_file` function). Returns ------- `list`, `list`, `list`: The pointing coordinates as a `~astropy.coordinates.Longitude` in a `list` and it's time. """ hdulist = fits.open(file) timesinutc = [] for tim in hdulist[1].data["START"]: timesinutc.append(met_to_utc(tim)) startind = np.searchsorted(timesinutc, timerange.start) endind = np.searchsorted(timesinutc, timerange.end) scx_radec = [] scz_radec = [] for i in range(startind, endind): scx_radec.append( ( Longitude(hdulist[1].data["RA_SCX"][i] * u.deg), Latitude(hdulist[1].data["DEC_SCX"][i] * u.deg), ) ) scz_radec.append( ( Longitude(hdulist[1].data["RA_SCZ"][i] * u.deg), Latitude(hdulist[1].data["DEC_SCZ"][i] * u.deg), ) ) return scx_radec, scz_radec, timesinutc[startind:endind] def nai_detector_angles(): """ Returns the dictionary of Fermi/GBM NAI detector zenith and azimuth angles, in spacecraft coordinates. Zenith angle is measured from "+z" (along the LAT boresight), azimuth is measured from "+x". References ---------- Meegan, Charles, et al. "The Fermi gamma-ray burst monitor." The Astrophysical Journal 702.1 (2009): 791. """ # angles listed as [azimuth, zenith] detectors = { "n0": [45.89 * u.deg, 20.58 * u.deg], "n1": [45.11 * u.deg, 45.31 * u.deg], "n2": [58.44 * u.deg, 90.21 * u.deg], "n3": [314.87 * u.deg, 45.24 * u.deg], "n4": [303.15 * u.deg, 90.27 * u.deg], "n5": [3.35 * u.deg, 89.79 * u.deg], "n6": [224.93 * u.deg, 20.43 * u.deg], "n7": [224.62 * u.deg, 46.18 * u.deg], "n8": [236.61 * u.deg, 89.97 * u.deg], "n9": [135.19 * u.deg, 45.55 * u.deg], "n10": [123.73 * u.deg, 90.42 * u.deg], "n11": [183.74 * u.deg, 90.32 * u.deg], } return detectors @add_common_docstring(**_variables_for_parse_time_docstring()) def nai_detector_radecs(detectors, scx, scz, time): """ calculates the "RA/DEC" for each NaI detector given spacecraft "z" and "x" "RA/DEC" positions. This routine is based on code found in `GTBURST <https://fermi.gsfc.nasa.gov/ssc/data/analysis/scitools/gtburst.html>`__, originally written by Dr Giacamo Vianello for the Fermi Science Tools. Parameters ---------- detectors : `dict` A dictionary containing the Fermi/GBM detector pointing angles relative to the spacecraft axes. Obtained from the `sunkit_instruments.fermi.nai_detector_angles` function. scx : array-like Two-element tuple containing the "RA/DEC" information of the Fermi spacecraft X-axis scz : array-like Two-element tuple containing the "RA/DEC" information of the Fermi spacecraft Z-axis time : {parse_time_types} A time specified as a parse_time-compatible time string, number, or a datetime object. This will correspond to the input ``scx`` and ``scz`` values. Returns ------- `dict` A dictionary containing the "RA/DEC" for each Fermi/GBM NaI detector at the given input time. """ scx_vector = np.array( [ np.cos(scx[0].to("rad").value) * np.cos(scx[1].to("rad").value), np.sin(scx[0].to("rad").value) * np.cos(scx[1].to("rad").value), np.sin(scx[1].to("rad").value), ] ) scz_vector = np.array( [ np.cos(scz[0].to("rad").value) * np.cos(scz[1].to("rad").value), np.sin(scz[0].to("rad").value) * np.cos(scz[1].to("rad").value), np.sin(scz[1].to("rad").value), ] ) # For each detector, do the rotation depending on the detector zenith and # azimuth angles. detector_radecs = copy.deepcopy(detectors) for l, d in detectors.items(): phi = d[0].value theta = d[1].value # rotate about spacecraft z-axis first vx_primed = rotate_vector(scx_vector, scz_vector, np.deg2rad(phi)) # now find spacecraft y-axis using cross product vy_primed = np.cross(scz_vector, vx_primed) # do the second part of the rotation around vy vz_primed = rotate_vector(scz_vector, vy_primed, np.deg2rad(theta)) # now we should be pointing at the new RA/DEC. ra = Longitude(np.degrees(np.arctan2(vz_primed[1], vz_primed[0])) * u.deg) dec = Latitude(np.degrees(np.arcsin(vz_primed[2])) * u.deg) # save the RA/DEC in a dictionary detector_radecs[l] = [ra, dec] detector_radecs["time"] = time return detector_radecs def rotate_vector(vector, axis, theta): """ The Euler-Rodrigues formula for rotating vectors. Parameters ---------- vector : `numpy.ndarray` A three-element vector to be rotated. axis : `numpy.ndarray` The three-element vector to rotate around. theta : `float` The angle (in radians) by which to rotate vector around axis. Reference --------- https://en.wikipedia.org/wiki/Euler-Rodrigues_parameters#Rotation_angle_and_rotation_axis """ axis = axis / np.sqrt(np.dot(axis, axis)) a = np.cos(theta / 2) b, c, d = -axis * np.sin(theta / 2) rot_matrix = np.array( [ [a * a + b * b - c * c - d * d, 2 * (b * c + a * d), 2 * (b * d - a * c)], [2 * (b * c - a * d), a * a + c * c - b * b - d * d, 2 * (c * d + a * b)], [2 * (b * d + a * c), 2 * (c * d - a * b), a * a + d * d - b * b - c * c], ] ) return np.dot(rot_matrix, vector) def get_detector_separation_angles(detector_radecs, sunpos): """ Finds the separation angle between the Sun and each NaI detector, given a dictionary of detector "RA/DEC"s. Parameters ---------- detector_radecs : `dict` The "RA/DEC" for each NaI detector as Astropy quantities. Obtained from the `sunkit_instruments.fermi.nai_detector_radecs` function sunpos : `list` Two-element list containing the "RA/DEC" of the Sun position as `astropy.unit.quantity`, e.g., ``[<Longitude 73.94 deg>, <Latitude 22.66 deg>]`` """ angles = copy.deepcopy(detector_radecs) for l, d in detector_radecs.items(): if not l == "time": angle = separation_angle(d, sunpos) angles[l] = angle return angles def separation_angle(radec1, radec2): """ Use the law of spherical cosines to calculate the separation angle between two "RA/DEC" positions. Parameters ---------- radec1 : `list` Two-element list containing the "RA/DEC" position as `astropy.unit.quantity`, e.g., ``[<Longitude 73.94 deg>, <Latitude 22.66 deg>]`` radec2 : `list` Two-element list containing the "RA/DEC" position as `astropy.unit.quantity`, e.g., ``[<Longitude 73.94 deg>, <Latitude 22.66 deg>]`` """ cosine_of_angle = ( np.cos(((90 * u.deg) - radec1[1].to("degree")).to("rad")) * np.cos((90 * u.deg - radec2[1].to("degree")).to("rad")) ) + ( np.sin(((90 * u.deg) - radec1[1].to("degree")).to("rad")) * np.sin(((90 * u.deg) - radec2[1].to("degree")).to("rad")) * np.cos((radec1[0].to("degree") - radec2[0].to("degree")).to("rad")) ) angle = (np.arccos(cosine_of_angle)).to("degree") return angle
[docs]def met_to_utc(timeinsec): """ Converts Fermi Mission Elapsed Time (MET) in seconds to a `~astropy.time.Time` object. Parameters ---------- timeinsec : `float` Time in seconds since "00:00 UT" on 1st January 2001 (the Fermi MET format). Returns ------- `astropy.time.Time` The input Fermi Mission Elapsed Time converted to a `~astropy.time.Time` object. """ # Times for GBM are in Mission Elapsed Time (MET). # The reference time for this is 2001-Jan-01 00:00. met_ref_time = parse_time("2001-01-01 00:00") return met_ref_time + timeinsec * u.second
@add_common_docstring(**_variables_for_parse_time_docstring()) def utc_to_met(time_ut): """ Converts a UT to a Fermi Mission Elapsed Time (MET) float. Parameters ---------- time_ut : {parse_time_types} A time specified as a parse_time-compatible time string, number, or a datetime object. Returns ------- `astropy.units.Quantity` The Fermi Mission Elapsed Time corresponding to the input UT """ met_ref_time = parse_time("2001-01-01 00:00") return (time_ut - met_ref_time).to(u.second)