Source code for ssapy.compute

import numpy as np
from astropy.time import Time
import astropy.units as u

from .body import get_body
from .constants import RGEO, LD, EARTH_RADIUS, EARTH_MU, MOON_RADIUS, MOON_MU
from .propagator import KeplerianPropagator
from .utils import (
    norm, normed, unitAngle3, get_angle, LRU_Cache, lb_to_unit, sunPos, _gpsToTT,
    iers_interp, rotation_matrix_from_vectors, angle_between_vectors, gcrf_to_itrf
)
from .orbit import Orbit
from .ellipsoid import Ellipsoid

try:
    import erfa
except ImportError:
    # Let this raise
    import astropy._erfa as erfa


def _doSqueeze(squeezeOrbit, squeezeTime, *args):
    out = []
    for i, arg in enumerate(args):
        if squeezeTime:
            if squeezeOrbit:
                out.append(arg[0, 0])  # scalar
            else:
                out.append(np.squeeze(arg, axis=1))  # (nOrbit,)
        else:
            if squeezeOrbit:
                out.append(np.squeeze(arg, axis=0))  # (nTime,)
            else:
                out.append(arg)  # (nOrbit, nTime)
    if len(out) == 1:
        return out[0]
    return tuple(out)


def _countOrbit(orbit):
    # orbit is one of:
    # 1) scalar Orbit
    # 2) vector Orbit
    # 3) list of scalar Orbit
    # convert to (2), set nOrbit, squeezeOrbit, and orbit.

    # check 1) and 2)
    if isinstance(orbit, Orbit):
        # check shape of r to decide scalar vs vector
        if orbit.r.ndim == 1:  # scalar Orbit
            nOrbit = 1
            squeezeOrbit = True
            if 'kozaiMeanKeplerianElements' in orbit.__dict__:
                kMKE = orbit.__dict__['kozaiMeanKeplerianElements']
            else:
                kMKE = None
            orbit = Orbit(
                np.atleast_2d(orbit.r),
                np.atleast_2d(orbit.v),
                np.atleast_1d(orbit.t),
                mu=orbit.mu,
                propkw={k: np.atleast_1d(v) for k, v in orbit.propkw.items()},
            )
            # Copying just kozaiMeanKeplerianElements for now, though maybe
            # ought to retain other (all?) lazy_properties attributes too?
            if kMKE is not None:
                orbit.kozaiMeanKeplerianElements = kMKE
        else:  # vector Orbit
            nOrbit = orbit.r.shape[0]
            squeezeOrbit = False
    else:
        import warnings
        warnings.warn(
            "list of Orbit syntax is deprecated.  "
            "Please use a vector Orbit instead.",
            DeprecationWarning
        )
        # list of scalar orbit
        nOrbit = len(orbit)
        squeezeOrbit = False
        # assumes all orbits have the same keywords and mu!
        propkeys = [] if nOrbit == 0 else orbit[0].propkw.keys()
        mu = None if nOrbit == 0 else orbit[0].mu
        orbit = Orbit(
            np.array([o.r for o in orbit]),
            np.array([o.v for o in orbit]),
            np.array([o.t for o in orbit]),
            mu=mu,
            propkw={k: np.array([o.propkw[k] for o in orbit])
                    for k in propkeys}
        )
    return nOrbit, squeezeOrbit, orbit


def _countTime(time):
    if isinstance(time, Time):
        time = time.gps
    squeezeTime = False
    try:
        nTime = len(time)
    except TypeError:
        time = np.array([time])
        nTime = 1
        squeezeTime = True
    return nTime, squeezeTime, time


def _countR(r):
    # orbit is one of:
    # 1) scalar r
    # 2) vector r
    # 3) list of scalar Orbit
    # convert to (2), set nOrbit, squeezeOrbit, and orbit.
    squeezeR = False
    if np.shape(r)[-1] == 3:
        pass
    else:
        raise ValueError(f"Incorrect r dimensions. Expected shape (n, 3), but got {np.shape(r)}.")
    # check 1) and 2)
    if r.ndim < 3:  # scalar r
        nR = 1
        r = np.reshape(np.atleast_3d(r), (nR, np.shape(r)[0], np.shape(r)[1]))
        squeezeR = True
    else:
        nR = np.shape(r)[0]
    return nR, squeezeR, r


def _processObserver(
    observer, obsPos, obsVel, nTime, squeezeTime, time, doObsVel=False
):
    # We always need obsPos output, but only sometimes obsVel.  So only require
    # obsVel if doObsVel is True.  (but always return it, and always compute it
    # if what we got as input was in an observer).
    squeezeObsPos = False
    if observer is None:
        if obsPos is None:
            raise ValueError(
                "Exactly one of obsPos and observer must be specified"
            )
        if doObsVel and obsVel is None:
            raise ValueError(
                "Exactly one of obsVel and observer required for doObsVel"
            )

    if observer is not None:
        if obsPos is not None:
            raise ValueError(
                "Exactly one of obsPos and observer must be specified"
            )
        if doObsVel and obsVel is not None:
            raise ValueError(
                "Exactly one of obsVel and observer required for doObsVel"
            )
        try:
            nObservers = len(observer)
        except TypeError:
            observer = [observer]
            nObservers = 1
            squeezeObsPos = squeezeTime
        # Now generate obsPos
        if nObservers == nTime:
            obsPos = np.empty((nObservers, 3), dtype=float)
            obsVel = np.empty((nObservers, 3), dtype=float)
            for i, (obs, t) in enumerate(zip(observer, time)):
                obsPos[i], obsVel[i] = obs.getRV(t)
        else:
            if nObservers == 1:
                obsPos, obsVel = observer[0].getRV(time)
            elif nTime == 1:
                obsPos = np.empty((nObservers, 3), dtype=float)
                obsVel = np.empty((nObservers, 3), dtype=float)
                for i, obs in enumerate(observer):
                    obsPos[i], obsVel[i] = obs.getRV(time[0])
            else:
                raise ValueError("observer and time must be broadcastable")

    if doObsVel:
        obsPos, obsVel = np.broadcast_arrays(obsPos, obsVel)
    obsPos = np.atleast_1d(obsPos)
    if obsPos.ndim == 2:
        nObsPos = obsPos.shape[0]
    else:
        obsPos = np.atleast_2d(obsPos)
        nObsPos = 1
        squeezeObsPos = True
    # times and obsPos must align
    if nObsPos != nTime:
        if nObsPos == 1:
            obsPos = np.broadcast_to(obsPos, (nTime, 3))
            if doObsVel:
                obsVel = np.broadcast_to(obsVel, (nTime, 3))
            nObsPos = nTime
        elif nTime == 1:
            time = np.broadcast_to(time, (nObsPos,))
            nTime = nObsPos
        else:
            raise ValueError("obsPos and time must be broadcastable")
    squeezeTime &= squeezeObsPos
    return obsPos, obsVel, time, squeezeTime


class HashableArrayContainer:
    def __init__(self, arr):
        self.arr = arr
        self.arr.flags.writeable = False

    def __hash__(self):
        return hash(self.arr.data.tobytes())

    def __eq__(self, rhs):
        return np.all(self.arr == rhs.arr)


