# -*- coding: utf-8 -*-
import hashlib
import os
import pickle
from inspect import getfullargspec as getargspec
from pathlib import Path
from urllib.request import urlretrieve
import astropy.constants as const
import astropy.units as u
import numpy as np
from astropy.time import Time
from scipy.interpolate import CubicSpline
from tqdm import tqdm
import EXOSIMS.Prototypes.Observatory
from EXOSIMS.util._numpy_compat import copy_if_needed
from EXOSIMS.util.eccanom import eccanom
from EXOSIMS.util.get_dirs import get_cache_dir, get_downloads_dir
from EXOSIMS.util.vprint import vprint
[docs]
class Observatory(object):
""":ref:`Observatory` Prototype
Args:
SRP (bool):
Toggle solar radiation pressure. Defaults True.
koAngles_SolarPanel (list(float)):
[Min, Max] keepout angles (in degrees) due to solar panels.
Defaults to [0,180].
ko_dtStep (float):
Step size to use when calculating keepout maps (in days). Defaults to 1.
settlingTime (float):
Observatory settling time after retargeting (in days). Defaults to 1. This
time is added to every observation and counts againts the total integration
time allocation.
thrust (float):
Slew thrust mangitude (in mN). Defaults to 450 mN.
slewIsp (float):
Slew specific impulse (in seconds). Defaults to 4160 s.
scMass (float):
Maneuvering spacecraft initial wet mass (in kg). Nominally this is the
starshade, but may also be the observatory if the starshade is kept on the
stable orbit. Defaults to 6000 kg.
.. warning::
If ``twotanks`` is true, this input will be ignored and attribute
``scMass`` will be initially set to the sum of ``slewMass`` and
``skMass``.
slewMass (float):
Initial fuel mass of slewing propulsion system (in kg). Defaults to 0. Only
used if twotanks is True.
skMass (float):
Initial fuel mass of stationkeeping propulsion system (in kg).
Defaults to 0. Only used if twotanks is True.
twotanks (bool):
Determines whether stationkeeping and slewing propulsion systems use
separate tanks. If False, it is assumed that all onboard fuel is fungible.
Defaults False.
skEff (float):
Stationkeeping propulsion system efficiency. Must be between 0 and 1.
Defaults to 0.7098 (approximately 45 deg cosine losses).
slewEff (float):
Slewing propulsion system efficiency. Must be between 0 and 1.
Defaults to 1.
dryMass (float):
Maneuvering spacecraft dry mass (in kg). Defaults to 3400 kg.
Must be smaller than scMass.
coMass (float):
Non-manuevering spacecraft (nominally the observatory) initial wet mass
(in kg). Defaults to 5800 kg.
occulterSep (float):
Initial occulter separation (in km). Defaults to 55000.
skIsp (float):
Stationkeeping propulsion system specific impulse (in seconds).
Defaults to 220 s.
defburnPortion (float):
Default burn portion for simple model slews. Must be between 0 and 1.
Defaults to 0.05.
constTOF (float):
Constant time of flight value (in days). Defaults to 14. DEPRECATED
maxdVpct (float):
Maximum delta V percentage allowed for any maneuver.
Must be between 0 and 1. Defaults to 0.02.
spkpath (str, optional):
Path to SPK file on disk.
If not set, defaults to de432s.bsp in :ref:`EXOSIMSDOWNLOADS`.
checkKeepoutEnd (bool):
Check keepout conditions at end of observation. Defaults True.
TODO: Move to SurveySimulation
forceStaticEphem (bool):
Use static ephemerides for solar system objects instead of jplephem.
Defaults False.
occ_dtmin (float):
Minimum slew time (in days). Defaults to 0.055
occ_dtmax (float):
Maximum slew time (in days). Defaults to 61
sk_Tmin (float):
Minimum time after mission start to compute stationkeeping (in days).
Defaults to 0.
sk_Tmax (float):
Maximum time after mission start to compute stationkeeping (in days).
Defaults to 365
non_lambertian_coefficient_front (float):
Non-Lambertion reflectivity coefficient of front face of manuevering
spacecraft. Used for SRP calculations. Defaults to 0.038.
non_lambertian_coefficient_back (float):
Non-Lambertion reflectivity coefficient of back face of manuevering
spacecraft. Used for SRP calculations. Defaults to 0.004.
specular_reflection_factor (float):
Specular reflectivity of maneuvering spacecraft. Used for SRP calculations.
Defaults to 0.975.
nreflection_coefficient (float):
non-specular reflectivity of maneuvering spacecraft.
Used for SRP calculations. Defaults to 0.999
emission_coefficient_front (float):
Emission coefficient of front face of maneuvering spacecraft.
Used for SRP calculations. Defaults to 0.8
emission_coefficient_back (float):
Emission coefficient of rear face of maneuvering spacecraft.
Used for SRP calculations. Defaults to 0.2
allowRefueling (bool):
Fuel tanks can be topped off when they reach empty from external fuel source
whose capacity is defined by the ``external_fuel_mass`` input.
Defaults False.
external_fuel_mass (float):
Initial mass of external fuel supply (in kg). Ignored if ``allowRefueling``
is False. Defaults to 0.
cachedir (str, optional):
Full path to cachedir.
If None (default) use default (see :ref:`EXOSIMSCACHE`)
**specs:
:ref:`sec:inputspec`
Attributes:
_outspec (dict):
:ref:`sec:outspec`
ao (astropy.units.quantity.Quantity):
Thurst acceleration (current thrust/spacecraft mass). Acceleration units.
cachedir (str):
Path to the EXOSIMS cache directory (see :ref:`EXOSIMSCACHE`)
checkKeepoutEnd (bool):
Toggle checking of keepout at end of observations (as well as the
beginning). TODO: need to depricate in favor of continuous visibility.
coMass (astropy.units.quantity.Quantity):
Non-manuevering spacecraft (nominally the observatory) wet mass. Mass units.
constTOF (astropy.units.quantity.Quantity):
Constant time of flight for single occulter slew. DEPRECATED
defburnPortion (float):
Default burn portion for simple slew model.
dryMass (astropy.units.quantity.Quantity):
Maneuvering spacecraft dry mass. Mass units.
dVmax (astropy.units.quantity.Quantity):
Maximum single-slew allowable delta V (as determined by maxdVpct input.
Units of velocity.
dVtot (astropy.units.quantity.Quantity):
Total possible slew delta V as determined by ideal rocket equation applied
to the initial fuel mass.
emission_coefficient_back (float):
Emission coefficient of back face of maneuvering spacecraft.
Used for SRP calculations.
emission_coefficient_front (float):
Emission coefficient of front face of maneuvering spacecraft.
Used for SRP calculations.
flowRate (astropy.units.quantity.Quantity):
Slew ropulsion system mass flow rate. Units: mass/time
forceStaticEphem (bool):
Use static ephemerides for solar system objects instead of jplephem.
havejplephem (bool):
jplephem module installed and SPK available.
kernel (jplephem.spk.SPK):
jplephem kernel used for ephemeris calculations of solar system bodies
ko_dtStep (astropy.units.quantity.Quantity):
Step size to use when calculating keepout maps. Time units.
koAngles_SolarPanel (astropy.units.quantity.Quantity):
[Min, Max] keepout angles due to solar panels.
maxdVpct (float):
Maximum delta V percentage allowed for any maneuver.
maxFuelMass (astropy.units.quantity.Quantity):
Total capacity of all fuel. This parameter is constant and should never
be modified externally.
non_lambertian_coefficient_back (float):
Non-Lambertion reflectivity coefficient of back face of manuevering
spacecraft. Used for SRP calculations.
non_lambertian_coefficient_front (float):
Non-Lambertion reflectivity coefficient of front face of manuevering
spacecraft. Used for SRP calculations.
nreflection_coefficient (float):
Non-specular reflectivity of maneuvering spacecraft.
Used for SRP calculations.
occ_dtmax (astropy.units.quantity.Quantity):
Maximum allowable slew time
occ_dtmin (astropy.units.quantity.Quantity):
Minimum allowable slew time.
occulterSep (astropy.units.quantity.Quantity):
Current occulter separation distance.
scMass (astropy.units.quantity.Quantity):
Current maneuvering spacecraft mass.
settlingTime (astropy.units.quantity.Quantity):
Observatory settling time after every retargeting.
sk_Tmax (astropy.units.quantity.Quantity):
Maximum time after mission start to compute stationkeeping
sk_Tmin (astropy.units.quantity.Quantity):
Minimum time after mission start to compute stationkeeping
skEff (float):
Stationkeeping propulsion system efficience.
skIsp (astropy.units.quantity.Quantity):
Stationkeeping propulsion system specific impulse. Time units.
skMass (astropy.units.quantity.Quantity):
Stationkeeping propulsion system fuel mass.
skMaxFuelMass (astropy.units.quantity.Quantity):
Total capacity of stationkeeping fuel tank. This parameter is constant
and should never be modified externally. Set to 0 if allowRefueling is
False.
slewEff (float):
Slew propulsion system efficiencey.
slewIsp (astropy.units.quantity.Quantity):
Slew propulsion system specific impulse. Time units.
slewMass (astropy.units.quantity.Quantity):
Slew propulsion system fuel mass.
slesMaxFuelMass (astropy.units.quantity.Quantity):
Total capacity of slewing fuel tank. This parameter is constant
and should never be modified externally. Set to 0 if allowRefueling is
False.
specular_reflection_factor (float):
Specular reflectivity of maneuvering spacecraft. Used for SRP calculations.
spkpath (str):
Full path to SPK file used by jplephem.
SRP (bool):
Toggles whether solar radiation pressure is included in calculations.
thrust (astropy.units.quantity.Quantity):
Slew propulsion system thrust. Force units.
twotanks (bool):
Toggles whether stationkeeping and slewing fuel are bookkept separately. If
false, all fuel is fungible.
.. note::
For finding positions of solar system bodies, this routine will attempt to
use the jplephem module and a local SPK file on disk. The module can be
installed via pip or from source. The default SPK file (which the code
attempts to automatically download) can be downloaded manually from:
http://naif.jpl.nasa.gov/pub/naif/generic_kernels/spk/planets/de432s.bsp
and should be placed in the :ref:`EXOSIMSDOWNLOADS` (or another path, specified
by the ``spkpath`` input).
"""
_modtype = "Observatory"
def __init__(
self,
SRP=True,
koAngles_SolarPanel=[0, 180],
ko_dtStep=1,
settlingTime=0.042,
thrust=450,
slewIsp=4160.0,
scMass=6000.0,
slewMass=0.0,
skMass=0.0,
twotanks=False,
skEff=0.7098,
slewEff=1.0,
dryMass=3400.0,
coMass=5800.0,
occulterSep=55000.0,
skIsp=220.0,
defburnPortion=0.05,
constTOF=14,
maxdVpct=0.02,
spkpath=None,
checkKeepoutEnd=True,
forceStaticEphem=False,
occ_dtmin=0.055,
occ_dtmax=61.0,
sk_Tmin=0.0,
sk_Tmax=365.0,
non_lambertian_coefficient_front=0.038,
non_lambertian_coefficient_back=0.004,
specular_reflection_factor=0.975,
nreflection_coefficient=0.999,
emission_coefficient_front=0.8,
emission_coefficient_back=0.2,
allowRefueling=False,
external_fuel_mass=0,
cachedir=None,
**specs,
):
# start the outspec
self._outspec = {}
# Defining constants
self.km2au = (1 * u.km).to_value(u.AU)
self.au2km = (1 * u.AU).to_value(u.km)
self.orbital_height = 42164.0 * self.km2au
self.orbital_inclination = np.radians(28.5)
self.orbital_right_ascension = np.radians(228.0)
self.orbital_frequency = 2 * np.pi
self.julian_century = 36525.0
# load the vprint function (same line in all prototype module constructors)
self.vprint = vprint(specs.get("verbose", True))
# validate inputs
assert isinstance(checkKeepoutEnd, bool), "checkKeepoutEnd must be a boolean."
assert isinstance(forceStaticEphem, bool), "forceStaticEphem must be a boolean."
# default Observatory values
self.SRP = SRP
# solar panel keepout angles
self.koAngles_SolarPanel = [float(x) for x in koAngles_SolarPanel] << u.deg
# time step for generating koMap of stars (day)
self.ko_dtStep = float(ko_dtStep) << u.d
# Observatory settling time after repoint
self.settlingTime = float(settlingTime) << u.d
# occulter slew thrust (mN)
self.thrust = float(thrust) << u.mN
# occulter slew specific impulse (s)
self.slewIsp = float(slewIsp) << u.s
# occulter initial (wet) mass (kg)
self.scMass = float(scMass) << u.kg
# slew fuel initial mass (kg)
self.slewMass = float(slewMass) << u.kg
# station keeping fuel initial mass (kg)
self.skMass = float(skMass) << u.kg
# boolean used to seperate manuevering fuel
self.twotanks = bool(twotanks)
# slew efficiency factor
self.slewEff = float(slewEff)
# station-keeping efficiency factor
self.skEff = float(skEff)
# occulter dry mass (kg)
self.dryMass = float(dryMass) << u.kg
# telescope mass (kg)
self.coMass = float(coMass) << u.kg
# occulter-telescope distance (km)
self.occulterSep = float(occulterSep) << u.km
# station-keeping Isp (s)
self.skIsp = float(skIsp) << u.s
# default burn portion
self.defburnPortion = float(defburnPortion)
# true if keepout called at obs end
self.checkKeepoutEnd = bool(checkKeepoutEnd)
# boolean used to force static ephemerides
self.forceStaticEphem = bool(forceStaticEphem)
# starshade constant slew time (days)
self.constTOF = np.array(constTOF, ndmin=1) << u.d
# Minimum occulter slew time (days)
self.occ_dtmin = float(occ_dtmin) << u.d
# Maximum occulter slew time (days)
self.occ_dtmax = float(occ_dtmax) << u.d
# Minimum days after missionstart to calculate stationkeeping (days)
self.sk_Tmin = float(sk_Tmin) << u.d
# Maximum days after missionstart to calculate stationkeeping (days)
self.sk_Tmax = float(sk_Tmax) << u.d
# Maximum deltaV percent
self.maxdVpct = float(maxdVpct)
# non-Lambertian coefficient (front)
self.non_lambertian_coefficient_front = float(non_lambertian_coefficient_front)
# non-Lambertian coefficient (back)
self.non_lambertian_coefficient_back = float(non_lambertian_coefficient_back)
# specular reflection factor
self.specular_reflection_factor = float(specular_reflection_factor)
# nreflection coefficient
self.nreflection_coefficient = float(nreflection_coefficient)
# emission coefficient (front)
self.emission_coefficient_front = float(emission_coefficient_front)
# emission coefficient (back)
self.emission_coefficient_back = float(emission_coefficient_back)
# Allow Refueling from external tank
self.allowRefueling = bool(allowRefueling)
# initial mass of external tank
self.external_fuel_mass = float(external_fuel_mass) << u.kg
# check that twotanks and dry mass add up to total mass
if self.twotanks:
assert (self.slewMass > 0) and (
self.skMass > 0
), "Tank mass must be positive"
self.scMass = self.slewMass + self.skMass + self.dryMass
# record tank capacities in case you need to refuel
self.skMaxFuelMass = self.skMass.copy()
self.slewMaxFuelMass = self.slewMass.copy()
else:
self.skMaxFuelMass = 0 << u.kg
self.slewMaxFuelMass = 0 << u.kg
# total tank capacity:
self.maxFuelMass = self.scMass - self.dryMass
assert (
self.maxFuelMass > 0 << u.kg
), "Initial spacecraft wet mass must be greater than dry mass."
# Acceleration
self.ao = self.thrust / self.scMass
# find the cache directory
self.cachedir = get_cache_dir(cachedir)
self._outspec["cachedir"] = self.cachedir
specs["cachedir"] = self.cachedir
# find amount of fuel on board starshade and an upper bound for single slew dV
self.dVtot = self.slewIsp * const.g0 * np.log(self.scMass / self.dryMass)
self.dVmax = self.dVtot * self.maxdVpct
# set values derived from quantities above
# slew flow rate (kg/day)
self.flowRate = (self.thrust / self.slewEff / const.g0 / self.slewIsp).to(
"kg/day"
)
# if jplephem is available, we'll use that for propagating solar system bodies
# otherwise, use static ephemerides
if self.forceStaticEphem is False:
try:
from jplephem.spk import SPK
self.havejplephem = True
except ImportError:
self.vprint(
"WARNING: Module jplephem not found, "
+ "using static solar system ephemerides."
)
self.havejplephem = False
else:
self.havejplephem = False
self.vprint("Using static solar system ephemerides.")
# define function for calculating obliquity of the ecliptic
# (arg Julian centuries from J2000)
self.obe = (
lambda TDB: 23.439279
- 0.0130102 * TDB
- 5.086e-8 * (TDB**2.0)
+ 5.565e-7 * (TDB**3.0)
+ 1.6e-10 * (TDB**4.0)
+ 1.21e-11 * (TDB**5.0)
)
# if you have jplephem, load spice file, otherwise load static ephem
if self.havejplephem:
if (spkpath is None) or not (os.path.exists(spkpath)):
# if the path does not exist, load the default de432s.bsp
filename = "de432s.bsp"
downloadsdir = get_downloads_dir()
spkpath = os.path.join(downloadsdir, filename)
# attempt to fetch ephemeris from NAIF
if not os.path.exists(spkpath) and os.access(
downloadsdir, os.W_OK | os.X_OK
):
spk_on_web = (
"https://naif.jpl.nasa.gov/pub/naif/generic_kernels/"
"spk/planets/de432s.bsp"
)
self.vprint(
"Fetching planetary ephemeris from %s to %s"
% (spk_on_web, spkpath)
)
urlretrieve(spk_on_web, spkpath)
self.kernel = SPK.open(spkpath)
else:
# All ephemeride data from Vallado Appendix D.4
# Values are: a = sma (AU), e = eccentricity, I = inclination (deg),
# O = long. ascending node (deg), w = long. perihelion (deg),
# lM = mean longitude (deg)
# store ephemerides data in heliocentric true ecliptic frame
a = 0.387098310
e = [0.20563175, 0.000020406, -0.0000000284, -0.00000000017]
I = [7.004986, -0.0059516, 0.00000081, 0.000000041] # noqa: E741
O = [48.330893, -0.1254229, -0.00008833, -0.000000196] # noqa: E741
w = [77.456119, 0.1588643, -0.00001343, 0.000000039]
lM = [252.250906, 149472.6746358, -0.00000535, 0.000000002]
Mercury = self.SolarEph(a, e, I, O, w, lM)
a = 0.723329820
e = [0.00677188, -0.000047766, 0.0000000975, 0.00000000044]
I = [3.394662, -0.0008568, -0.00003244, 0.000000010] # noqa: E741
O = [76.679920, -0.2780080, -0.00014256, -0.000000198] # noqa: E741
w = [131.563707, 0.0048646, -0.00138232, -0.000005332]
lM = [181.979801, 58517.8156760, 0.00000165, -0.000000002]
Venus = self.SolarEph(a, e, I, O, w, lM)
a = 1.000001018
e = [0.01670862, -0.000042037, -0.0000001236, 0.00000000004]
I = [0.0, 0.0130546, -0.00000931, -0.000000034] # noqa: E741
O = [174.873174, -0.2410908, 0.00004067, -0.000001327] # noqa: E741
w = [102.937348, 0.3225557, 0.00015026, 0.000000478]
lM = [100.466449, 35999.3728519, -0.00000568, 0.0]
Earth = self.SolarEph(a, e, I, O, w, lM)
a = 1.523679342
e = [0.09340062, 0.000090483, -0.0000000806, -0.00000000035]
I = [1.849726, -0.0081479, -0.00002255, -0.000000027] # noqa: E741
O = [49.558093, -0.2949846, -0.00063993, -0.000002143] # noqa: E741
w = [336.060234, 0.4438898, -0.00017321, 0.000000300]
lM = [355.433275, 19140.2993313, 0.00000261, -0.000000003]
Mars = self.SolarEph(a, e, I, O, w, lM)
a = [5.202603191, 0.0000001913]
e = [0.04849485, 0.000163244, -0.0000004719, -0.00000000197]
I = [1.303270, -0.0019872, 0.00003318, 0.000000092] # noqa: E741
O = [100.464441, 0.1766828, 0.00090387, -0.000007032] # noqa: E741
w = [14.331309, 0.2155525, 0.00072252, -0.000004590]
lM = [34.351484, 3034.9056746, -0.00008501, 0.000000004]
Jupiter = self.SolarEph(a, e, I, O, w, lM)
a = [9.554909596, -0.0000021389]
e = [0.05550862, -0.000346818, -0.0000006456, 0.00000000338]
I = [2.488878, 0.0025515, -0.00004903, 0.000000018] # noqa: E741
O = [113.665524, -0.2566649, -0.00018345, 0.000000357] # noqa: E741
w = [93.056787, 0.5665496, 0.00052809, 0.000004882]
lM = [50.077471, 1222.1137943, 0.00021004, -0.000000019]
Saturn = self.SolarEph(a, e, I, O, w, lM)
a = [19.218446062, -0.0000000372, 0.00000000098]
e = [0.04629590, -0.000027337, 0.0000000790, 0.00000000025]
I = [0.773196, -0.0016869, 0.00000349, 0.000000016] # noqa: E741
O = [74.005947, 0.0741461, 0.00040540, 0.000000104] # noqa: E741
w = [173.005159, 0.0893206, -0.00009470, 0.000000413]
lM = [314.055005, 428.4669983, -0.00000486, 0.000000006]
Uranus = self.SolarEph(a, e, I, O, w, lM)
a = [30.110386869, -0.0000001663, 0.00000000069]
e = [0.00898809, 0.000006408, -0.0000000008]
I = [1.769952, 0.0002257, 0.00000023, -0.000000000] # noqa: E741
O = [131.784057, -0.0061651, -0.00000219, -0.000000078] # noqa: E741
w = [48.123691, 0.0291587, 0.00007051, 0.0]
lM = [304.348665, 218.4862002, 0.00000059, -0.000000002]
Neptune = self.SolarEph(a, e, I, O, w, lM)
a = [39.48168677, -0.00076912]
e = [0.24880766, 0.00006465]
I = [17.14175, 0.003075] # noqa: E741
O = [110.30347, -0.01036944] # noqa: E741
w = [224.06676, -0.03673611]
lM = [238.92881, 145.2078]
Pluto = self.SolarEph(a, e, I, O, w, lM)
# store all as dictionary:
self.planets = {
"Mercury": Mercury,
"Venus": Venus,
"Earth": Earth,
"Mars": Mars,
"Jupiter": Jupiter,
"Saturn": Saturn,
"Uranus": Uranus,
"Neptune": Neptune,
"Pluto": Pluto,
}
self.spkpath = spkpath
# populate outspec
inputatts = getargspec(EXOSIMS.Prototypes.Observatory.Observatory.__init__)[0]
if "self" in inputatts:
inputatts.remove("self")
for att in inputatts:
dat = self.__dict__[att]
self._outspec[att] = dat.value if isinstance(dat, u.Quantity) else dat
self.j2000_jd = Time(2000.0, format="jyear", scale="tai").jd
self.has_earth_pos_interp = False
self.earth_pos_interp_range_mjd = [None, None]
def __del__(self):
"""destructor method. only here to clean up SPK kernel if it exists."""
if "kernel" in self.__dict__:
if self.kernel:
self.kernel.close()
def __str__(self):
"""String representation of the Observatory object
When the command 'print' is used on the Observatory object, this method
will print the attribute values contained in the object
"""
for att in self.__dict__:
print("%s: %r" % (att, getattr(self, att)))
return "Observatory class object attributes"
[docs]
def equat2eclip(self, r_equat, currentTime, rotsign=1):
"""Rotates heliocentric coordinates from equatorial to ecliptic frame.
Args:
r_equat (~astropy.units.Quantity(~numpy.ndarray(float))):
Positions vector in heliocentric equatorial frame in units of AU. nx3
currentTime (~astropy.time.Time):
Current absolute mission time in MJD
rotsign (int):
Optional flag, default 1, set -1 to reverse the rotation
Returns:
~astropy.units.Quantity(~numpy.ndarray(float)):
Positions vector in heliocentric ecliptic frame in units of AU. nx3
"""
# check size of arrays
assert currentTime.size == 1 or currentTime.size == len(
r_equat
), "If multiple times and positions, currentTime and r_equat sizes must match"
# find Julian centuries from J2000
TDB = self.cent(currentTime)
# find obliquity of the ecliptic
obe = rotsign * np.array(np.radians(self.obe(TDB)), ndmin=1)
# positions vector in heliocentric ecliptic frame
_r_equat = r_equat.to_value(u.AU)
if currentTime.size == 1:
r_eclip = (
np.array(
[
np.dot(self.rot(obe[0], 1), _r_equat[x, :])
for x in range(len(_r_equat))
]
)
<< u.AU
)
else:
rots = {_obe: self.rot(_obe, 1) for _obe in np.unique(obe)}
r_eclip = (
np.array(
[np.dot(rots[obe[x]], _r_equat[x, :]) for x in range(len(_r_equat))]
)
<< u.AU
)
return r_eclip
[docs]
def eclip2equat(self, r_eclip, currentTime):
"""Rotates heliocentric coordinates from ecliptic to equatorial frame.
Args:
r_eclip (~astropy.units.Quantity(~numpy.ndarray(float))):
Positions vector in heliocentric ecliptic frame in units of AU
currentTime (astropy.time.Time):
Current absolute mission time in MJD
Returns:
~astropy.units.Quantity(~numpy.ndarray(float)):
Positions vector in heliocentric equatorial frame in units of AU. nx3
"""
r_equat = self.equat2eclip(r_eclip, currentTime, rotsign=-1)
return r_equat
[docs]
def orbit(self, currentTime, eclip=False):
"""Finds observatory orbit positions vector in heliocentric equatorial (default)
or ecliptic frame for current time (MJD).
This method returns the telescope geosynchronous circular orbit position vector.
Args:
currentTime (~astropy.time.Time):
Current absolute mission time in MJD
eclip (bool):
Boolean used to switch to heliocentric ecliptic frame. Defaults to
False, corresponding to heliocentric equatorial frame.
Returns:
~astropy.units.Quantity(~numpy.ndarray(float)):
Observatory orbit positions vector in heliocentric equatorial (default)
or ecliptic frame in units of AU. nx3
Note:
Use eclip=True to get ecliptic coordinates.
"""
mjdtime = np.array(currentTime.mjd, ndmin=1) # modified julian day time
t = mjdtime % 1 # gives percent of day
f = self.orbital_frequency # orbital frequency (2*pi/sideral day)
r = self.orbital_height # orbital height in AU
I = self.orbital_inclination # noqa: E741 # orbital inclination in degrees
O = (
self.orbital_right_ascension
) # noqa: E741 # right ascension of the ascending node
# observatory positions vector wrt Earth in orbital plane
_ft = f * t
r_orb = r * np.array([np.cos(_ft), np.sin(_ft), np.zeros(t.size)])
# observatory positions vector wrt Earth in equatorial frame
r_obs_earth = np.dot(np.dot(self.rot(-O, 3), self.rot(-I, 1)), r_orb).T
# Earth positions vector in heliocentric equatorial frame
r_earth = self.get_Earth_position(currentTime.mjd)
# observatory positions vector in heliocentric equatorial frame
r_obs = (r_obs_earth + r_earth) << u.AU
assert np.all(
np.isfinite(r_obs)
), "Observatory positions vector r_obs has infinite value."
if eclip:
# observatory positions vector in heliocentric ecliptic frame
r_obs = self.equat2eclip(r_obs, currentTime)
return r_obs
[docs]
def keepout(self, TL, sInds, currentTime, koangles, returnExtra=False):
"""Finds keepout Boolean values for stars of interest.
This method returns the keepout Boolean values for stars of interest, where
True is an observable star.
Args:
TL (:ref:`TargetList`):
TargetList class object
sInds (~numpy.ndarray(int)):
Integer indices of the stars of interest
currentTime (~astropy.time.Time):
Current absolute mission time in MJD
MAY ONLY BE ONE VALUE OR DUPLICATES OF THE SAME VALUE
koangles (~astropy.units.Quantity(~numpy.ndarray(float))):
s x 4 x 2 array where s is the number of starlight suppression systems
as defined in the Optical System. Each of the remaining 4 x 2 arrays
are system specific koAngles for the Sun, Moon, Earth, and small bodies
(4), each with a minimum and maximum value (2) in units of deg.
returnExtra (bool):
Optional flag, default False, set True to return additional information.
Returns:
tuple or ~numpy.ndarray(bool):
kogood (~numpy.ndarray(bool)):
kogood s x n x m array of boolean values. True is a target
unobstructed and observable, and False is a target unobservable due
to obstructions in the keepout zone.
r_body (~astropy.units.Quantity(~numpy.ndarray(float))):
Only returned if returnExtra is True
11 x m x 3 array where m is len(currentTime) of heliocentric
equatorial Cartesian elements of the Sun, Moon, Earth and
Mercury->Pluto
r_targ (~astropy.units.Quantity(~numpy.ndarray(float))):
Only returned if returnExtra is True
m x n x 3 array where m is len(currentTime) or 1 if staticStars is
true in TargetList of heliocentric equatorial Cartesian coords of
target and n is the len(sInds)
culprit (numpy.ndarray(int)):
Only returned if returnExtra is True
s x n x m x 12 array of boolean integer values identifying which
body is responsible for keepout (when equal to 1). m is number of
targets and n is len(currentTime). Last dimension is ordered same
as r_body, with an extra line for solar panels being the culprit
koangleArray (~astropy.units.Quantity(~numpy.ndarray(float))):
Only returned if returnExtra is True
s x 11 x 2 element array of minimum and maximum keepouts used for
each body. Same ordering as r_body.
"""
# if multiple time values, check they are different otherwise reduce to scalar
if currentTime.size > 1:
if np.all(currentTime == currentTime[0]):
currentTime = currentTime[0]
# cast sInds to array
sInds = np.array(sInds, ndmin=1, copy=copy_if_needed)
# get all array sizes
nStars = sInds.size
nTimes = currentTime.size
nSystems = koangles.shape[0]
nBodies = 11
# observatory positions vector in heliocentric equatorial frame
r_obs = self.orbit(currentTime) # (m x 3)
# traget star positions vector in heliocentric equatorial frame
r_targ = TL.starprop(sInds, currentTime)
# r_targ = r_targ.reshape(nTimes,nStars,3) # (m x n x 3).
if TL.staticStars and (nStars == 1):
# When the stars are not moving the position vectors are always the same
# so they are tiled because r_targ returns only a 1xnStars array
r_targ = np.tile(r_targ, (nTimes, 1))
else:
r_targ = r_targ.reshape(nTimes, nStars, 3) # (m x n x 3).
# body positions vector in heliocentric equatorial frame
r_body = (
np.array(
[
self.solarSystem_body_position(currentTime, "Sun").to_value("AU"),
self.solarSystem_body_position(currentTime, "Moon").to_value("AU"),
self.solarSystem_body_position(currentTime, "Earth").to_value("AU"),
self.solarSystem_body_position(currentTime, "Mercury").to_value(
"AU"
),
self.solarSystem_body_position(currentTime, "Venus").to_value("AU"),
self.solarSystem_body_position(currentTime, "Mars").to_value("AU"),
self.solarSystem_body_position(currentTime, "Jupiter").to_value(
"AU"
),
self.solarSystem_body_position(currentTime, "Saturn").to_value(
"AU"
),
self.solarSystem_body_position(currentTime, "Uranus").to_value(
"AU"
),
self.solarSystem_body_position(currentTime, "Neptune").to_value(
"AU"
),
self.solarSystem_body_position(currentTime, "Pluto").to_value("AU"),
]
)
<< u.AU
)
# position vectors wrt spacecraft
r_targ = (r_targ - r_obs.reshape(nTimes, 1, 3)).to("pc") # (m x n x 3)
r_body = (r_body - r_obs.reshape(1, nTimes, 3)).to("AU") # (11 x m x 3)
# unit vectors wrt spacecraft
u_targ = (r_targ / np.linalg.norm(r_targ, axis=-1, keepdims=True)).value
u_body = (r_body / np.linalg.norm(r_body, axis=-1, keepdims=True)).value
# create array of koangles for all bodies, using minimum and maximum keepout
# angles of each starlight suppression system in the telescope for
# bright objects (Sun, Moon, Earth, other small bodies)
koangleArray = np.zeros([nSystems, nBodies, 2])
koangleArray[:, 0:3, :] = koangles[:, 0:3, :]
koangleArray[:, 3:, :] = koangles[:, 3, :].reshape(
nSystems, 1, 2
) # small bodies have same values
koangleArray = np.deg2rad(koangleArray)
# find angles and make angle comparisons to build kogood array:
# if bright objects have an angle with the target vector less than koangle
# (e.g. pi/4) they are in the field of view and the target star may not be
# observed, thus ko associated with this target becomes False.
koAngles_SolarPanel = self.koAngles_SolarPanel.to_value(u.rad)
kogood = np.ones([nSystems, nStars, nTimes], dtype=bool) # (s x n x m)
culprit = np.zeros([nSystems, nStars, nTimes, nBodies + 1]) # (s x n x m x 12)
# Compute keepout for each starlight suppression system
for s in tqdm(
np.arange(nSystems), desc="Starlight Suppression System", position=0
):
# Get keepout constraints for this system once
# Lower bound
ko_lower = koangleArray[s, :, 0]
# Upper bound
ko_upper = koangleArray[s, :, 1]
# Compute keepout for each star
for n in tqdm(
np.arange(nStars), desc="Star Keepout", position=1, leave=False
):
# unit vectors for the 11 bodies and the nth target at all times
u_t_all = u_targ[:, n, :]
# Dot product for all bodies at all times
# u_body shape: (11, nTimes, 3), u_t_all shape: (nTimes, 3)
# Reshape for broadcasting: u_t_all to (1, nTimes, 3)
dot_products = np.sum(u_body * u_t_all.reshape(1, nTimes, 3), axis=2).T
# dot_products shape: (nTimes, 11)
# relative angle between the target and bright body look vectors
angles = np.arccos(np.clip(dot_products, -1, 1))
# create array of "culprits" that prevent a target from being
# observed
culprit[s, n, :, :-1] = (angles < ko_lower) | (angles > ko_upper)
# adding solar panel restrictions as a final culprit
culprit[s, n, :, -1] = (angles[:, 0] < koAngles_SolarPanel[0]) | (
angles[:, 0] > koAngles_SolarPanel[1]
)
# if any bright body obstructs, kogood becomes False
kogood[s, n, :] = ~np.any(culprit[s, n, :, :], axis=1)
# checking that all ko elements are boolean
trues = [isinstance(element, np.bool_) for element in kogood.flatten()]
assert all(trues), "An element of kogood is not Boolean"
if returnExtra:
return kogood, r_body, r_targ, culprit, koangleArray
else:
return kogood
[docs]
def generate_koMap(self, TL, missionStart, missionFinishAbs, koangles):
"""Creates keepout map for all targets throughout mission lifetime.
This method returns a binary map showing when all stars in the given
target list are in or out of the keepout zone (i.e. when they are not
observable) from mission start to mission finish.
Args:
TL (:ref:`TargetList`):
TargetList class object
missionStart (~astropy.time.Time):
Absolute start of mission time in MJD
missionFinishAbs (~astropy.time.Time):
Absolute end of mission time in MJD
koangles (~astropy.units.Quantity(~numpy.ndarray(float))):
s x 4 x 2 array where s is the number of starlight suppression systems
as defined in the Optical System. Each of the remaining 4 x 2 arrays
are system specific koAngles for the Sun, Moon, Earth, and small bodies
(4), each with a minimum and maximum value (2) in units of deg.
Returns:
tuple:
koMap (~numpy.ndarray(bool)):
True is a target unobstructed and observable, and False is a
target unobservable due to obstructions in the keepout zone.
koTimes (~astropy.time.Time):
Absolute MJD mission times from start to end in steps of 1 d
"""
# generating hash name
filename = "koMap_"
atts = ["koAngles_SolarPanel", "ko_dtStep"]
extstr = ""
for att in sorted(atts, key=str.lower):
if not callable(getattr(self, att)):
extstr += "%s: " % att + str(getattr(self, att)) + " "
extstr += "%s: " % "missionStart" + str(missionStart) + " "
extstr += "%s: " % "missionFinishAbs" + str(missionFinishAbs) + " "
extstr += "%s: " % "koangles" + str(koangles) + " "
extstr += "%s: " % "Name" + str(getattr(TL, "Name")) + " "
# TODO: is this needed?
# extstr += '%s: ' % 'Name' + TL.StarCatalog.__class__.__name__ + ' '
extstr += "%s: " % "nStars" + str(getattr(TL, "nStars")) + " "
ext = hashlib.md5(extstr.encode("utf-8")).hexdigest()
filename += ext
koPath = os.path.join(self.cachedir, filename + ".komap")
# global times when keepout is checked for all stars
koTimes = np.arange(
missionStart.value, missionFinishAbs.value, self.ko_dtStep.value
)
koTimes = Time(
koTimes, format="mjd", scale="tai"
) # scale must be tai to account for leap seconds
if os.path.exists(koPath):
# keepout map already exists for parameters
self.vprint("Loading cached keepout map file from %s" % koPath)
try:
with open(koPath, "rb") as ff:
A = pickle.load(ff)
except UnicodeDecodeError:
with open(koPath, "rb") as ff:
A = pickle.load(ff, encoding="latin1")
self.vprint("Keepout Map loaded from cache.")
koMap = A["koMap"]
else:
self.vprint('Cached keepout map file not found at "%s".' % koPath)
# looping over all stars to generate map of when all stars are observable
self.vprint("Starting keepout calculations for %s stars." % TL.nStars)
koMap = self.keepout(TL, np.arange(TL.nStars), koTimes, koangles, False)
A = {"koMap": koMap}
with open(koPath, "wb") as f:
pickle.dump(A, f)
self.vprint("Keepout Map calculations finished")
self.vprint("Keepout Map array stored in %r" % koPath)
return koMap, koTimes
[docs]
def calculate_observableTimes(self, TL, sInds, currentTime, koMaps, koTimes, mode):
"""Returns the next window of time during which targets are observable
This method returns a nx2 ndarray of times for every star given in the
target list. The two entries for every star are the next times (after
current time) when the star exits and enters keepout (i.e. the start
and end times of the next window of observability).
Args:
TL (:ref:`TargetList`):
TargetList class object
sInds (numpy.ndarray(int)):
Integer indices of the stars of interest
currentTime (~astropy.time.Time):
Current absolute mission time in MJD
koMaps (dict):
Keepout values for n stars throughout time range of length m,
key names being the system names specified in mode.
True is a target unobstructed and observable, and False is a
target unobservable due to obstructions in the keepout zone.
koTimes (~astropy.time.Time):
Absolute MJD mission times from start to end in steps of 1 d
mode (dict):
Selected observing mode
Returns:
~astropy.time.Time:
Start and end times of next observability time window in
absolute time MJD. n is length of sInds
"""
# creating time arrays to use in the keepout method (# stars == # times)
# minimum slew time for occulter to align with new star
if mode["syst"]["occulter"]:
nextObTimes = np.ones(len(sInds)) * currentTime.value + self.occ_dtmin.value
else:
nextObTimes = np.ones(len(sInds)) * currentTime.value
nextObTimes = Time(
nextObTimes, format="mjd", scale="tai"
) # converting to astropy MJD time array
# find appropriate koMap
systName = mode["syst"]["name"]
koMap = koMaps[systName]
# finding observable times
observableTimes = self.find_nextObsWindow(
TL, sInds, nextObTimes, koMap, koTimes
).value
observableTimesNorm = (
observableTimes - nextObTimes.value
) # days since currentTime
# in case of an occulter, correct for short windows
if mode["syst"]["occulter"]:
# find length of observable range in days
observable_range = np.diff(observableTimesNorm, axis=0)[0]
# re-do calculations for observable windows that are less than
# dt_min days long
reDo = np.where(observable_range < self.occ_dtmin.value)[0]
if reDo.size:
correctedObTimes = (
nextObTimes[reDo].value + observableTimesNorm[1, reDo]
)
correctedObTimes = Time(correctedObTimes, format="mjd", scale="tai")
observableTimes[:, reDo] = self.find_nextObsWindow(
TL, reDo, correctedObTimes, koMap, koTimes
).value
return Time(observableTimes, format="mjd", scale="tai")
[docs]
def find_nextObsWindow(self, TL, sInds, currentTimes, koMap, koTimes):
"""Method used by calculate_observableTimes for finding next observable windows
This method returns a nx2 ndarray of times for every star given in the
target list. The two entries for every star are the next times (after
current time) when the star exits and enters keepout (i.e. the start
and end times of the next window of observability).
Args:
TL (:ref:`TargetList`):
TargetList class object
sInds (~numpy.ndarray(int)):
Integer indices of the stars of interest
currentTimes (~astropy.time.Time):
Current absolute mission time in MJD same length as sInds
koMap (~numpy.ndarray(int)):
Keepout values for n stars throughout time range of length m (mxn)
koTimes (astropy.time.Time):
Absolute MJD mission times from start to end in steps of 1 d
mode (dict):
Selected observing mode
Returns:
astropy.time.Time(~numpy.ndarray):
Start and end times of next observability time window in MJD
"""
# create arrays
nLoops = len(sInds)
nextExitTime = np.zeros(nLoops)
nextEntryTime = np.zeros(nLoops)
# getting saved time closest to currentTime
currentTimes_mjd = currentTimes.mjd
koTimes_mjd = koTimes.mjd
# shape: (len(koTimes), nLoops)
time_diffs = np.abs(koTimes_mjd[:, np.newaxis] - currentTimes_mjd)
T = np.argmin(time_diffs, axis=0)
# checking to see if stars are in keepout at currentTime
kogoodStart = [bool(koMap[x, S]) for x, S in zip(sInds, T)]
kobadStart = [bool(not koMap[x, S]) for x, S in zip(sInds, T)]
nextExitTime[kogoodStart] = currentTimes_mjd[kogoodStart]
# finding next entry into keepout for currently observable stars
for n, S in zip(sInds[kogoodStart], T[kogoodStart]):
idxG_E = np.where(~koMap[n, S:])
# enters KO after missionEnd
if not idxG_E[0].tolist():
nEnd = np.where(sInds == n)
nextEntryTime[nEnd] = koTimes_mjd[-1]
else:
nextEntry = idxG_E[0][0] + S
# enters KO after missionEnd (missed these)
if nextEntry > len(koTimes):
nEnd = np.where(sInds == n)
nextEntryTime[nEnd] = koTimes_mjd[-1]
# enters KO before missionEnd
else:
nGood = np.where(sInds == n)
nextEntryTime[nGood] = koTimes_mjd[nextEntry]
# finding next exit and entry of keepout for unobservable stars (in keepout)
for n, S in zip(sInds[kobadStart], T[kobadStart]):
idx_X = np.where(koMap[n, S:])
# exit KO after missionEnd (enter after as well)
if not idx_X[0].tolist():
nEnd = np.where(sInds == n)
nextExitTime[nEnd] = koTimes_mjd[-1]
nextEntryTime[nEnd] = koTimes_mjd[-1]
else:
nextExit = idx_X[0][0] + S
# exit KO after missionEnd (missed these)
if nextExit > len(koTimes):
nEnd = np.where(sInds == n)
nextExitTime[nEnd] = koTimes_mjd[-1]
nextEntryTime[nEnd] = koTimes_mjd[-1]
# exit KO before missionEnd
else:
nBad = np.where(sInds == n)
nextExitTime[nBad] = koTimes_mjd[nextExit]
idx_E = np.where(~koMap[n, nextExit:])
# enters KO again after missionEnd
if not idx_E[0].tolist():
nEnd = np.where(sInds == n)
nextEntryTime[nEnd] = koTimes_mjd[-1]
else:
nextEntry = idx_E[0][0] + nextExit
# enters KO again after missionEnd (missed these)
if nextEntry > len(koTimes):
nEnd = np.where(sInds == n)
nextEntryTime[nEnd] = koTimes_mjd[-1]
# enters KO before missionEnd
else:
nextEntryTime[nBad] = koTimes_mjd[nextEntry]
observableTimes = np.vstack([nextExitTime, nextEntryTime]) << u.d
return observableTimes
[docs]
def star_angularSep(self, TL, old_sInd, sInds, currentTime):
"""Finds angular separation from old star to given list of stars
This method returns the angular separation from the last observed
star to all others on the given list at the currentTime.
Args:
TL (:ref:`TargetList`):
TargetList class object
old_sInd (int):
Integer index of the last star of interest
sInds (~numpy.ndarray(int)):
Integer indices of the stars of interest
currentTime (~astropy.time.Time):
Current absolute mission time in MJD
Returns:
float:
Angular separation between two target stars
"""
if old_sInd is None:
sd = np.zeros(len(sInds)) << u.rad
else:
# position vector of previous target star
r_old = TL.starprop(old_sInd, currentTime)[0]
u_old = r_old.to_value("AU") / np.linalg.norm(r_old.to_value("AU"))
# position vector of new target stars
r_new = TL.starprop(sInds, currentTime)
u_new = (
r_new.to_value("AU").T / np.linalg.norm(r_new.to_value("AU"), axis=1)
).T
# angle between old and new stars
sd = np.arccos(np.clip(np.dot(u_old, u_new.T), -1, 1)) << u.rad
# A-frame
a1 = u_old / np.linalg.norm(u_old) # normalized old look vector
a2 = np.array([a1[1], -a1[0], 0]) # normal to a1
a3 = np.cross(a1, a2) # last part of the A basis vectors
# finding sign of angle
# The star angular separation can be negative
u2_Az = np.dot(a3, u_new.T)
sgn = np.sign(u2_Az)
sgn[np.where(sgn == 0)] = 1
sd = sgn * sd
return sd
[docs]
def solarSystem_body_position(self, currentTime, bodyname, eclip=False):
"""Finds solar system body positions vector in heliocentric equatorial (default)
or ecliptic frame for current time (MJD).
This passes all arguments to one of spk_body or keplerplanet, depending
on the value of self.havejplephem.
Args:
currentTime (~astropy.time.Time):
Current absolute mission time in MJD
bodyname (str):
Solar system object name
eclip (bool):
Boolean used to switch to heliocentric ecliptic frame. Defaults to
False, corresponding to heliocentric equatorial frame.
Returns:
~astropy.units.Quantity(~numpy.ndarray(float)):
Solar system body positions in heliocentric equatorial (default)
or ecliptic frame in units of AU. nx3
Note:
Use eclip=True to get ecliptic coordinates.
"""
# heliocentric
if bodyname == "Sun":
return np.zeros((currentTime.size, 3)) << u.AU
# choose JPL or static ephemerides
if self.havejplephem:
r_body = self.spk_body(currentTime, bodyname, eclip=eclip).to("AU")
else:
r_body = self.keplerplanet(currentTime, bodyname, eclip=eclip).to("AU")
return r_body
[docs]
def spk_body(self, currentTime, bodyname, eclip=False):
"""Finds solar system body positions vector in heliocentric equatorial (default)
or ecliptic frame for current time (MJD).
This method uses spice kernel from NAIF to find heliocentric
equatorial position vectors for solar system objects.
Args:
currentTime (~astropy.time.Time):
Current absolute mission time in MJD
bodyname (str):
Solar system object name
eclip (bool):
Boolean used to switch to heliocentric ecliptic frame. Defaults to
False, corresponding to heliocentric equatorial frame.
Returns:
~astropy.units.Quantity(~numpy.ndarray(float)):
Solar system body positions in heliocentric equatorial (default)
or ecliptic frame in units of AU. nx3
Note: Use eclip=True to get ecliptic coordinates.
"""
# dictionary of solar system bodies available in spice kernel (in km)
bodies = {
"Mercury": 199,
"Venus": 299,
"Earth": 399,
"Mars": 4,
"Jupiter": 5,
"Saturn": 6,
"Uranus": 7,
"Neptune": 8,
"Pluto": 9,
"Sun": 10,
"Moon": 301,
}
assert bodyname in bodies, "%s is not a recognized body name." % (bodyname)
# julian day time
jdtime = np.array(currentTime.jd, ndmin=1)
# body positions vector in heliocentric equatorial frame
if bodies[bodyname] == 199:
r_body = (
self.kernel[0, 1].compute(jdtime)
+ self.kernel[1, 199].compute(jdtime)
- self.kernel[0, 10].compute(jdtime)
)
elif bodies[bodyname] == 299:
r_body = (
self.kernel[0, 2].compute(jdtime)
+ self.kernel[2, 299].compute(jdtime)
- self.kernel[0, 10].compute(jdtime)
)
elif bodies[bodyname] == 399:
r_body = (
self.kernel[0, 3].compute(jdtime)
+ self.kernel[3, 399].compute(jdtime)
- self.kernel[0, 10].compute(jdtime)
)
elif bodies[bodyname] == 301:
r_body = (
self.kernel[0, 3].compute(jdtime)
+ self.kernel[3, 301].compute(jdtime)
- self.kernel[0, 10].compute(jdtime)
)
else:
r_body = self.kernel[0, bodies[bodyname]].compute(jdtime) - self.kernel[
0, 10
].compute(jdtime)
# reshape and convert units
# The constant here is the km to AU conversion factor
r_body = r_body.T * self.km2au << u.AU
if eclip:
# body positions vector in heliocentric ecliptic frame
r_body = self.equat2eclip(r_body, currentTime)
return r_body
[docs]
def generate_Earth_interpolator(self, startTime, endTime, dt):
"""Creates an interpolator for Earth position over the full mission and
to save time during orbit calculations.
Args:
startTime (~astropy.time.Time):
Start time of the time range
endTime (~astropy.time.Time):
End time of the time range
dt (float):
Time step in days
"""
times = Time(np.arange(startTime.mjd, endTime.mjd, dt), format="mjd")
times_mjd = times.to_value("mjd")
# Get the Earth position in heliocentric ecliptic frame
earth_pos_eclip_path = Path(
self.cachedir,
f"Earth_positions_eclip_{startTime.mjd}_{endTime.mjd}_{dt}.npy",
)
if earth_pos_eclip_path.exists():
earth_pos_eclip = np.load(earth_pos_eclip_path)
else:
earth_pos_eclip = self.solarSystem_body_position(
times, "Earth", eclip=True
).to_value("AU")
np.save(earth_pos_eclip_path, earth_pos_eclip)
self.earth_pos_interp = {}
self.earth_pos_interp[True] = CubicSpline(
times_mjd, earth_pos_eclip, extrapolate=False
)
# Get the Earth position in heliocentric equatorial frame
earth_pos_equat_path = Path(
self.cachedir,
f"Earth_positions_equat_{startTime.mjd}_{endTime.mjd}_{dt}.npy",
)
if earth_pos_equat_path.exists():
earth_pos_equat = np.load(earth_pos_equat_path)
else:
earth_pos_equat = self.solarSystem_body_position(
times, "Earth", eclip=False
).to_value("AU")
np.save(earth_pos_equat_path, earth_pos_equat)
self.earth_pos_interp[False] = CubicSpline(
times_mjd, earth_pos_equat, extrapolate=False
)
self.has_earth_pos_interp = True
self.earth_pos_interp_range_mjd = [startTime.mjd, endTime.mjd]
[docs]
def get_Earth_position(self, currentTime_mjd, eclip=False):
"""Get the Earth position in a given heliocentric frame at a given time.
Preferentially uses the CubicSpline interpolant but will calculate it
directly as a fall-back. This was set up for test cases where the SurveySim
module is not initialized. Default is heliocentric equatorial frame.
Args:
currentTime_mjd (float):
Current absolute mission time in MJD
eclip (bool):
Boolean used to switch to heliocentric ecliptic frame. Defaults to
False, corresponding to heliocentric equatorial frame.
Returns:
numpy.ndarray(float):
Earth position in heliocentric frame in units of AU
"""
if self.has_earth_pos_interp and (
np.all(currentTime_mjd >= self.earth_pos_interp_range_mjd[0])
and np.all(currentTime_mjd <= self.earth_pos_interp_range_mjd[1])
):
# Use the interpolator if it exists
return self.earth_pos_interp[eclip](currentTime_mjd).reshape(-1, 3)
else:
# Otherwise, calculate the position directly
return self.solarSystem_body_position(
Time(currentTime_mjd, format="mjd", scale="tai"), "Earth", eclip=eclip
).to_value("AU")
[docs]
def keplerplanet(self, currentTime, bodyname, eclip=False):
"""Finds solar system body positions vector in heliocentric equatorial (default)
or ecliptic frame for current time (MJD).
This method uses algorithms 2 and 10 from Vallado 2013 to find
heliocentric equatorial position vectors for solar system objects.
Args:
currentTime (~astropy.time.Time):
Current absolute mission time in MJD
bodyname (str):
Solar system object name
eclip (bool):
Boolean used to switch to heliocentric ecliptic frame. Defaults to
False, corresponding to heliocentric equatorial frame.
Returns:
~astropy.units.Quantity(~numpy.ndarray(float)):
Solar system body positions in heliocentric equatorial (default)
or ecliptic frame in units of AU
Note:
Use eclip=True to get ecliptic coordinates.
"""
# Moon positions based on Earth positions
if bodyname == "Moon":
r_Earth = self.keplerplanet(currentTime, "Earth")
return r_Earth + self.moon_earth(currentTime)
assert bodyname in self.planets, "%s is not a recognized body name." % (
bodyname
)
# find Julian centuries from J2000
TDB = self.cent(currentTime)
# update ephemerides data (convert sma from km to AU)
planet = self.planets[bodyname]
a = self.propeph(planet.a, TDB) * self.km2au << u.AU
e = self.propeph(planet.e, TDB)
I = np.radians(self.propeph(planet.I, TDB)) # noqa: E741
O = np.radians(self.propeph(planet.O, TDB)) # noqa: E741
w = np.radians(self.propeph(planet.w, TDB))
lM = np.radians(self.propeph(planet.lM, TDB))
# find mean anomaly and argument of perigee
M = (lM - w) % (2.0 * np.pi)
wp = (w - O) % (2.0 * np.pi)
# find eccentric anomaly
E = eccanom(M, e)[0]
# find true anomaly
nu = np.arctan2(np.sin(E) * np.sqrt(1.0 - e**2.0), np.cos(E) - e)
# find semiparameter
p = a * (1.0 - e**2.0)
# body positions vector in orbital plane
rx = p * np.cos(nu) / (1.0 + e * np.cos(nu))
ry = p * np.sin(nu) / (1.0 + e * np.cos(nu))
rz = np.zeros(currentTime.size)
r_orb = np.array([rx, ry, rz])
# body positions vector in heliocentric ecliptic plane
r_body = (
np.array(
[
np.dot(
np.dot(self.rot(-O[x], 3), self.rot(-I[x], 1)),
np.dot(self.rot(-wp[x], 3), r_orb[:, x]),
)
for x in range(currentTime.size)
]
)
<< u.AU
)
if not eclip:
# body positions vector in heliocentric equatorial frame
r_body = self.eclip2equat(r_body, currentTime)
return r_body
[docs]
def moon_earth(self, currentTime):
"""Finds geocentric equatorial positions vector for Earth's moon
This method uses Algorithm 31 from Vallado 2013 to find the geocentric
equatorial positions vector for Earth's moon.
Args:
currentTime (~astropy.time.Time):
Current absolute mission time in MJD
Returns:
~astropy.units.Quantity(~numpy.ndarray(float)):
Geocentric equatorial position vector in units of AU
"""
TDB = np.array(self.cent(currentTime), ndmin=1)
la = np.radians(
218.32
+ 481267.8813 * TDB
+ 6.29 * np.sin(np.radians(134.9 + 477198.85 * TDB))
- 1.27 * np.sin(np.radians(259.2 - 413335.38 * TDB))
+ 0.66 * np.sin(np.radians(235.7 + 890534.23 * TDB))
+ 0.21 * np.sin(np.radians(269.9 + 954397.70 * TDB))
- 0.19 * np.sin(np.radians(357.5 + 35999.05 * TDB))
- 0.11 * np.sin(np.radians(186.6 + 966404.05 * TDB))
)
phi = np.radians(
5.13 * np.sin(np.radians(93.3 + 483202.03 * TDB))
+ 0.28 * np.sin(np.radians(228.2 + 960400.87 * TDB))
- 0.28 * np.sin(np.radians(318.3 + 6003.18 * TDB))
- 0.17 * np.sin(np.radians(217.6 - 407332.20 * TDB))
)
P = np.radians(
0.9508
+ 0.0518 * np.cos(np.radians(134.9 + 477198.85 * TDB))
+ 0.0095 * np.cos(np.radians(259.2 - 413335.38 * TDB))
+ 0.0078 * np.cos(np.radians(235.7 + 890534.23 * TDB))
+ 0.0028 * np.cos(np.radians(269.9 + 954397.70 * TDB))
)
e = np.radians(
23.439291 - 0.0130042 * TDB - 1.64e-7 * TDB**2 + 5.04e-7 * TDB**3
)
r = 1.0 / np.sin(P) * 6378.137 # km
r_moon = r * np.array(
[
np.cos(phi) * np.cos(la),
np.cos(e) * np.cos(phi) * np.sin(la) - np.sin(e) * np.sin(phi),
np.sin(e) * np.cos(phi) * np.sin(la) + np.cos(e) * np.sin(phi),
]
)
# set format and units
r_moon = (r_moon * u.km).T.to("AU")
return r_moon
[docs]
def cent(self, currentTime):
"""Finds time in Julian centuries since J2000 epoch
This quantity is needed for many algorithms from Vallado 2013.
Args:
currentTime (~astropy.time.Time):
Current absolute mission time in MJD
Returns:
~numpy.ndarray(float):
time in Julian centuries since the J2000 epoch
"""
TDB = (currentTime.jd - self.j2000_jd) / self.julian_century
return TDB
[docs]
def propeph(self, x, TDB):
"""Propagates an ephemeris from Vallado 2013 to current time.
Args:
x (list):
ephemeride list (maximum of 4 elements)
TDB (float):
time in Julian centuries since the J2000 epoch
Returns:
numpy.darray(float):
ephemerides value at current time
"""
if isinstance(x, list):
if len(x) < 4:
q = 4 - len(x)
i = 0
while i < q:
x.append(0.0)
i += 1
elif isinstance(x, float) or isinstance(x, int):
x = [float(x)]
i = 0
while i < 3:
x.append(0.0)
i += 1
# propagated ephem
y = x[0] + x[1] * TDB + x[2] * (TDB**2) + x[3] * (TDB**3)
# cast to array
y = np.array(y, ndmin=1, copy=copy_if_needed)
return y
[docs]
def rot(self, th, axis):
"""Finds the rotation matrix of angle th about the axis value
Args:
th (float):
Rotation angle in radians
axis (int):
Integer value denoting rotation axis (1,2, or 3)
Returns:
~numpy.ndarray(float):
Rotation matrix
"""
_sin = np.sin(th)
_cos = np.cos(th)
if axis == 1:
rot_th = np.array(
[
[1.0, 0.0, 0.0],
[0.0, _cos, _sin],
[0.0, -_sin, _cos],
]
)
elif axis == 2:
rot_th = np.array(
[
[_cos, 0.0, -_sin],
[0.0, 1.0, 0.0],
[_sin, 0.0, _cos],
]
)
elif axis == 3:
rot_th = np.array(
[
[_cos, _sin, 0.0],
[-_sin, _cos, 0.0],
[0.0, 0.0, 1.0],
]
)
return rot_th
[docs]
def distForces(self, TL, sInd, currentTime):
"""Finds lateral and axial disturbance forces on an occulter
Args:
TL (:ref:`TargetList`):
TargetList class object
sInd (int):
Integer index of the star of interest
currentTime (~astropy.time.Time):
Current absolute mission time in MJD
Returns:
tuple:
:obj:`~astropy.units.Quantity`:
dF_lateral: Lateral disturbance force in units of N
:obj:`~astropy.units.Quantity`:
dF_axial: Axial disturbance force in units of N
"""
# get spacecraft position vector
r_obs = self.orbit(currentTime)[0]
# sun -> earth position vector
r_Es = self.solarSystem_body_position(currentTime, "Earth")[0]
# Telescope -> target vector and unit vector
r_targ = TL.starprop(sInd, currentTime)[0] - r_obs
u_targ = r_targ.to("AU").value / np.linalg.norm(r_targ.to("AU").value)
# sun -> occulter vector
r_Os = r_obs.to("AU") + self.occulterSep.to("AU") * u_targ
# Earth-Moon barycenter -> spacecraft vectors
r_TE = r_obs - r_Es
r_OE = r_Os - r_Es
# force on occulter
Mfactor = -self.scMass * const.M_sun * const.G
F_sO = r_Os / (np.linalg.norm(r_Os.to("AU").value) * r_Os.unit) ** 3.0 * Mfactor
F_EO = (
r_OE
/ (np.linalg.norm(r_OE.to("AU").value) * r_OE.unit) ** 3.0
* Mfactor
/ 328900.56
)
F_O = F_sO + F_EO
# force on telescope
Mfactor = -self.coMass * const.M_sun * const.G
F_sT = (
r_obs / (np.linalg.norm(r_obs.to("AU").value) * r_obs.unit) ** 3.0 * Mfactor
)
F_ET = (
r_TE
/ (np.linalg.norm(r_TE.to("AU").value) * r_TE.unit) ** 3.0
* Mfactor
/ 328900.56
)
F_T = F_sT + F_ET
# differential forces
dF = F_O - F_T * self.scMass / self.coMass
dF_axial = (dF.dot(u_targ)).to("N")
dF_lateral = (dF - dF_axial * u_targ).to("N")
dF_lateral = np.linalg.norm(dF_lateral.to("N").value) * dF_lateral.unit
dF_axial = np.abs(dF_axial)
return dF_lateral, dF_axial
[docs]
def mass_dec(self, dF_lateral, t_int):
"""Returns mass_used and deltaV
The values returned by this method are used to decrement spacecraft
mass for station-keeping.
Args:
dF_lateral (astropy.units.Quantity):
Lateral disturbance force in units of N
t_int (astropy.units.Quantity):
Integration time in units of day
Returns:
tuple:
intMdot (astropy.units.Quantity):
Mass flow rate in units of kg/s
mass_used (astropy.units.Quantity):
Mass used in station-keeping units of kg
deltaV (astropy.units.Quantity):
Change in velocity required for station-keeping in units of km/s
"""
intMdot = (dF_lateral / self.skEff / const.g0 / self.skIsp).to("kg/s")
mass_used = (intMdot * t_int).to("kg")
deltaV = (dF_lateral / self.scMass * t_int).to("km/s")
return intMdot, mass_used, deltaV
[docs]
def mass_dec_sk(self, TL, sInd, currentTime, t_int):
"""Returns mass_used, deltaV and disturbance forces
This method calculates all values needed to decrement spacecraft mass
for station-keeping.
Args:
TL (:ref:`TargetList`):
TargetList class object
sInd (int):
Integer index of the star of interest
currentTime (astropy.time.Time):
Current absolute mission time in MJD
t_int (astropy.units.Quantity):
Integration time in units of day
Returns:
tuple:
dF_lateral (astropy.units.Quantity):
Lateral disturbance force in units of N
dF_axial (astropy.units.Quantity):
Axial disturbance force in units of N
intMdot (astropy.units.Quantity):
Mass flow rate in units of kg/s
mass_used (astropy.units.Quantity):
Mass used in station-keeping units of kg
deltaV (astropy.units.Quantity):
Change in velocity required for station-keeping in units of km/s
"""
dF_lateral, dF_axial = self.distForces(TL, sInd, currentTime)
intMdot, mass_used, deltaV = self.mass_dec(dF_lateral, t_int)
return dF_lateral, dF_axial, intMdot, mass_used, deltaV
[docs]
def calculate_dV(self, TL, old_sInd, sInds, sd, slewTimes, tmpCurrentTimeAbs):
"""Finds the change in velocity needed to transfer to a new star line of sight
This method sums the total delta-V needed to transfer from one star
line of sight to another. It determines the change in velocity to move from
one station-keeping orbit to a transfer orbit at the current time, then from
the transfer orbit to the next station-keeping orbit at currentTime + dt.
Station-keeping orbits are modeled as discrete boundary value problems.
This method can handle multiple indeces for the next target stars and calculates
the dVs of each trajectory from the same starting star.
The prototype implementation does not perform any real calculations and returns
all zero values.
Args:
TL (:ref:`TargetList`):
TargetList class object
old_sInd (int):
Index of the current star
sInds (~numpy.ndarray(int)):
Integer index of the next star(s) of interest
sd (~astropy.units.Quantity(~numpy.ndarray(float))):
Angular separation between stars in rad
slewTimes (~astropy.time.Time(~numpy.ndarray)):
Slew times.
tmpCurrentTimeAbs (~astropy.time.Time):
Current absolute mission time in MJD
Returns:
~astropy.units.Quantity(~numpy.ndarray(float)):
Delta-V values in units of length/time
"""
dV = np.zeros(len(sInds))
return dV * u.m / u.s
[docs]
def calculate_slewTimes(self, TL, old_sInd, sInds, sd, obsTimes, currentTime):
"""Finds slew times and separation angles between target stars
This method determines the slew times of an occulter spacecraft needed
to transfer from one star's line of sight to all others in a given
target list.
Args:
TL (:ref:`TargetList`):
TargetList class object
old_sInd (int):
Integer index of the most recently observed star
sInds (~numpy.ndarray(int)):
Integer indices of the star of interest
sd (~astropy.units.Quantity):
Angular separation between stars in rad
obsTimes (~astropy.time.Time(~numpy.ndarray)):
Observation times for targets.
currentTime (~astropy.time.Time(~numpy.ndarray)):
Current absolute mission time in MJD
Returns:
~astropy.units.Quantity:
Time to transfer to new star line of sight in units of days
"""
self.ao = self.thrust / self.scMass
slewTime_fac = (
(
2.0
* self.occulterSep
/ np.abs(self.ao)
/ (self.defburnPortion / 2.0 - self.defburnPortion**2.0 / 4.0)
)
.decompose()
.to("d2")
)
if old_sInd is None:
slewTimes = np.zeros(TL.nStars) * u.d
else:
# calculate slew time
slewTimes = np.sqrt(
slewTime_fac * np.sin(abs(sd) / 2.0)
) # an issue exists if sd is negative
# The following are debugging
assert (
np.where(np.isnan(slewTimes))[0].shape[0] == 0
), "At least one slewTime is nan"
return slewTimes
[docs]
def log_occulterResults(self, DRM, slewTimes, sInd, sd, dV):
"""Updates the given DRM to include occulter values and results
Args:
DRM (dict):
Design Reference Mission, contains the results of one complete
observation (detection and characterization)
slewTimes (astropy.units.Quantity):
Time to transfer to new star line of sight in units of days
sInd (int):
Integer index of the star of interest
sd (astropy.units.Quantity):
Angular separation between stars in rad
dV (astropy.units.Quantity):
Delta-V used to transfer to new star line of sight in units of m/s
Returns:
dict:
Design Reference Mission dictionary, contains the results of one
complete observation (detection and characterization)
"""
DRM["slew_time"] = slewTimes.to("day")
DRM["slew_angle"] = sd.to("deg")
slew_mass_used = slewTimes * self.defburnPortion * self.flowRate
DRM["slew_dV"] = (slewTimes * self.ao * self.defburnPortion).to("m/s")
DRM["slew_mass_used"] = slew_mass_used.to("kg")
self.scMass = self.scMass - slew_mass_used
DRM["scMass"] = self.scMass.to("kg")
if self.twotanks:
self.slewMass = self.slewMass - slew_mass_used
DRM["slewMass"] = self.slewMass.to("kg")
return DRM
[docs]
def refuel_tank(self, TK, tank=None):
"""Attempt to refuel a fuel tank and report status
Args:
TK (:ref:`TimeKeeping`):
TimeKeeping object. Not used in prototype but an input for any
implementations that wish to do time-aware operations.
tank (str, optional):
Either 'sk' or 'slew' when ``twotanks`` is True. Otherwise, None.
Defaults None
Returns:
bool:
True represents successful refeuling. False means refueling is not
possible for selected tank.
"""
if not (self.allowRefueling):
return False
if self.external_fuel_mass <= 0 * u.kg:
return False
if tank is not None:
assert tank.lower() in ["sk", "slew"], "Tank must be 'sk' or 'slew'."
assert self.twotanks, "You may only specify a tank when twotanks is True."
if tank == "sk":
tank_mass = self.skMass
tank_capacity = self.skMaxFuelMass
tank_name = "stationkeeping"
else:
tank_mass = self.slewMass
tank_capacity = self.slewMaxFuelMass
tank_name = "slew"
else:
tank_mass = self.scMass
tank_capacity = self.maxFuelMass + self.dryMass
tank_name = ""
# Add as much fuel as can fit in the tank (plus any currently carried negative
# value, or whatever remains in the external tank)
topoff = (
np.min(
[
self.external_fuel_mass.to(u.kg).value,
(tank_capacity - tank_mass).to(u.kg).value,
]
)
* u.kg
)
assert topoff >= 0 * u.kg, "Topoff calculation produced negative result."
self.external_fuel_mass -= topoff
tank_mass += topoff
if tank is not None:
self.scMass += topoff
self.vprint("{} {} fuel added".format(topoff, tank_name))
self.vprint("{} remaining in external tank.".format(self.external_fuel_mass))
return True
[docs]
class SolarEph:
"""Solar system ephemerides class
This class takes the constants in Appendix D.4 of Vallado as inputs
and stores them for use in defining solar system ephemerides at a
given time.
Args:
a (list):
semimajor axis list (in AU)
e (list):
eccentricity list
I (list):
inclination list
O (list):
right ascension of the ascending node list
w (list):
longitude of periapsis list
lM (list):
mean longitude list
Each of these lists has a maximum of 4 elements. The values in
these lists are used to propagate the solar system planetary
ephemerides for a specific solar system planet.
Attributes:
a (list):
list of semimajor axis (in AU)
e (list):
list of eccentricity
I (list):
list of inclination
O (list):
list of right ascension of the ascending node
w (list):
list of longitude of periapsis
lM (list):
list of mean longitude values
Each of these lists has a maximum of 4 elements. The values in
these lists are used to propagate the solar system planetary
ephemerides for a specific solar system planet.
"""
def __init__(self, a, e, I, O, w, lM): # noqa: E741
# store list of semimajor axis values (convert from AU to km)
self.a = (a * u.AU).to_value(u.km)
if not isinstance(self.a, float):
self.a = self.a.tolist()
# store list of dimensionless eccentricity values
self.e = e
# store list of inclination values (degrees)
self.I = I # noqa: E741
# store list of right ascension of ascending node values (degrees)
self.O = O # noqa: E741
# store list of longitude of periapsis values (degrees)
self.w = w
# store list of mean longitude values (degrees)
self.lM = lM
def __str__(self):
"""String representation of the SolarEph object
When the command 'print' is used on the SolarEph object, this
method will print the attribute values contained in the object
"""
for att in self.__dict__:
print("%s: %r" % (att, getattr(self, att)))
return "SolarEph class object attributes"