[docs] def rv(orbit, time, propagator=KeplerianPropagator()): """Calculate positions and velocities on the outer product of all supplied orbits and times. Parameters ---------- orbit : Orbit or list of Orbit (n,) Desired orbit(s) time : array_like or astropy.time.Time (m,) If float (array), then should correspond to GPS seconds; i.e., seconds since 1980-01-06 00:00:00 UTC propagator : Propagator, optional The propagator instance to use. Notes ----- If orbit or time is scalar valued (as opposed to a list of Orbit, e.g.) then the corresponding dimension of the output will be squeezed out. For Keplerian orbit propagation it is more efficient to use a "vector Orbit" instead of a list of single scalar Orbits. Returns ------- r : array_like (n, m, 3) Position in meters. v : array_like (n, m, 3) Velocity in meters per second. """ nOrbit, squeezeOrbit, orbit = _countOrbit(orbit) nTime, squeezeTime, time = _countTime(time) outR, outV = _rv(orbit, HashableArrayContainer(time), propagator) return _doSqueeze(squeezeOrbit, squeezeTime, outR, outV)
def __rv(orbit, time, propagator): outR, outV = propagator._getRVMany(orbit, time.arr) return outR, outV _rv = LRU_Cache(__rv, maxsize=16)
[docs] def groundTrack(orbit, time, propagator=KeplerianPropagator(), format='geodetic'): """Calculate satellite ground track on the outer product of all supplied times and state vectors or orbits. Parameters ---------- r : array_like (n,3) Position (m) time : array_like or astropy.time.Time (m,) If float (array), then should correspond to GPS seconds; i.e., seconds since 1980-01-06 00:00:00 UTC or orbit : r : array_like (n,3) Position (m) or Orbit or list of Orbit (n,) Desired orbit(s) propagator : Propagator, optional The propagator instance to use. format : 'geodetic' or 'cartesian' If 'geodetic', then returns longitude, latitude, height. If 'cartesian', then returns xyz in ITRF frame. Notes ----- If orbit or time is scalar valued (as opposed to a list of Orbit, e.g.) then the corresponding dimension of the output will be squeezed out. For Keplerian orbit propagation it is more efficient to use a "vector Orbit" instead of a list of single scalar Orbits. Returns ------- lon, lat, height : array_like (n, m, 3) Radians and meters. or x, y, z : array_like(n, m, 3) Meters. """ if format not in ['cartesian', 'geodetic']: raise ValueError("Format must be either 'cartesian' or 'geodetic'") nTime, squeezeTime, time = _countTime(time) if isinstance(orbit, Orbit): nOrbit, squeezeOrbit, orbit = _countOrbit(orbit) r, v = rv(orbit, time, propagator=propagator) # (n, m, 3) else: nOrbit, squeezeOrbit, r = _countR(orbit) # (n, m, 3) # Reverse the math in EarthObserver.getRV mjd_tt = _gpsToTT(time) d_ut1_tt_mjd, pmx, pmy = iers_interp(time) pn = erfa.pnm80(2400000.5, mjd_tt) gst = erfa.gst94(2400000.5, mjd_tt + d_ut1_tt_mjd) cg, sg = np.cos(gst), np.sin(gst) gstMat = np.zeros([len(time), 3, 3], dtype=float) gstMat[:, 0, 0] = cg gstMat[:, 0, 1] = sg gstMat[:, 1, 0] = -sg gstMat[:, 1, 1] = cg gstMat[:, 2, 2] = 1.0 polar = np.zeros([len(time), 3, 3], dtype=float) polar[:] = np.eye(3) polar[:, 0, 2] = pmx polar[:, 1, 2] = -pmy polar[:, 2, 0] = -pmx polar[:, 2, 1] = pmy U = gstMat @ pn itrs = np.einsum( "tab,ntb->nta", (polar @ U), r ) x = itrs[..., 0] y = itrs[..., 1] z = itrs[..., 2] if format == 'cartesian': return _doSqueeze(squeezeOrbit, squeezeTime, x, y, z) elif format == 'geodetic': ellipsoid = Ellipsoid() lon, lat, height = ellipsoid.cartToSphere(x, y, z) return _doSqueeze(squeezeOrbit, squeezeTime, lon, lat, height)
def _obsAngleCorrection( r, v, obsPos, obsVel, orbit, time, propagator, correctionType, max_iter=10 ): if correctionType == "linear": # First order correction is to back up satellite by v * dt # where dt is the approximate time delay from the 0th order distance. # Important: *don't* do anything to the observer position, since we # really did do the observation at the originally given time and place. dt = norm(r - obsPos) / 299792458 r -= np.einsum("nmi,nm->nmi", v, dt) # no change to v since linear assumption _is_ that v is constant for the # light time duration. elif correctionType == "exact": # Solve |r(t-dt) - r_obs(t)| = c dt for dt. dist = norm(r - obsPos) dt, dt_previous = dist / 299792458, np.inf iter = 0 while np.any(np.abs(dt - dt_previous)) > 1e-12: # picosecond accurate if iter > max_iter: raise RuntimeError( "Exact light time correction did not converge in " "{} iterations".format(iter) ) # Unfortunately, can't just do # r, v = rv(orbit, time-dt, propagator=propagator) # because ssapy.rv isn't parallelized for n orbits and n, m times. # so do one orbit at a time. r = np.empty((len(orbit), len(time), 3)) v = np.empty((len(orbit), len(time), 3)) for i, o in enumerate(orbit): r[i], v[i] = rv(o, time - dt[i], propagator=propagator) dist = norm(r - obsPos) dt, dt_previous = dist / 299792458, dt iter += 1 else: raise ValueError( "Invalid value for correctionType: {}" .format(correctionType) ) # Apply velocity correction (i.e., diurnal aberration for an Earth-based # observer, and whatever you call the analogous correction for an orbital # observer). dr = norm(r - obsPos) r += obsVel * dr[..., None] / 299792458 return r, v
[docs] def dircos( orbit, time, obsPos=None, obsVel=None, observer=None, propagator=KeplerianPropagator(), obsAngleCorrection=None ): """ Calculate observed direction-cosines of orbiting objects as viewed at specified times and positions. The direction cosines are the cosines of the angles between the vector pointing towards the orbiting object and the x, y, z axes. An equivalent description is that they are the components of the unit vector pointing towards the orbiting object. Parameters ---------- orbit : Orbit or list of Orbit (n,) Orbit(s) for which to calculate direction cosines. time : array_like or astropy.time.Time (m,) If float (array), then should correspond to GPS seconds; i.e., seconds since 1980-01-06 00:00:00 UTC obsPos : array_like (m, 3), optional Position of observer at given time. obsVel : array_like (m, 3), optional Velocity of observer at given time. Only required if correcting for diurnal aberration. observer : Observer or list of Observers (m,), optional Observer(s) used to calculate obsPos and obsVel. propagator : Propagator, optional The propagator instance to use. obsAngleCorrection : {None, "linear", "exact"}, optional Correct actual angle to observed angle, meaning account for light-time delay, and aberration due to the observer's velocity. None means don't do any correction. "linear" means do an aberration correction and first order light-time correction. "exact" means do an aberration correction and iteratively solve for the exact light-time correction. (The "linear" correction is almost always sufficiently accurate). Notes ----- Exactly 1 of `obsPos` and `observer` must be supplied. `observer` and `obsPos` follow similar broadcasting rules as detailed below explicitly only for `obsPos`. The length of `time` and `obsPos` must match or be broadcastable to match. If `orbit` is scalar-valued (an Orbit instead of a list of Orbit), then that dimension will be squeezed out in the return value. Likewise, if both `time` and `obsPos` are scalar, that dimension will be squeezed out. For Keplerian orbit propagation it is more efficient to use a "vector Orbit" instead of a list of single scalar Orbits. When doing light time corrections, the time argument is the arrival time of the photons at the observer, as opposed to the emission time at the satellite. Returns ------- dircos : array_like (n, m, 3) Direction cosines on the outer product of orbit(s) and time/obsPos. """ nTime, squeezeTime, time = _countTime(time) nOrbit, squeezeOrbit, orbit = _countOrbit(orbit) obsPos, obsVel, time, squeezeTime = _processObserver( observer, obsPos, obsVel, nTime, squeezeTime, time, doObsVel=obsAngleCorrection ) # At this point, should have len(orbit) == n # and len(obsPos) == len(time) == nTime # So get positions of orbits... r, v = rv(orbit, time, propagator=propagator) # (n, m, 3) if obsAngleCorrection: r, v = _obsAngleCorrection( r, v, obsPos, obsVel, orbit, time, propagator, obsAngleCorrection ) dc = normed(r - obsPos) return _doSqueeze(squeezeOrbit, squeezeTime, dc)
[docs] def radec( orbit, time, obsPos=None, obsVel=None, observer=None, propagator=KeplerianPropagator(), obsAngleCorrection=None, rate=False ): """ Calculate observed right ascension, declination, and slant range of orbiting objects as viewed at specified times and positions. Parameters ---------- orbit : Orbit or list of Orbit (n,) Orbit(s) for which to calculate ra and dec. time : array_like or astropy.time.Time (m,) If float (array), then should correspond to GPS seconds; i.e., seconds since 1980-01-06 00:00:00 UTC obsPos : array_like (m, 3), optional Positions of observers at given times. obsVel : array_like (m, 3), optional Velocity of observers at given times. observer : Observer or list of Observers (m,), optional Observer(s) used to calculate obsPos. propagator : Propagator, optional The propagator instance to use. rate : bool If True, return additionally the time derivatives of the quantity, times cos(dec) in the case of right ascension. obsAngleCorrection : {None, "linear", "exact"}, optional Correct actual angle to observed angle, meaning account for light-time delay, and aberration due to the observer's velocity. None means don't do any correction. "linear" means do an aberration correction and first order light-time correction. "exact" means do an aberration correction and iteratively solve for the exact light-time correction. (The "linear" correction is almost always sufficiently accurate). Notes ----- Exactly 1 of `obsPos` and `observer` must be supplied. `observer` and `obsPos` follow similar broadcasting rules as detailed below explicitly only for `obsPos`. The length of `time` and `obsPos` must match or be broadcastable to match. If `orbit` is scalar-valued (an Orbit instead of a list of Orbit), then that dimension will be squeezed out in the return value. Likewise, if both `time` and `obsPos` are scalar, that dimension will be squeezed out. For Keplerian orbit propagation it is more efficient to use a "vector Orbit" instead of a list of single scalar Orbits. When doing light time corrections, the time argument is the arrival time of the photons at the observer, as opposed to the emission time at the satellite. Returns ------- ra, dec : array_like (n, m) Right ascension and declination in radians. range : array_like (n, m) (Slant) range in meters. If rate, also: raRate : array_like (n, m) Time derivatives of right ascension times cos(dec), rad / sec. decRate : array_like (n, m) Time derivative of declination, rad / sec. rangeRate : array_like (n, m) Time derivative of range, meter / s. """ nTime, squeezeTime, time = _countTime(time) nOrbit, squeezeOrbit, orbit = _countOrbit(orbit) obsPos, obsVel, time, squeezeTime = _processObserver( observer, obsPos, obsVel, nTime, squeezeTime, time, doObsVel=obsAngleCorrection ) if rate and obsVel is None: raise ValueError('obsVel must not be None if rate is True.') # At this point, should have len(orbit) == n # and len(obsPos) == len(time) == nTime # So get positions of orbits... r, v = rv(orbit, time, propagator=propagator) # (n, m, 3) if obsAngleCorrection: r, v = _obsAngleCorrection( r, v, obsPos, obsVel, orbit, time, propagator, obsAngleCorrection ) ra, dec, slantRange, raRate, decRate, rangeRate = rvObsToRaDecRate( r, v, obsPos, obsVel ) res = (ra, dec, slantRange) if rate: res = res + (raRate, decRate, rangeRate) return _doSqueeze(*((squeezeOrbit, squeezeTime) + res))
[docs] def altaz( orbit, time, observer, propagator=KeplerianPropagator(), obsAngleCorrection=None ): """ Calculate observed altitude and azimuth of orbiting objects as viewed at specified times and locations. Parameters ---------- orbit : Orbit or list of Orbit (n,) Orbits for which to calculate direction cosines. time : array_like or astropy.time.Time (m,) If float (array), then should correspond to GPS seconds; i.e., seconds since 1980-01-06 00:00:00 UTC observer : EarthObserver Where on Earth to compute alt and az. propagator : Propagator, optional The propagator instance to use. obsAngleCorrection : {None, "linear", "exact"}, optional Correct actual angle to observed angle, meaning account for light-time delay, and aberration due to the observer's velocity. None means don't do any correction. "linear" means do an aberration correction and first order light-time correction. "exact" means do an aberration correction and iteratively solve for the exact light-time correction. (The "linear" correction is almost always sufficiently accurate). Notes ----- If `orbit` is scalar-valued (an Orbit instead of a list of Orbit), then that dimension will be squeezed out in the return value. Likewise, if `time` is scalar, that dimension will be squeezed out. For Keplerian orbit propagation it is more efficient to use a "vector Orbit" instead of a list of single scalar Orbits. When doing light time corrections, the time argument is the arrival time of the photons at the observer, as opposed to the emission time at the satellite. Returns ------- alt, az : array_like (n, m) Altitude and azimuth in radians. """ from astropy.coordinates import SkyCoord, AltAz, GCRS nTime, squeezeTime, time = _countTime(time) nOrbit, squeezeOrbit, orbit = _countOrbit(orbit) obsPos, obsVel = observer.getRV(time) r, v = rv(orbit, time, propagator=propagator) # (n, m, 3) if obsAngleCorrection: r, v = _obsAngleCorrection( r, v, obsPos, obsVel, orbit, time, propagator, obsAngleCorrection ) # Using astropy backend for AltAz computation sc = SkyCoord( x=r[..., 0], y=r[..., 1], z=r[..., 2], unit='m', representation_type='cartesian', frame=GCRS(obstime=Time(time, format='gps')) ) aa = sc.transform_to(AltAz(location=observer._location)) return _doSqueeze(squeezeOrbit, squeezeTime, aa.alt.radian, aa.az.radian)
[docs] def quickAltAz( orbit, time, observer, propagator=KeplerianPropagator(), obsAngleCorrection=None ): """ Quickly estimate observed altitude and azimuth of orbiting objects as viewed at specified times and locations. This algorithm approximates "up" as pointing directly away from the center of the Earth, instead of normal to the reference ellipsoid. Use `altAz` if you want values wrt the reference ellipsoid. Parameters ---------- orbit : Orbit or list of Orbit (n,) Orbits for which to calculate direction cosines. time : array_like or astropy.time.Time (m,) If float (array), then should correspond to GPS seconds; i.e., seconds since 1980-01-06 00:00:00 UTC observer : EarthObserver Where on Earth to compute alt and az. propagator : Propagator, optional The propagator instance to use. obsAngleCorrection : {None, "linear", "exact"}, optional Correct actual angle to observed angle, meaning account for light-time delay, and aberration due to the observer's velocity. None means don't do any correction. "linear" means do an aberration correction and first order light-time correction. "exact" means do an aberration correction and iteratively solve for the exact light-time correction. (The "linear" correction is almost always sufficiently accurate). Notes ----- If `orbit` is scalar-valued (an Orbit instead of a list of Orbit), then that dimension will be squeezed out in the return value. Likewise, if `time` is scalar, that dimension will be squeezed out. For Keplerian orbit propagation it is more efficient to use a "vector Orbit" instead of a list of single scalar Orbits. When doing light time corrections, the time argument is the arrival time of the photons at the observer, as opposed to the emission time at the satellite. Returns ------- alt, az : array_like (n, m) Altitude and azimuth in radians. """ nTime, squeezeTime, time = _countTime(time) nOrbit, squeezeOrbit, orbit = _countOrbit(orbit) ro, vo = observer.getRV(time) # (m, 3) r, v = rv(orbit, time, propagator=propagator) # (n, m, 3) if obsAngleCorrection: r, v = _obsAngleCorrection( r, v, ro, vo, orbit, time, propagator, obsAngleCorrection ) up = normed(ro) east = normed(vo) north = np.cross(up, east) dr = r - ro northing = np.einsum("nmi,mi->nm", dr, north) easting = np.einsum("nmi,mi->nm", dr, east) alt = np.pi / 2 - unitAngle3( normed(np.broadcast_to(ro, r.shape)), normed(dr) ) az = np.arctan2(easting, northing) az %= 2 * np.pi return _doSqueeze(squeezeOrbit, squeezeTime, alt, az)
[docs] def radecRate( orbit, time, obsPos=None, obsVel=None, observer=None, propagator=KeplerianPropagator(), obsAngleCorrection=None ): """Calculate ra/dec rate and slant range rate of orbit at specified times and observer positions and velocities. DEPRECATED. Use radec(..., rate=True) in new code. Parameters ---------- orbit : Orbit or list of Orbit (n,) Orbit(s) for which to calculate slant range rate. time : array_like or astropy.time.Time (m,) If float (array), then should correspond to GPS seconds; i.e., seconds since 1980-01-06 00:00:00 UTC obsPos : array_like (m, 3), optional Positions of observers at given times. obsVel : array_like (m, 3), optional Velocity of observers at given times. observer : Observer or list of Observers (m,), optional Observer(s) used to calculate obsPos and obsVel. propagator : Propagator, optional The propagator instance to use. obsAngleCorrection : {None, "linear", "exact"}, optional Correct actual angle to observed angle, meaning account for light-time delay, and aberration due to the observer's velocity. None means don't do any correction. "linear" means do an aberration correction and first order light-time correction. "exact" means do an aberration correction and iteratively solve for the exact light-time correction. (The "linear" correction is almost always sufficiently accurate). Notes ----- Exactly 1 of `obsPos` and `observer` must be supplied. `observer` and `obsPos` follow similar broadcasting rules as detailed below explicitly only for `obsPos`. If `obsPos` is specified, `obsVel` must also be specified and congruent to `obsPos`. The length of `time` and `obsPos` must match or be broadcastable to match. If `orbit` is scalar-valued (an Orbit instead of a list of Orbit), then that dimension will be squeezed out in the return value. Likewise, if both `time` and `obsPos` are scalar, that dimension will be squeezed out. For Keplerian orbit propagation it is more efficient to use a "vector Orbit" instead of a list of single scalar Orbits. When doing light time corrections, the time argument is the arrival time of the photons at the observer, as opposed to the emission time at the satellite. Returns ------- ra : array_like (n, m) right ascension in radians raRate : array_like (n, m) Rate of change of right ascension*cos(dec) in radians per second. dec : array_link (n, m) declination in radians decRate : array_like (n, m) Rate of change of declination in radians per second. slantRange : array_like (n, m) Range in meters slantRangeRate : array_like (n, m) Slant range rate in meters per second. """ import warnings warnings.warn("This function is deprecated; use " "ssapy.compute.radec(..., rate=True)") ra, dec, slant, raRate, decRate, slantRate = radec( orbit, time, obsPos=obsPos, obsVel=obsVel, observer=observer, propagator=propagator, obsAngleCorrection=obsAngleCorrection ) return raRate, decRate, slantRate
[docs] def rvObsToRaDecRate(r, v, obsPos=None, obsVel=None): """Convert object and observer position and velocity to angles. This only does the geometric part; it ignores light travel time, which may be applied to the object location before input. Assumes that r, v, obsPos, and obsVel all have common shape. Parameters ---------- r : array_like (..., 3) object position in meters v : array_like (..., 3) object velocity in meters per second obsPos : array_like (..., 3) observer position in meters obsVel : array_like (..., 3), optional observer velocity in meters per second Returns ------- ra : array_like (...) right ascension in radians dec : array_link (...) declination in radians slantRange : array_like (...) Range in meters raRate : array_like (...) Rate of change of right ascension*cos(dec) in radians per second. decRate : array_like (...) Rate of change of declination in radians per second. slantRangeRate : array_like (...) Slant range rate in meters per second. """ if obsPos is None: obsPos = np.zeros_like(r) dr = r - obsPos slantRange = norm(dr) ra = np.arctan2(dr[..., 1], dr[..., 0]) dec = np.arcsin(dr[..., 2] / slantRange) if obsVel is None: obsVel = np.zeros_like(v) dv = v - obsVel # Now need to rotate to ra/dec coords sd, cd = np.sin(dec), np.cos(dec) sa, ca = np.sin(ra), np.cos(ra) rHat = normed(dr) # (n, m, 3) raHat = np.zeros_like(rHat, dtype=float) raHat[..., 0] = -sa raHat[..., 1] = ca decHat = np.zeros_like(raHat) decHat[..., 0] = -sd * ca decHat[..., 1] = -sd * sa decHat[..., 2] = cd raRate = np.einsum("...i,...i->...", dv, raHat) / slantRange decRate = np.einsum("...i,...i->...", dv, decHat) / slantRange rangeRate = np.einsum("...i,...i->...", dv, rHat) return ra, dec, slantRange, raRate, decRate, rangeRate
[docs] def radecRateObsToRV(ra, dec, slantRange, raRate=None, decRate=None, slantRangeRate=None, obsPos=None, obsVel=None): """Convert object angles and observer position to 3D observer position This only does the geometric part; it ignores light travel time. This is the inverse of rvObsToRaDecRate. If obsVel is None, then the returned velocity will also be None. Parameters ---------- ra : array_like (...) right ascension in radians dec : array_like (...) declination in radians slantRange : array_like (...) Range in meters raRate : array_like (...) Rate of change of right ascension*cos(dec) in radians per second. decRate : array_like (...) Rate of change of declination in radians per second. slantRangeRate : array_like (...) Slant range rate in meters per second. obsPos : array_like (..., 3) Observer position in meters obsVel : array_like (..., 3) Observer velocity in meters Returns ------- r : array_like (..., 3) object position in meters v : array_like (..., 3) object velocity in meters per second observer velocity in meters per second v is None if obsVel is None. """ if obsPos is None: return ValueError('obsPos must be set!') rHat = lb_to_unit(ra, dec) sd, cd = np.sin(dec), np.cos(dec) sa, ca = np.sin(ra), np.cos(ra) raHat = np.zeros_like(rHat, dtype=float) raHat[..., 0] = -sa raHat[..., 1] = ca decHat = np.zeros_like(raHat) decHat[..., 0] = -sd * ca decHat[..., 1] = -sd * sa decHat[..., 2] = cd if len(rHat.shape) > 1: slantRange = np.array(slantRange)[..., None] raRate = np.array(raRate)[..., None] decRate = np.array(decRate)[..., None] slantRangeRate = np.array(slantRangeRate)[..., None] r = obsPos + rHat * slantRange if obsVel is None: v = None else: v = (obsVel + rHat * slantRangeRate + slantRange * (raHat * raRate + decHat * decRate)) return r, v
[docs] def earthShadowCoords(r, time): """Determine components of position `r` parallel and perpendicular to sun unit vector. The sun unit vector points from the center of the sun through the center of the Earth. Decomposing a satellite position into these coordinates yields a simple model of whether or not the satellite is in Earth's shadow: r_par, r_perp = earthShadowCoords(r, time) inShadow = r_par > 0 and r_perp < EARTH_RADIUS Parameters ---------- r : ndarray (3,) Position (GCRF) in meters time : float or astropy.time.Time If float, then should correspond to GPS seconds; i.e., seconds since 1980-01-06 00:00:00 UTC Returns ------- r_par : float Position of satellite projected onto the sun unit vector in meters. r_perp : float Distance of satellite from Earth-Sun line in meters. """ if isinstance(time, Time): time = time.gps r_sun = sunPos(time, fast=False) if r.ndim == 1: n_sun = - r_sun / norm(r_sun) r_par = np.dot(r, n_sun) r_perp = norm(r - n_sun * r_par) else: n_sun = - r_sun / norm(r_sun)[:, None] r_par = np.sum(r * n_sun, axis=1) r_perp = norm(r - n_sun * r_par[:, None]) return r_par, r_perp
[docs] def find_passes( orbit, observers, tStart, tSpan, dt, propagator=KeplerianPropagator(), horizon=np.deg2rad(20) ): """Find satellite overhead passes for a collection of observers. Uses a brute force test of a grid of time points from tStart to tStart+tSpan separated by dt. Returns passes even if they occur during the daytime or if the satellite is not illuminated by the sun. The only criterion for a successful "pass" is for the topocentric altitude of the satellite to be above the input `horizon`. More details about a pass can subsequently be obtained by running the `refine_passes` function. Note this function is only suitable for `EarthObserver`s and not `OrbitalObserver`s. Parameters ---------- orbit : Orbit Satellite orbit in question. observers : List of EarthObserver. Earth observers for which to check satellite visibility. tStart : float or astropy.time.Time Beginning of search window. If float, then should correspond to GPS seconds; i.e., seconds since 1980-01-06 00:00:00 UTC tSpan : float or Quantity Time span in which to search for passes. If float, then seconds. dt : float or Quantity Time increment to use during search. Satellite visibility will be computed every dt increment. Smaller values will decrease the probability that a short duration pass will be missed, but will make the search take longer to complete. If float, then seconds. propagator : Propagator, optional The propagator instance to use. horizon : float or Quantity, optional Minimum altitude for which to consider a satellite "visible". If float, then should be in radians. Returns ------- passDict : dict keys are `EarthObserver`s. values are lists (possibly empty) of `astropy.Time` corresponding to visible passes of the satellite. Only one time is returned per pass, so multiple times in the return list indicate multiple distinct passes. """ import astropy.units as u if isinstance(tStart, Time): tStart = tStart.gps if isinstance(tSpan, u.Quantity): tSpan = tSpan.to(u.s).value if isinstance(dt, u.Quantity): dt = dt.to(u.s).value if isinstance(horizon, u.Quantity): horizon = horizon.to(u.rad).value times = np.arange(tStart, tStart + tSpan, dt) out = {} for observer in observers: out[observer] = [] alt, _ = quickAltAz(orbit, times, observer, propagator=propagator) w = np.nonzero(alt > horizon)[0] for wi in w: if (wi - 1) in w: # already caught this pass continue out[observer].append(Time(times[wi], format='gps')) return out
[docs] def refine_pass( orbit, observer, time, propagator=KeplerianPropagator(), horizon=np.deg2rad(20), maxSpan=86400.0 ): """Refine a satellite overhead pass. Parameters ---------- orbit : Orbit Orbit in question. observer : EarthObserver. Observer for which to refine pass. time : float or astropy.time.Time A time when satellite is visible to observer. (Found using find_passes, for instance) If float, then should correspond to GPS seconds; i.e., seconds since 1980-01-06 00:00:00 UTC propagator : Propagator, optional The propagator instance to use. horizon : float or Quantity, optional Minimum altitude for which to consider a satellite "visible". If float, then should be in radians. maxSpan : float or Quantity, optional Maximum amount of time before or after `time` to search for rise/set times, or time of max altitude. If float, then seconds. Returns ------- dict Key/values are: tStart : Time Time at which Orbit rises above horizon tEnd : Time Time at which Orbit rises above horizon tMaxAlt : Time Time Orbit passes through maximum altitude maxAlt : Quantity Maximum altitude duration : Quantity Duration of pass. illumAtStart : bool Is satellite illuminated at `tStart`? illumAtEnd : bool Is satellite illuminated at `tEnd`? tTerminator : Time or None If illumAtStart != illumAtEnd, then time satellite passes through cylindrical terminator shadow. Otherwise, None. sunAltStart : Quantity Altitude of sun at `tStart` sunAltEnd : Quantity Altitude of sun at `tEnd` """ from scipy.optimize import bisect, minimize_scalar import astropy.units as u import warnings if isinstance(time, Time): time = time.gps if isinstance(horizon, u.Quantity): horizon = horizon.to(u.rad).value if isinstance(maxSpan, u.Quantity): maxSpan = maxSpan.to(u.s).value def dalt(t): return ( horizon - quickAltAz(orbit, t, observer, propagator=propagator)[0] ) # Height "outside" of cylindrical Earth shadow model. def dshadow(t): r, _ = rv(orbit, t, propagator=propagator) r_par, r_perp = earthShadowCoords(r, t) height = r_perp - EARTH_RADIUS if r_par < 0: # a bit of a hack, if the sat is on the sun side, then we'll # return abs(dshadow) so that the value is positive, which indicates # illumination. return abs(height) else: return height # Bracket and then bisect to find rise/set times def bracket(f, x, dx, xmax): x0, x1 = x, x + dx f0, f1 = f(x0), f(x1) while f0 * f1 > 0 and abs(x1 - x0) < xmax: x1 += dx f1 = f(x1) return x1 tLow = bracket(dalt, time, -300, maxSpan) if dalt(tLow) * dalt(time) > 0: # didn't bracket, use tLow and issue warning warnings.warn("Failed to bracket. tStart is not rise time!") tStart = tLow else: tStart = bisect(dalt, tLow, time) tStart = Time(tStart, format='gps') tStart.format = 'iso' tHigh = bracket(dalt, time, 300, maxSpan) if dalt(tHigh) * dalt(time) > 0: warnings.warn("Failed to bracket. tEnd is not set time!") tEnd = tHigh else: tEnd = bisect(dalt, time, tHigh) tEnd = Time(tEnd, format='gps') tEnd.format = 'iso' # Find maximum altitude dalt1 = dalt(tStart.gps) dalt2 = dalt(time) dalt3 = dalt(tEnd.gps) if dalt2 < dalt1 and dalt2 < dalt3: result = minimize_scalar(dalt, (tStart.gps, time, tEnd.gps)) tMaxAlt = Time(result.x, format='gps') tMaxAlt.format = 'iso' maxAlt = -result.fun + horizon else: ama = np.argmin([dalt1, dalt2, dalt3]) tMaxAlt = Time([tStart.gps, time, tEnd.gps][ama], format='gps') tMaxAlt.format = 'iso' maxAlt = -[dalt1, dalt2, dalt3][ama] + horizon duration = (tEnd - tStart).to(u.min) # Find illumination state illumAtStart = dshadow(tStart.gps) > 0 sunAltStart = np.rad2deg(observer.sunAlt(tStart.gps)) * u.deg illumAtEnd = dshadow(tEnd.gps) > 0 sunAltEnd = np.rad2deg(observer.sunAlt(tEnd.gps)) * u.deg # If we transition between shadow/illuminated, then find the time of # transition (time of terminator pass) if illumAtStart != illumAtEnd: tTerminator = bisect(dshadow, tStart.gps, tEnd.gps) tTerminator = Time(tTerminator, format='gps') tTerminator.format = 'iso' else: tTerminator = None return { "orbit": orbit, "observer": observer, "propagator": propagator, "horizon": horizon, "tStart": tStart, "tEnd": tEnd, "altStart": np.rad2deg(-dalt(tStart.gps) + horizon) * u.deg, "altEnd": np.rad2deg(-dalt(tEnd.gps) + horizon) * u.deg, "tMaxAlt": tMaxAlt, "maxAlt": np.rad2deg(maxAlt) * u.deg, "duration": duration, "illumAtStart": illumAtStart, "illumAtEnd": illumAtEnd, "tTerminator": tTerminator, "sunAltStart": sunAltStart, "sunAltEnd": sunAltEnd }
def nby3shape(arr_): if arr_.ndim == 1: return np.reshape(arr_, (1, 3)) if arr_.ndim == 2: if np.shape(arr_)[1] == 3: return arr_ else: return arr_.T
[docs] def calculate_orbital_elements(r_, v_, mu_barycenter=EARTH_MU): """ Calculate the Keplerian orbital elements from position and velocity vectors. This function computes the Keplerian orbital elements (semi-major axis, eccentricity, inclination, true longitude, argument of periapsis, longitude of ascending node, true anomaly, and specific angular momentum) for one or more celestial objects given their position and velocity vectors. Parameters: ---------- r_ : (n, 3) numpy.ndarray Array of position vectors (in meters) of the celestial objects. Each row represents a position vector. v_ : (n, 3) numpy.ndarray Array of velocity vectors (in meters per second) of the celestial objects. Each row represents a velocity vector. mu_barycenter : float, optional Gravitational parameter (standard gravitational constant times mass) of the central body (default is `EARTH_MU`). This parameter defines the gravitational influence of the central body on the orbiting object. Returns: ------- dict A dictionary containing the orbital elements for each celestial object: - 'a': Semi-major axis (in meters). - 'e': Eccentricity (dimensionless). - 'i': Inclination (in radians). - 'tl': True longitude (in radians). - 'ap': Argument of periapsis (in radians). - 'raan': Longitude of ascending node (in radians). - 'ta': True anomaly (in radians). - 'L': Specific angular momentum (in meters squared per second). Notes: ------ - Position and velocity vectors should be provided in the same units. - The function assumes that the input vectors are provided in an array where each row corresponds to a different celestial object. - Orbital elements are computed using standard orbital mechanics formulas. - The inclination is measured from the reference plane, and the argument of periapsis and true anomaly are measured in the orbital plane. Example: -------- >>> r = np.array([[1e7, 1e7, 1e7], [1e8, 1e8, 1e8]]) >>> v = np.array([[1e3, 2e3, 3e3], [4e3, 5e3, 6e3]]) >>> calculate_orbital_elements(r, v, mu_barycenter=3.986e14) {'a': [1.5707e7, 2.234e8], 'e': [0.123, 0.456], 'i': [0.785, 0.654], 'tl': [2.345, 3.456], 'ap': [0.123, 0.456], 'raan': [1.234, 2.345], 'ta': [0.567, 1.678], 'L': [1.234e10, 2.345e11]} """ # mu_barycenter - all bodies interior to Earth # 1.0013415732186798 #All bodies of solar system mu_ = mu_barycenter rarr = nby3shape(r_) varr = nby3shape(v_) aarr = [] earr = [] incarr = [] true_longitudearr = [] argument_of_periapsisarr = [] longitude_of_ascending_nodearr = [] true_anomalyarr = [] hmagarr = [] for r, v in zip(rarr, varr): r = np.array(r) # print(f'r: {r}') v = np.array(v) # print(f'v: {v}') rmag = np.sqrt(r.dot(r)) vmag = np.sqrt(v.dot(v)) h = np.cross(r, v) hmag = np.sqrt(h.dot(h)) n = np.cross(np.array([0, 0, 1]), h) a = 1 / ((2 / rmag) - (vmag ** 2) / mu_) evector = np.cross(v, h) / (mu_) - r / rmag e = np.sqrt(evector.dot(evector)) inc = np.arccos(h[2] / hmag) if np.dot(r, v) > 0: true_anomaly = np.arccos(np.dot(evector, r) / (e * rmag)) else: true_anomaly = 2 * np.pi - np.arccos(np.dot(evector, r) / (e * rmag)) if evector[2] >= 0: argument_of_periapsis = np.arccos(np.dot(n, evector) / (e * np.sqrt(n.dot(n)))) else: argument_of_periapsis = 2 * np.pi - np.arccos(np.dot(n, evector) / (e * np.sqrt(n.dot(n)))) if n[1] >= 0: longitude_of_ascending_node = np.arccos(n[0] / np.sqrt(n.dot(n))) else: longitude_of_ascending_node = 2 * np.pi - np.arccos(n[0] / np.sqrt(n.dot(n))) true_longitude = true_anomaly + argument_of_periapsis + longitude_of_ascending_node aarr.append(a) earr.append(e) incarr.append(inc) true_longitudearr.append(true_longitude) argument_of_periapsisarr.append(argument_of_periapsis) longitude_of_ascending_nodearr.append(longitude_of_ascending_node) true_anomalyarr.append(true_anomaly) hmagarr.append(hmag) return { 'a': aarr, 'e': earr, 'i': incarr, 'tl': true_longitudearr, 'ap': argument_of_periapsisarr, 'raan': longitude_of_ascending_nodearr, 'ta': true_anomalyarr, 'L': hmagarr }
###################################################################################### # Lambertian brightness functions ######################################################################################
[docs] def moon_shine(r_moon, r_sat, r_earth, r_sun, radius, albedo, albedo_moon, albedo_back, albedo_front, area_panels): # In SI units, takes single values or arrays returns a fractional flux """ Calculate the fractional flux of sunlight reflected from the Moon to the satellite. This function computes the flux of sunlight reflected from the Moon to the satellite, including contributions from both the front and back surfaces of the satellite's solar panels. Parameters: ---------- r_moon : (n, 3) numpy.ndarray Array of coordinates representing the position of the Moon. r_sat : (n, 3) numpy.ndarray Array of coordinates representing the position of the satellite. r_earth : (n, 3) numpy.ndarray Array of coordinates representing the position of the Earth. r_sun : (n, 3) numpy.ndarray Array of coordinates representing the position of the Sun. radius : float Radius of the satellite in meters. albedo : float Albedo of the satellite's surface. albedo_moon : float Albedo of the Moon. albedo_back : float Albedo of the back surface of the satellite's solar panels. albedo_front : float Albedo of the front surface of the satellite's solar panels. area_panels : float Area of the satellite's solar panels in square meters. Returns: ------- dict Dictionary containing the flux contributions from the Moon to the satellite: - 'moon_bus': Fraction of light reflected off the satellite's bus from the Moon. - 'moon_panels': Fraction of light reflected off the satellite's panels from the Moon. Notes: ------ - The function assumes that the solar panels are always facing the Sun and calculates flux based on the phase angles. - Flux contributions from both the front and back surfaces of the solar panels are computed. Example: -------- >>> r_moon = np.array([[1e8, 1e8, 1e8]]) >>> r_sat = np.array([[1e7, 1e7, 1e7]]) >>> r_earth = np.array([[0, 0, 0]]) >>> r_sun = np.array([[1e11, 0, 0]]) >>> moon_shine(r_moon, r_sat, r_earth, r_sun, radius=0.4, albedo=0.20, albedo_moon=0.12, albedo_back=0.50, albedo_front=0.05, area_panels=100) {'moon_bus': array([...]), 'moon_panels': array([...])} """ # https://amostech.com/TechnicalPapers/2013/POSTER/COGNION.pdf moon_phase_angle = get_angle(r_sun, r_moon, r_sat) # Phase of the moon as viewed from the sat. sun_angle = get_angle(r_sun, r_sat, r_moon) # angle from Sun to object to Earth moon_to_earth_angle = get_angle(r_moon, r_sat, r_earth) r_moon_sat = np.linalg.norm(r_sat - r_moon, axis=-1) r_earth_sat = np.linalg.norm(r_sat - r_earth, axis=-1) # Earth is the observer. flux_moon_to_sat = 2 / 3 * albedo_moon * MOON_RADIUS**2 / (np.pi * (r_moon_sat)**2) * (np.sin(moon_phase_angle) + (np.pi - moon_phase_angle) * np.cos(moon_phase_angle)) # Fraction of sunlight reflected from the Moon to satellite # Fraction of light from back of solar panel flux_back = np.zeros_like(sun_angle) flux_back[sun_angle > np.pi / 2] = np.abs(albedo_back * area_panels / (np.pi * r_earth_sat[sun_angle > np.pi / 2]**2) * np.cos(np.pi - moon_to_earth_angle[sun_angle > np.pi / 2]) * flux_moon_to_sat[sun_angle > np.pi / 2]) # Fraction of Moon light reflected off back of solar panels - which are assumed to be always facing the Sun. Angle: Sun - Observer - Sat flux_front = np.zeros_like(sun_angle) flux_front[sun_angle < np.pi / 2] = np.abs(albedo_front * area_panels / (np.pi * r_earth_sat[sun_angle < np.pi / 2]**2) * np.cos(moon_to_earth_angle[sun_angle < np.pi / 2]) * flux_moon_to_sat[sun_angle < np.pi / 2]) # Fraction of Sun light scattered off front of the solar panels - which are assumed to be always facing the Sun. Angle: Sun - Sat - Observer flux_panels = flux_back + flux_front flux_bus = 2 / 3 * albedo * radius**2 / (np.pi * r_earth_sat**2) * flux_moon_to_sat return {'moon_bus': flux_bus, 'moon_panels': flux_panels}
[docs] def earth_shine(r_sat, r_earth, r_sun, radius, albedo, albedo_earth, albedo_back, area_panels): # In SI units, takes single values or arrays returns a flux """ Calculate the fractional flux of sunlight reflected from the Earth to the satellite. This function computes the flux of sunlight reflected from the Earth to the satellite, including contributions from the back surface of the satellite's solar panels. Parameters: ---------- r_sat : (n, 3) numpy.ndarray Array of coordinates representing the position of the satellite. r_earth : (n, 3) numpy.ndarray Array of coordinates representing the position of the Earth. r_sun : (n, 3) numpy.ndarray Array of coordinates representing the position of the Sun. radius : float Radius of the satellite in meters. albedo : float Albedo of the satellite's surface. albedo_earth : float Albedo of the Earth. albedo_back : float Albedo of the back surface of the satellite's solar panels. area_panels : float Area of the satellite's solar panels in square meters. Returns: ------- dict Dictionary containing the flux contributions from the Earth to the satellite: - 'earth_bus': Fraction of light reflected off the satellite's bus from the Earth. - 'earth_panels': Fraction of light reflected off the satellite's panels from the Earth. Notes: ------ - The function assumes that the solar panels are always facing the Sun and calculates flux based on the phase angle. - Flux contributions from the back surface of the solar panels are computed. Example: -------- >>> r_sat = np.array([[1e7, 1e7, 1e7]]) >>> r_earth = np.array([[1.496e11, 0, 0]]) >>> r_sun = np.array([[0, 0, 0]]) >>> earth_shine(r_sat, r_earth, r_sun, radius=0.4, albedo=0.20, albedo_earth=0.30, albedo_back=0.50, area_panels=100) {'earth_bus': array([...]), 'earth_panels': array([...])} """ # https://amostech.com/TechnicalPapers/2013/POSTER/COGNION.pdf phase_angle = get_angle(r_sun, r_sat, r_earth) # angle from Sun to object to Earth earth_angle = np.pi - phase_angle # Sun to Earth to object. r_earth_sat = np.linalg.norm(r_sat - r_earth, axis=-1) # Earth is the observer. flux_earth_to_sat = 2 / 3 * albedo_earth * EARTH_RADIUS**2 / (np.pi * (r_earth_sat)**2) * (np.sin(earth_angle) + (np.pi - earth_angle) * np.cos(earth_angle)) # Fraction of sunlight reflected from the Earth to satellite # Fraction of light from back of solar panel flux_back = np.zeros_like(phase_angle) flux_back[phase_angle > np.pi / 2] = albedo_back * area_panels / (np.pi * r_earth_sat[phase_angle > np.pi / 2]**2) * np.cos(np.pi - phase_angle[phase_angle > np.pi / 2]) * flux_earth_to_sat[phase_angle > np.pi / 2] # Fraction of Earth light reflected off back of solar panels - which are assumed to be always facing the Sun. Angle: Sun - Observer - Sat flux_bus = 2 / 3 * albedo * radius**2 / (np.pi * r_earth_sat**2) * flux_earth_to_sat return {'earth_bus': flux_bus, 'earth_panels': flux_back}
[docs] def sun_shine(r_sat, r_earth, r_sun, radius, albedo, albedo_front, area_panels): # In SI units, takes single values or arrays returns a fractional flux """ Calculate the fractional flux of sunlight reflected from the Sun to the satellite. This function computes the flux of sunlight reflected from the Sun to the satellite, including contributions from the front surface of the satellite's solar panels. Parameters: ---------- r_sat : (n, 3) numpy.ndarray Array of coordinates representing the position of the satellite. r_earth : (n, 3) numpy.ndarray Array of coordinates representing the position of the Earth. r_sun : (n, 3) numpy.ndarray Array of coordinates representing the position of the Sun. radius : float Radius of the satellite in meters. albedo : float Albedo of the satellite's surface. albedo_front : float Albedo of the front surface of the satellite's solar panels. area_panels : float Area of the satellite's solar panels in square meters. Returns: ------- dict Dictionary containing the flux contributions from the Sun to the satellite: - 'sun_bus': Fraction of light reflected off the satellite's bus from the Sun. - 'sun_panels': Fraction of light reflected off the satellite's panels from the Sun. Notes: ------ - The function assumes that the solar panels are always facing the Sun and calculates flux based on the phase angle. - Flux contributions from the front surface of the solar panels are computed. Example: -------- >>> r_sat = np.array([[1e7, 1e7, 1e7]]) >>> r_earth = np.array([[0, 0, 0]]) >>> r_sun = np.array([[1e11, 0, 0]]) >>> sun_shine(r_sat, r_earth, r_sun, radius=0.4, albedo=0.20, albedo_front=0.05, area_panels=100) {'sun_bus': array([...]), 'sun_panels': array([...])} """ # https://amostech.com/TechnicalPapers/2013/POSTER/COGNION.pdf phase_angle = get_angle(r_sun, r_sat, r_earth) # angle from Sun to object to Earth r_earth_sat = np.linalg.norm(r_sat - r_earth, axis=-1) # Earth is the observer. flux_front = np.zeros_like(phase_angle) flux_front[phase_angle < np.pi / 2] = albedo_front * area_panels / (np.pi * r_earth_sat[phase_angle < np.pi / 2]**2) * np.cos(phase_angle[phase_angle < np.pi / 2]) # Fraction of Sun light scattered off front of the solar panels - which are assumed to be always facing the Sun. Angle: Sun - Sat - Observer flux_bus = 2 / 3 * albedo * radius**2 / (np.pi * (r_earth_sat)**2) * (np.sin(phase_angle) + (np.pi - phase_angle) * np.cos(phase_angle)) # Fraction of light reflected off satellite from Sun return {'sun_bus': flux_bus, 'sun_panels': flux_front}
[docs] def calc_M_v(r_sat, r_earth, r_sun, r_moon=False, radius=0.4, albedo=0.20, sun_Mag=4.80, albedo_earth=0.30, albedo_moon=0.12, albedo_back=0.50, albedo_front=0.05, area_panels=100, return_components=False): """ Calculate the apparent magnitude (M_v) of a satellite due to reflections from the Sun, Earth, and optionally the Moon. This function computes the apparent magnitude of a satellite based on its reflected light from the Sun, Earth, and optionally the Moon. It uses separate functions to calculate the flux contributions from each of these sources and combines them to determine the overall apparent magnitude. Parameters: ---------- r_sat : (n, 3) numpy.ndarray Position of the satellite in meters. r_earth : (n, 3) numpy.ndarray Position of the Earth in meters. r_sun : (n, 3) numpy.ndarray Position of the Sun in meters. r_moon : (n, 3) numpy.ndarray or False, optional Position of the Moon in meters. If False, the Moon's contribution is ignored (default is False). radius : float, optional Radius of the satellite in meters (default is 0.4 m). albedo : float, optional Albedo of the satellite's surface, representing its reflectivity (default is 0.20). sun_Mag : float, optional Solar magnitude (apparent magnitude of the Sun) used in magnitude calculations (default is 4.80). albedo_earth : float, optional Albedo of the Earth, representing its reflectivity (default is 0.30). albedo_moon : float, optional Albedo of the Moon, representing its reflectivity (default is 0.12). albedo_back : float, optional Albedo of the back surface of the satellite (default is 0.50). albedo_front : float, optional Albedo of the front surface of the satellite (default is 0.05). area_panels : float, optional Area of the satellite's panels in square meters (default is 100 m^2). return_components : bool, optional If True, returns the magnitude as well as the flux components from the Sun, Earth, and Moon (default is False). Returns: ------- float The apparent magnitude (M_v) of the satellite as observed from Earth. dict, optional If `return_components` is True, a dictionary containing the flux components from the Sun, Earth, and Moon. Notes: ------ - The function uses separate calculations for flux contributions from the Sun, Earth, and Moon: - `sun_shine` calculates the flux from the Sun. - `earth_shine` calculates the flux from the Earth. - `moon_shine` calculates the flux from the Moon (if applicable). - The apparent magnitude is calculated based on the distances between the satellite, Sun, Earth, and optionally the Moon, as well as their respective albedos and other parameters. Example usage: -------------- >>> r_sat = np.array([[1e7, 2e7, 3e7]]) >>> r_earth = np.array([1.496e11, 0, 0]) >>> r_sun = np.array([0, 0, 0]) >>> Mag_v = calc_M_v(r_sat, r_earth, r_sun, return_components=True) >>> Mag_v (15.63, {'sun_bus': 0.1, 'sun_panels': 0.2, 'earth_bus': 0.05, 'earth_panels': 0.1, 'moon_bus': 0.03, 'moon_panels': 0.07}) """ r_sun_sat = np.linalg.norm(r_sat - r_sun, axis=-1) frac_flux_sun = {'sun_bus': 0, 'sun_panels': 0} frac_flux_earth = {'earth_bus': 0, 'earth_panels': 0} frac_flux_moon = {'moon_bus': 0, 'moon_panels': 0} frac_flux_sun = sun_shine(r_sat, r_earth, r_sun, radius, albedo, albedo_front, area_panels) frac_flux_earth = earth_shine(r_sat, r_earth, r_sun, radius, albedo, albedo_earth, albedo_back, area_panels) if r_moon is not False: frac_flux_moon = moon_shine(r_moon, r_sat, r_earth, r_sun, radius, albedo, albedo_moon, albedo_back, albedo_front, area_panels) merged_dict = {**frac_flux_sun, **frac_flux_earth, **frac_flux_moon} total_frac_flux = sum(merged_dict.values()) Mag_v = (2.5 * np.log10((r_sun_sat / (10 * u.Unit('parsec').to(u.Unit('m'))))**2) + sun_Mag) - 2.5 * np.log10(total_frac_flux) if return_components: return Mag_v, merged_dict else: return Mag_v
[docs] def M_v_lambertian(r_sat, times, radius=1.0, albedo=0.20, sun_Mag=4.80, albedo_earth=0.30, albedo_moon=0.12, plot=False): """ Calculate the apparent magnitude (M_v) of a satellite due to reflections from the Sun, Earth, and Moon. This function computes the apparent magnitude of a satellite based on its reflected light from the Sun, Earth, and Moon, using the Lambertian reflection model. It optionally generates plots to visualize the results. Parameters: ---------- r_sat : (n, 3) numpy.ndarray Position of the satellite in meters. times : Time or array_like The times corresponding to the satellite and celestial body positions. Used to obtain the positions of the Sun and Moon. radius : float, optional Radius of the satellite in meters (default is 1.0 m). albedo : float, optional Albedo of the satellite's surface, representing its reflectivity (default is 0.20). sun_Mag : float, optional Solar magnitude (apparent magnitude of the Sun) used in magnitude calculations (default is 4.80). albedo_earth : float, optional Albedo of the Earth, representing its reflectivity (default is 0.30). albedo_moon : float, optional Albedo of the Moon, representing its reflectivity (default is 0.12). plot : bool, optional If True, generates plots to visualize solar phase angle and magnitudes (default is False). Returns: ------- numpy.ndarray The apparent magnitude (M_v) of the satellite as observed from Earth, considering reflections from the Sun, Earth, and Moon. Notes: ------ - The function uses a Lambertian reflection model to compute the fraction of sunlight reflected by the satellite, Earth, and Moon. - The apparent magnitude is calculated based on the distances between the satellite, Sun, Earth, and Moon, as well as their respective albedos. - The function generates four subplots if `plot` is set to True: 1. Solar phase angle of the satellite. 2. Solar magnitude (M_v) of the satellite due to the Sun. 3. Magnitude (M_v) of the satellite due to reflections from the Earth. 4. Magnitude (M_v) of the satellite due to reflections from the Moon. Example usage: -------------- >>> r_sat = np.array([[1e7, 2e7, 3e7]]) >>> times = Time("2024-01-01") >>> M_v = M_v_lambertian(r_sat, times, plot=True) >>> M_v array([15.63]) """ pc_to_m = 3.085677581491367e+16 r_sun = get_body('Sun').position(times).T r_moon = get_body('Moon').position(times).T r_earth = np.zeros_like(r_sun) r_sun_sat = np.linalg.norm(r_sat - r_sun, axis=-1) r_earth_sat = np.linalg.norm(r_sat, axis=-1) r_moon_sat = np.linalg.norm(r_sat - r_moon, axis=-1) sun_angle = get_angle(r_sun, r_sat, r_earth) earth_angle = np.pi - sun_angle moon_phase_angle = get_angle(r_sun, r_moon, r_sat) # Phase of the moon as viewed from the sat. moon_to_earth_angle = get_angle(r_moon, r_sat, r_earth) flux_moon_to_sat = 2 / 3 * albedo_moon * MOON_RADIUS**2 / (np.pi * (r_moon_sat)**2) * (np.sin(moon_phase_angle) + (np.pi - moon_phase_angle) * np.cos(moon_phase_angle)) # Fraction of sunlight reflected from the Moon to satellite flux_earth_to_sat = 2 / 3 * albedo_earth * EARTH_RADIUS**2 / (np.pi * (r_earth_sat)**2) * (np.sin(earth_angle) + (np.pi - earth_angle) * np.cos(earth_angle)) # Fraction of sunlight reflected from the Earth to satellite frac_flux_sun = 2 / 3 * albedo * radius**2 / (np.pi * (r_earth_sat)**2) * (np.sin(sun_angle) + (np.pi - sun_angle) * np.cos(sun_angle)) # Fraction of light reflected off satellite from Sun frac_flux_earth = 2 / 3 * albedo * radius**2 / (np.pi * r_earth_sat**2) * flux_earth_to_sat frac_flux_moon = 2 / 3 * albedo * radius**2 / (np.pi * r_earth_sat**2) * flux_moon_to_sat Mag_v = (2.5 * np.log10((r_sun_sat / (10 * pc_to_m))**2) + sun_Mag) - 2.5 * np.log10(frac_flux_sun + frac_flux_earth + frac_flux_moon) if plot: import matplotlib.pyplot as plt sun_scale = 149597870700.0 * (RGEO / np.max(r_earth_sat) ) * 0.75 color_map ='inferno_r' fig = plt.figure(figsize=(18, 4)) ax = fig.add_subplot(1, 4, 1) ax.scatter(r_earth[:, 0], r_earth[:, 1], c='Blue', s=10) scatter = ax.scatter(r_sat[:, 0] / RGEO, r_sat[:, 1] / RGEO, c=sun_angle, cmap=color_map) colorbar = plt.colorbar(scatter) ax.scatter(r_sun[:, 0] / sun_scale, r_sun[:, 1] / sun_scale, c=plt.cm.Oranges(np.linspace(0.25, 0.75, len(r_sat[:, 0]))), s=10) ax.set_title('Solar Phase') ax.set_xlabel('X [GEO]') ax.set_ylabel('Y [GEO]') ax.axis('equal') ax = fig.add_subplot(1, 4, 2) ax.scatter(r_earth[0], r_earth[1], c='Blue', s=10) scatter = ax.scatter(r_sat[:, 0] / RGEO, r_sat[:, 1] / RGEO, c=(2.5 * np.log10((r_sun_sat / (10 * pc_to_m))**2) + sun_Mag) - 2.5 * np.log10(frac_flux_sun), cmap=color_map) colorbar = plt.colorbar(scatter) ax.scatter(r_sun[:, 0] / sun_scale, r_sun[:, 1] / sun_scale, c=plt.cm.Oranges(np.linspace(0.25, 0.75, len(r_sat[:, 0]))), s=10) ax.set_title('Solar M_v') ax.axis('equal') ax = fig.add_subplot(1, 4, 3) ax.scatter(r_earth[:, 0], r_earth[:, 1], c='Blue', s=10) scatter = ax.scatter(r_sat[:, 0] / RGEO, r_sat[:, 1] / RGEO, c=(2.5 * np.log10((r_sun_sat / (10 * pc_to_m))**2) + sun_Mag) - 2.5 * np.log10(frac_flux_earth), cmap=color_map) colorbar = plt.colorbar(scatter) ax.scatter(r_sun[:, 0] / sun_scale, r_sun[:, 1] / sun_scale, c=plt.cm.Oranges(np.linspace(0.25, 0.75, len(r_sat[:, 0]))), s=10) ax.set_title('Earth M_v') ax.axis('equal') ax = fig.add_subplot(1, 4, 4) ax.scatter(r_earth[:, 0], r_earth[:, 1], c='Blue', s=10) scatter = ax.scatter(r_sat[:, 0] / RGEO, r_sat[:, 1] / RGEO, c=(2.5 * np.log10((r_sun_sat / (10 * pc_to_m))**2) + sun_Mag) - 2.5 * np.log10(frac_flux_moon), cmap=color_map) ax.scatter(r_moon[:, 0] / RGEO, r_moon[:, 1] / RGEO, c=plt.cm.Greys(np.linspace(0.5, 1, len(r_sat[:, 0]))), s=5) colorbar = plt.colorbar(scatter) ax.set_title('Lunar M_v') ax.axis('equal') plt.show() return Mag_v
[docs] def calc_gamma(r, t): """ Calculate the gamma angle between position and velocity vectors in the ITRF frame. Parameters: r (numpy.ndarray): The position vectors in the GCRF frame, shaped (n, 3), where n is the number of vectors. t (numpy.ndarray or astropy.time.Time): The times corresponding to the position vectors. Can be an array of GPS seconds or an Astropy Time object. Returns: numpy.ndarray: An array of gamma angles in degrees between the position and velocity vectors for each time point. Notes: - This function first converts the given position vectors from the GCRF frame to the ITRF frame. - It then calculates the angle between the position and velocity vectors at each time point in the ITRF frame. - If the input time array is an Astropy Time object, it converts it to GPS time before processing. """ r_itrf, v_itrf = gcrf_to_itrf(r, t, v=True) if isinstance(t[0], Time): t = t.gps gamma = np.degrees(np.apply_along_axis(lambda x: angle_between_vectors(x[:3], x[3:]), 1, np.concatenate((r_itrf, v_itrf), axis=1))) - 90 return gamma
[docs] def moon_normal_vector(t): """ Calculate the normal vector to the Moon's orbital plane at a given time. This function computes the normal vector to the Moon's orbital plane by taking the cross product of the Moon's position vectors at two different times: the given time `t` and one week later. The resulting vector is normalized to provide the direction of the orbital plane's normal. Parameters: ---------- t : Time The time at which to calculate the normal vector to the Moon's orbital plane. Returns: ------- numpy.ndarray A 3-element array representing the normal vector to the Moon's orbital plane at the given time. The vector is normalized to have a unit length. Notes: ------ - The function assumes a circular orbit for the Moon and uses a time step of one week (604800 seconds) to calculate the normal vector. - The normal vector is perpendicular to the plane defined by the Moon's position vectors at the given time and one week later. Example usage: -------------- >>> t = Time("2024-01-01") >>> normal_vector = moon_normal_vector(t) >>> normal_vector array([-0.093, 0.014, 0.995]) """ r = get_body("moon").position(t).T r_random = get_body("moon").position(t.gps + 604800).T return np.cross(r, r_random) / np.linalg.norm(r, axis=-1)
[docs] def lunar_lagrange_points(t): """ Calculate the positions of the lunar Lagrange points in the GCRF frame at a given time. This function computes the positions of the five Lagrange points (L1, L2, L3, L4, and L5) in the Earth-Moon system at a specific time `t`. It considers the positions of the Earth and Moon and uses the orbital period of the Moon to find the Lagrange points. Parameters: ---------- t : Time The time at which to calculate the Lagrange points. The position of the Moon at this time is used to compute the Lagrange points. Returns: ------- dict A dictionary containing the coordinates of the five Lagrange points: - "L1": Position of the first Lagrange point between the Earth and the Moon. - "L2": Position of the second Lagrange point beyond the Moon. - "L3": Position of the third Lagrange point directly opposite the Moon, relative to the Earth. - "L4": Position of the fourth Lagrange point, calculated as the Moon's position offset by one-sixth of the lunar period. - "L5": Position of the fifth Lagrange point, calculated as the Moon's position offset by negative one-sixth of the lunar period. Notes: ------ - The function assumes a circular orbit for the Moon. - The lunar period used is approximately 2.36 million seconds. - The gravitational parameters of the Earth and Moon are denoted as `EARTH_MU` and `MOON_MU`, respectively. Example usage: -------------- >>> t = Time("2024-01-01") >>> lagrange_points = lunar_lagrange_points(t) >>> lagrange_points["L1"] array([1.02e6, 0.0, 0.0]) >>> lagrange_points["L4"] array([1.5e6, 1.5e6, 0.0]) """ r = get_body("moon").position(t).T d = np.linalg.norm(r) # Distance between Earth and Moon unit_vector_moon = r / np.linalg.norm(r, axis=-1) # plane_vector = np.cross(r, r_random) lunar_period_seconds = 2.3605915968e6 # Coefficients of the quadratic equation a = EARTH_MU - MOON_MU b = 2 * MOON_MU * d c = -MOON_MU * d**2 # Solve the quadratic equation discriminant = b**2 - 4*a*c if discriminant >= 0: L1_from_moon = (-b - np.sqrt(discriminant)) / (2*a) * unit_vector_moon L2_from_moon = (-b + np.sqrt(discriminant)) / (2*a) * unit_vector_moon else: print("Discriminate is less than 0! THAT'S WEIRD FIX IT.") L1_from_moon = None L2_from_moon = None return { "L1": L1_from_moon + r, "L2": L2_from_moon + r, "L3": -r, "L4": get_body("moon").position(t.gps + lunar_period_seconds / 6).T, "L5": get_body("moon").position(t.gps - lunar_period_seconds / 6).T }
[docs] def lunar_lagrange_points_circular(t): """ Calculate the positions of the lunar Lagrange points in the GCRF frame for a given time. This function calculates the positions of the five Lagrange points (L1, L2, L3, L4, and L5) in the Earth-Moon system at a specific time `t`. It accounts for the rotation of the Moon's orbit around the Earth, providing the positions in a circular approximation of the Earth-Moon system. Parameters: ---------- t : Time The time at which to calculate the Lagrange points. The position of the Moon at this time is used to compute the Lagrange points. Returns: ------- dict A dictionary containing the coordinates of the five Lagrange points: - "L1": Position of the first Lagrange point between the Earth and the Moon. - "L2": Position of the second Lagrange point beyond the Moon. - "L3": Position of the third Lagrange point directly opposite the Moon, relative to the Earth. - "L4": Position of the fourth Lagrange point, forming an equilateral triangle with the Earth and Moon. - "L5": Position of the fifth Lagrange point, forming an equilateral triangle with the Earth and Moon, but on the opposite side. Notes: ------ - The function assumes a circular orbit for the Moon and uses a rotation matrix to align the z-axis with the Moon's normal vector. - The positions of L4 and L5 are calculated using a rotation matrix to align with the Moon's orientation. - The gravitational parameters of the Earth and Moon are denoted as `EARTH_MU` and `MOON_MU`, respectively. Example usage: -------------- >>> t = Time("2024-01-01") >>> lagrange_points = lunar_lagrange_points_circular(t) >>> lagrange_points["L1"] array([1.02e6, 0.0, 0.0]) >>> lagrange_points["L4"] array([1.5e6, 1.5e6, 0.0]) """ r = get_body("moon").position(t).T d = np.linalg.norm(r) # Distance between Earth and Moon unit_vector_moon = r / np.linalg.norm(r, axis=-1) # Coefficients of the quadratic equation a = EARTH_MU - MOON_MU b = 2 * MOON_MU * d c = -MOON_MU * d**2 # Solve the quadratic equation discriminant = b**2 - 4*a*c if discriminant >= 0: L1_from_moon = (-b - np.sqrt(discriminant)) / (2*a) * unit_vector_moon L2_from_moon = (-b + np.sqrt(discriminant)) / (2*a) * unit_vector_moon else: print("Discriminate is less than 0! THAT'S WEIRD FIX IT.") L1_from_moon = None L2_from_moon = None # L45 # Create the rotation matrix to align z-axis with the normal vector normal_vector = moon_normal_vector(t) rotation_matrix = rotation_matrix_from_vectors(np.array([0, 0, 1]), normal_vector) theta = np.radians(60) + np.arctan2(unit_vector_moon[1], unit_vector_moon[0]) L4 = np.vstack((d * np.cos(theta), d * np.sin(theta), np.zeros_like(theta))).T L4 = np.squeeze(L4 @ rotation_matrix.T) theta = -np.radians(60) + np.arctan2(unit_vector_moon[1], unit_vector_moon[0]) L5 = np.vstack((d * np.cos(theta), d * np.sin(theta), np.zeros_like(theta))).T L5 = np.squeeze(L5 @ rotation_matrix.T) return { "L1": L1_from_moon + r, "L2": L2_from_moon + r, "L3": -r, "L4": L4, "L5": L5 }
[docs] def lagrange_points_lunar_frame(): """ Calculate the positions of the lunar Lagrange points in the lunar frame, This frame is defined by the coordinate transformation in utils.py gcrf_to_lunar(). Returns: ------- dict A dictionary containing the coordinates of the five Lagrange points: - "L1": Position of the first Lagrange point between the Earth and the Moon. - "L2": Position of the second Lagrange point beyond the Moon. - "L3": Position of the third Lagrange point directly opposite the Moon, relative to the Earth. - "L4": Position of the fourth Lagrange point, forming an equilateral triangle with the Earth and Moon. - "L5": Position of the fifth Lagrange point, forming an equilateral triangle with the Earth and Moon, but on the opposite side. Notes: ------ - The function assumes that the Earth and Moon are the two primary bodies, with the Earth-Moon distance denoted as `LD`. - The gravitational parameters of the Earth and Moon are denoted as `EARTH_MU` and `MOON_MU`, respectively. - The positions of L4 and L5 are calculated using the fact that these points form an equilateral triangle with the Earth and Moon. Example usage: -------------- >>> lagrange_points = lagrange_points_lunar_frame() >>> lagrange_points["L1"] array([1.01e6, 0.0, 0.0]) >>> lagrange_points["L4"] array([1.5e6, 1.5e6, 0.0]) """ r = np.array([LD / RGEO, 0, 0]) d = np.linalg.norm(r) # Distance between Earth and Moon unit_vector_moon = r / np.linalg.norm(r, axis=-1) # Coefficients of the quadratic equation a = EARTH_MU - MOON_MU b = 2 * MOON_MU * d c = -MOON_MU * d**2 # Solve the quadratic equation discriminant = b**2 - 4*a*c if discriminant >= 0: L1_from_moon = (-b - np.sqrt(discriminant)) / (2*a) * unit_vector_moon L2_from_moon = (-b + np.sqrt(discriminant)) / (2*a) * unit_vector_moon else: print("Discriminate is less than 0! THAT'S WEIRD FIX IT.") L1_from_moon = None L2_from_moon = None # L45 theta = np.radians(60) L4 = np.squeeze(np.vstack((d * np.cos(theta), d * np.sin(theta), np.zeros_like(theta))).T) L5 = np.squeeze(np.vstack((d * np.cos(-theta), d * np.sin(-theta), np.zeros_like(-theta))).T) return { "L1": L1_from_moon + r, "L2": L2_from_moon + r, "L3": -r, "L4": L4, "L5": L5 }