11294204788

Committed 11 Oct 2024 02:30PM UTC coverage: 66.032%. First build

Build # 11294204788

Build Type

Pull #395

github

Committed by

web-flow

Commit Message

Merge 79a05365c into 18ce85989

Pull Request Pull Request #395: np.Inf and np.array(obj, copy = False) changed

Run Details

86 of 104 new or added lines in 26 files covered. (82.69%)

9543 of 14452 relevant lines covered (66.03%)

0.66 hits per line

Source File
Press 'n' to go to next uncovered line, 'b' for previous

71.11

/EXOSIMS/Observatory/SotoStarshade.py

from EXOSIMS.Observatory.ObservatoryL2Halo import ObservatoryL2Halo
from EXOSIMS.Prototypes.TargetList import TargetList
import numpy as np
import astropy.units as u
from scipy.integrate import solve_bvp
import astropy.constants as const
import hashlib
import scipy.optimize as optimize
import scipy.interpolate as interp
import time
import os
import pickle

EPS = np.finfo(float).eps


class SotoStarshade(ObservatoryL2Halo):
    """StarShade Observatory class
    This class is implemented at L2 and contains all variables, functions,
    and integrators to calculate occulter dynamics.
    """

    def __init__(self, orbit_datapath=None, f_nStars=10, **specs):

        ObservatoryL2Halo.__init__(self, **specs)
        self.f_nStars = int(f_nStars)

        # instantiating fake star catalog, used to generate good dVmap
        fTL = TargetList(
            **{
                "ntargs": self.f_nStars,
                "modules": {
                    "StarCatalog": "FakeCatalog",
                    "TargetList": " ",
                    "OpticalSystem": "Nemati",
                    "ZodiacalLight": "Stark",
                    "PostProcessing": " ",
                    "Completeness": " ",
                    "BackgroundSources": "GalaxiesFaintStars",
                    "PlanetPhysicalModel": " ",
                    "PlanetPopulation": "KeplerLike1",
                },
                "scienceInstruments": [{"name": "imager"}],
                "starlightSuppressionSystems": [{"name": "HLC-565"}],
            }
        )

        f_sInds = np.arange(0, fTL.nStars)
        dV, ang, dt = self.generate_dVMap(fTL, 0, f_sInds, self.equinox[0])

        # pick out unique angle values
        ang, unq = np.unique(ang, return_index=True)
        dV = dV[:, unq]

        # create dV 2D interpolant -- assuming further that x, y are scalars.
        # Matching linear interpolation.
        r = interp.RectBivariateSpline(dt, ang, dV, kx=1, ky=1)
        self.dV_interp = lambda x, y: r(x, y)[0]

    # self.dV_interp = interp.interp2d(dt, ang, dV.T)

    def generate_dVMap(self, TL, old_sInd, sInds, currentTime):
        """Creates dV map for an occulter slewing between targets.

        This method returns a 2D array of the dV needed for an occulter
        to slew between all the different stars on the target list. The dV
        map is calculated relative to a reference star and the stars are
        ordered by their angular separation from the reference star (X-axis).
        The Y-axis represents the time of flight ("slew time") for a
        trajectory between two stars.

        Args:
            TL (TargetList module):
                TargetList class object
            old_sInd (integer):
                Integer index of the last star of interest
            sInds (integer ndarray):
                Integer indices of the stars of interest
            currentTime (astropy Time array):
                Current absolute mission time in MJD

        Returns:
            tuple:
            dVMap (float ndarray):
                Map of dV needed to transfer from a reference star to another.
                Each ordered pair (psi,t) of the dV map corresponds to a
                trajectory to a star an angular distance psi away with flight
                time of t. units of (m/s)
            angles (float ndarray):
                Range of angles (in deg) used in dVMap as the X-axis
            dt (float ndarray):
                Range of slew times (in days) used in dVMap as the Y-axis
        """

        # generating hash name
        filename = "dVMap_"
        extstr = ""
        extstr += "%s: " % "occulterSep" + str(getattr(self, "occulterSep")) + " "
        extstr += "%s: " % "period_halo" + str(getattr(self, "period_halo")) + " "
        extstr += "%s: " % "f_nStars" + str(getattr(self, "f_nStars")) + " "
        extstr += "%s: " % "occ_dtmin" + str(getattr(self, "occ_dtmin")) + " "
        extstr += "%s: " % "occ_dtmax" + str(getattr(self, "occ_dtmax")) + " "
        ext = hashlib.md5(extstr.encode("utf-8")).hexdigest()
        filename += ext
        dVpath = os.path.join(self.cachedir, filename + ".dVmap")

        # initiating slew Times for starshade
        dt = np.arange(self.occ_dtmin.value, self.occ_dtmax.value, 1)

        # angular separation of stars in target list from old_sInd
        ang = self.star_angularSep(TL, old_sInd, sInds, currentTime)
        sInd_sorted = np.argsort(ang)
        angles = ang[sInd_sorted].to("deg").value

        # initializing dV map
        dVMap = np.zeros([len(dt), len(sInds)])

        # checking to see if map exists or needs to be calculated
        if os.path.exists(dVpath):
            # dV map already exists for given parameters
            self.vprint("Loading cached Starshade dV map file from %s" % dVpath)
            try:
                with open(dVpath, "rb") as ff:
                    A = pickle.load(ff)
            except UnicodeDecodeError:
                with open(dVpath, "rb") as ff:
                    A = pickle.load(ff, encoding="latin1")
            self.vprint("Starshade dV Map loaded from cache.")
            dVMap = A["dVMap"]
        else:
            self.vprint('Cached Starshade dV map file not found at "%s".' % dVpath)
            # looping over all target list and desired slew times to generate dV map
            self.vprint("Starting dV calculations for %s stars." % TL.nStars)
            tic = time.perf_counter()
            for i in range(len(dt)):
                dVMap[i, :] = self.impulsiveSlew_dV(
                    dt[i], TL, old_sInd, sInd_sorted, currentTime
                )  # sorted
                if not i % 5:
                    self.vprint("   [%s / %s] completed." % (i, len(dt)))
            toc = time.perf_counter()
            B = {"dVMap": dVMap}
            with open(dVpath, "wb") as ff:
                pickle.dump(B, ff)
            self.vprint("dV map computation completed in %s seconds." % (toc - tic))
            self.vprint("dV Map array stored in %r" % dVpath)

        return dVMap, angles, dt

    def boundary_conditions(self, rA, rB):
        """Creates boundary conditions for solving a boundary value problem

        This method returns the boundary conditions for the starshade transfer
        trajectory between the lines of sight of two different stars. Point A
        corresponds to the starshade alignment with star A; Point B, with star B.

        Args:
            rA (float 1x3 ndarray):
                Starshade position vector aligned with current star of interest
            rB (float 1x3 ndarray):
                Starshade position vector aligned with next star of interest

        Returns:
            float 1x6 ndarray:
                Star position vector in rotating frame in units of AU
        """

        BC1 = rA[0] - self.rA[0]
        BC2 = rA[1] - self.rA[1]
        BC3 = rA[2] - self.rA[2]

        BC4 = rB[0] - self.rB[0]
        BC5 = rB[1] - self.rB[1]
        BC6 = rB[2] - self.rB[2]

        BC = np.array([BC1, BC2, BC3, BC4, BC5, BC6])

        return BC

    def send_it(self, TL, nA, nB, tA, tB):
        """Solves boundary value problem between starshade star alignments

        This method solves the boundary value problem for starshade star alignments
        with two given stars at times tA and tB. It uses scipy's solve_bvp method.

        Args:
            TL (float 1x3 ndarray):
                TargetList class object
            nA (integer):
                Integer index of the current star of interest
            nB (integer):
                Integer index of the next star of interest
            tA (astropy Time array):
                Current absolute mission time in MJD
            tB (astropy Time array):
                Absolute mission time for next star alignment in MJD

        Returns:
            float nx6 ndarray:
                State vectors in rotating frame in normalized units
        """

        angle, uA, uB, r_tscp = self.lookVectors(TL, nA, nB, tA, tB)

        vA = self.haloVelocity(tA)[0].value / (2 * np.pi)
        vB = self.haloVelocity(tB)[0].value / (2 * np.pi)

        # position vector of occulter in heliocentric frame
        self.rA = uA * self.occulterSep.to("au").value + r_tscp[0]
        self.rB = uB * self.occulterSep.to("au").value + r_tscp[-1]

        a = ((np.mod(tA.value, self.equinox[0].value) * u.d)).to("yr").value * (
            2 * np.pi
        )
        b = ((np.mod(tB.value, self.equinox[0].value) * u.d)).to("yr").value * (
            2 * np.pi
        )

        # running shooting algorithm
        t = np.linspace(a, b, 2)

        sG = np.array(
            [
                [self.rA[0], self.rB[0]],
                [self.rA[1], self.rB[1]],
                [self.rA[2], self.rB[2]],
                [vA[0], vB[0]],
                [vA[1], vB[1]],
                [vA[2], vB[2]],
            ]
        )

        sol = solve_bvp(
            self.equationsOfMotion_CRTBP, self.boundary_conditions, t, sG, tol=1e-10
        )

        s = sol.y.T
        t_s = sol.x

        return s, t_s

    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.

        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
        """

        if old_sInd is None:
            dV = np.zeros(slewTimes.shape)
        else:
            dV = np.zeros(slewTimes.shape)
            badSlews_i, badSlew_j = np.where(slewTimes.value < self.occ_dtmin.value)
            for i in range(len(sInds)):
                for t in range(len(slewTimes.T)):
                    dV[i, t] = self.dV_interp(slewTimes[i, t], sd[i].to("deg"))
            dV[badSlews_i, badSlew_j] = np.inf
        return dV * u.m / u.s

    def impulsiveSlew_dV(self, dt, TL, nA, N, tA):
        """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.

        Args:
            dt (float 1x1 ndarray):
                Number of days corresponding to starshade slew time
            TL (float 1x3 ndarray):
                TargetList class object
            nA (integer):
                Integer index of the current star of interest
            N  (integer):
                Integer index of the next star(s) of interest
            tA (astropy Time array):
                Current absolute mission time in MJD

        Returns:
            float nx6 ndarray:
                State vectors in rotating frame in normalized units
        """

        if dt.shape:
            dt = dt[0]

        if nA is None:
            dV = np.zeros(len(N))
        else:
            # if only calculating one trajectory, this allows loop to run
            if N.size == 1:
                N = np.array([N])

            # time to reach star B's line of sight
            tB = tA + dt * u.d

            # initializing arrays for BVP state solutions
            sol_slew = np.zeros([2, len(N), 6])
            t_sol = np.zeros([2, len(N)])
            for x in range(len(N)):
                # simulating slew trajectory from star A at tA to star B at tB
                sol, t = self.send_it(TL, nA, N[x], tA, tB)
                sol_slew[:, x, :] = np.array([sol[0], sol[-1]])
                t_sol[:, x] = np.array([t[0], t[-1]])

            # starshade velocities at both endpoints of the slew trajectory
            r_slewA = sol_slew[0, :, 0:3]
            r_slewB = sol_slew[-1, :, 0:3]
            v_slewA = sol_slew[0, :, 3:6]
            v_slewB = sol_slew[-1, :, 3:6]

            if len(N) == 1:
                t_slewA = t_sol[0]
                t_slewB = t_sol[1]
            else:
                t_slewA = t_sol[0, 0]
                t_slewB = t_sol[1, 1]

            r_haloA = (self.haloPosition(tA) + self.L2_dist * np.array([1, 0, 0]))[
                0
            ] / u.AU
            r_haloB = (self.haloPosition(tB) + self.L2_dist * np.array([1, 0, 0]))[
                0
            ] / u.AU

            v_haloA = self.haloVelocity(tA)[0] / u.AU * u.year / (2 * np.pi)
            v_haloB = self.haloVelocity(tB)[0] / u.AU * u.year / (2 * np.pi)

            dvA = self.rot2inertV(r_slewA, v_slewA, t_slewA) - self.rot2inertV(
                r_haloA.value, v_haloA.value, t_slewA
            )
            dvB = self.rot2inertV(r_slewB, v_slewB, t_slewB) - self.rot2inertV(
                r_haloB.value, v_haloB.value, t_slewB
            )

            if len(dvA) == 1:
                dV = np.linalg.norm(dvA) * u.AU / u.year * (2 * np.pi) + np.linalg.norm(
                    dvB
                ) * u.AU / u.year * (2 * np.pi)
            else:
                dV = np.linalg.norm(dvA, axis=1) * u.AU / u.year * (
                    2 * np.pi
                ) + np.linalg.norm(dvB, axis=1) * u.AU / u.year * (2 * np.pi)

        return dV.to("m/s")

    def minimize_slewTimes(self, TL, nA, nB, tA):
        """Minimizes the slew time for a starshade transferring to a new star
        line of sight

        This method uses scipy's optimization module to minimize the slew time for
        a starshade transferring between one star's line of sight to another's under
        the constraint that the total change in velocity cannot exceed more than a
        certain percentage of the total fuel on board the starshade.

        Args:
            TL (float 1x3 ndarray):
                TargetList class object
            nA (integer):
                Integer index of the current star of interest
            nB (integer):
                Integer index of the next star of interest
            tA (astropy Time array):
                Current absolute mission time in MJD

        Returns:
            tuple:
                opt_slewTime (float):
                    Optimal slew time in days for starshade transfer to a new
                    line of sight
                opt_dV (float):
                    Optimal total change in velocity in m/s for starshade
                    line of sight transfer

        """

        def slewTime_objFun(dt):
            if dt.shape:
                dt = dt[0]

            return dt

        def slewTime_constraints(dt, TL, nA, nB, tA):
            dV = self.calculate_dV(dt, TL, nA, nB, tA)
            dV_max = self.dVmax

            return (dV_max - dV).value, dt - 1

        dt_guess = 20
        Tol = 1e-3

        t0 = [dt_guess]

        res = optimize.minimize(
            slewTime_objFun,
            t0,
            method="COBYLA",
            constraints={
                "type": "ineq",
                "fun": slewTime_constraints,
                "args": ([TL, nA, nB, tA]),
            },
            tol=Tol,
            options={"disp": False},
        )

        opt_slewTime = res.x
        opt_dV = self.calculate_dV(opt_slewTime, TL, nA, nB, tA)

        return opt_slewTime, opt_dV.value

    def minimize_fuelUsage(self, TL, nA, nB, tA):
        """Minimizes the fuel usage of a starshade transferring to a new star
        line of sight

        This method uses scipy's optimization module to minimize the fuel usage for
        a starshade transferring between one star's line of sight to another's. The
        total slew time for the transfer is bounded with some dt_min and dt_max.

        Args:
            TL (float 1x3 ndarray):
                TargetList class object
            nA (integer):
                Integer index of the current star of interest
            nB (integer):
                Integer index of the next star of interest
            tA (astropy Time array):
                Current absolute mission time in MJD

        Returns:
            tuple:
                opt_slewTime (float):
                    Optimal slew time in days for starshade transfer to a
                    new line of sight
                opt_dV (float):
                    Optimal total change in velocity in m/s for starshade
                    line of sight transfer

        """

        def fuelUsage_objFun(dt, TL, nA, N, tA):
            dV = self.calculate_dV(dt, TL, nA, N, tA)
            return dV.value

        def fuelUsage_constraints(dt, dt_min, dt_max):
            return dt_max - dt, dt - dt_min

        dt_guess = 20
        dt_min = 1
        dt_max = 45
        Tol = 1e-5

        t0 = [dt_guess]

        res = optimize.minimize(
            fuelUsage_objFun,
            t0,
            method="COBYLA",
            args=(TL, nA, nB, tA),
            constraints={
                "type": "ineq",
                "fun": fuelUsage_constraints,
                "args": ([dt_min, dt_max]),
            },
            tol=Tol,
            options={"disp": False},
        )
        opt_slewTime = res.x
        opt_dV = res.fun

        return opt_slewTime, opt_dV

    def calculate_slewTimes(self, TL, old_sInd, sInds, sd, obsTimes, tmpCurrentTimeAbs):
        """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.
            tmpCurrentTimeAbs (~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
        """

        if old_sInd is None:
            slewTimes = np.zeros(len(sInds)) * u.d
        else:
            obsTimeRangeNorm = (obsTimes - tmpCurrentTimeAbs).value
            slewTimes = obsTimeRangeNorm[0, :] * u.d

        return slewTimes

    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 Quantity):
                Time to transfer to new star line of sight in units of days
            sInd (integer):
                Integer index of the star of interest
            sd (astropy Quantity):
                Angular separation between stars in rad
            dV (astropy Quantity):
                Delta-V used to transfer to new star line of sight in units of m/s

        Returns:
            dict:
                Design Reference Mission dicitonary, contains the results of one
                complete observation (detection and characterization)
        """

        DRM["slew_time"] = slewTimes.to("day")
        DRM["slew_angle"] = sd.to("deg")

        dV = dV.to("m/s")
        slew_mass_used = self.scMass * (
            1 - np.exp(-dV.value / (self.slewIsp.value * const.g0.value))
        )
        DRM["slew_dV"] = dV
        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

1	from EXOSIMS.Observatory.ObservatoryL2Halo import ObservatoryL2Halo	1✔
2	from EXOSIMS.Prototypes.TargetList import TargetList	1✔
3	import numpy as np	1✔
4	import astropy.units as u	1✔
5	from scipy.integrate import solve_bvp	1✔
6	import astropy.constants as const	1✔
7	import hashlib	1✔
8	import scipy.optimize as optimize	1✔
9	import scipy.interpolate as interp	1✔
10	import time	1✔
11	import os	1✔
12	import pickle	1✔
13
14	EPS = np.finfo(float).eps	1✔
15
16
17	class SotoStarshade(ObservatoryL2Halo):	1✔
18	"""StarShade Observatory class
19	This class is implemented at L2 and contains all variables, functions,
20	and integrators to calculate occulter dynamics.
21	"""
22
23	def __init__(self, orbit_datapath=None, f_nStars=10, **specs):	1✔
24
25	ObservatoryL2Halo.__init__(self, **specs)	1✔
26	self.f_nStars = int(f_nStars)	1✔
27
28	# instantiating fake star catalog, used to generate good dVmap
29	fTL = TargetList(	1✔
30	**{
31	"ntargs": self.f_nStars,
32	"modules": {
33	"StarCatalog": "FakeCatalog",
34	"TargetList": " ",
35	"OpticalSystem": "Nemati",
36	"ZodiacalLight": "Stark",
37	"PostProcessing": " ",
38	"Completeness": " ",
39	"BackgroundSources": "GalaxiesFaintStars",
40	"PlanetPhysicalModel": " ",
41	"PlanetPopulation": "KeplerLike1",
42	},
43	"scienceInstruments": [{"name": "imager"}],
44	"starlightSuppressionSystems": [{"name": "HLC-565"}],
45	}
46	)
47
48	f_sInds = np.arange(0, fTL.nStars)	1✔
49	dV, ang, dt = self.generate_dVMap(fTL, 0, f_sInds, self.equinox[0])	1✔
50
51	# pick out unique angle values
52	ang, unq = np.unique(ang, return_index=True)	1✔
53	dV = dV[:, unq]	1✔
54
55	# create dV 2D interpolant -- assuming further that x, y are scalars.
56	# Matching linear interpolation.
57	r = interp.RectBivariateSpline(dt, ang, dV, kx=1, ky=1)	1✔
58	self.dV_interp = lambda x, y: r(x, y)[0]	1✔
59
60	# self.dV_interp = interp.interp2d(dt, ang, dV.T)
61
62	def generate_dVMap(self, TL, old_sInd, sInds, currentTime):	1✔
63	"""Creates dV map for an occulter slewing between targets.
64
65	This method returns a 2D array of the dV needed for an occulter
66	to slew between all the different stars on the target list. The dV
67	map is calculated relative to a reference star and the stars are
68	ordered by their angular separation from the reference star (X-axis).
69	The Y-axis represents the time of flight ("slew time") for a
70	trajectory between two stars.
71
72	Args:
73	TL (TargetList module):
74	TargetList class object
75	old_sInd (integer):
76	Integer index of the last star of interest
77	sInds (integer ndarray):
78	Integer indices of the stars of interest
79	currentTime (astropy Time array):
80	Current absolute mission time in MJD
81
82	Returns:
83	tuple:
84	dVMap (float ndarray):
85	Map of dV needed to transfer from a reference star to another.
86	Each ordered pair (psi,t) of the dV map corresponds to a
87	trajectory to a star an angular distance psi away with flight
88	time of t. units of (m/s)
89	angles (float ndarray):
90	Range of angles (in deg) used in dVMap as the X-axis
91	dt (float ndarray):
92	Range of slew times (in days) used in dVMap as the Y-axis
93	"""
94
95	# generating hash name
96	filename = "dVMap_"	1✔
97	extstr = ""	1✔
98	extstr += "%s: " % "occulterSep" + str(getattr(self, "occulterSep")) + " "	1✔
99	extstr += "%s: " % "period_halo" + str(getattr(self, "period_halo")) + " "	1✔
100	extstr += "%s: " % "f_nStars" + str(getattr(self, "f_nStars")) + " "	1✔
101	extstr += "%s: " % "occ_dtmin" + str(getattr(self, "occ_dtmin")) + " "	1✔
102	extstr += "%s: " % "occ_dtmax" + str(getattr(self, "occ_dtmax")) + " "	1✔
103	ext = hashlib.md5(extstr.encode("utf-8")).hexdigest()	1✔
104	filename += ext	1✔
105	dVpath = os.path.join(self.cachedir, filename + ".dVmap")	1✔
106
107	# initiating slew Times for starshade
108	dt = np.arange(self.occ_dtmin.value, self.occ_dtmax.value, 1)	1✔
109
110	# angular separation of stars in target list from old_sInd
111	ang = self.star_angularSep(TL, old_sInd, sInds, currentTime)	1✔
112	sInd_sorted = np.argsort(ang)	1✔
113	angles = ang[sInd_sorted].to("deg").value	1✔
114
115	# initializing dV map
116	dVMap = np.zeros([len(dt), len(sInds)])	1✔
117
118	# checking to see if map exists or needs to be calculated
119	if os.path.exists(dVpath):	1✔
120	# dV map already exists for given parameters
121	self.vprint("Loading cached Starshade dV map file from %s" % dVpath)	1✔
122	try:	1✔
123	with open(dVpath, "rb") as ff:	1✔
124	A = pickle.load(ff)	1✔
125	except UnicodeDecodeError:	×
126	with open(dVpath, "rb") as ff:	×
127	A = pickle.load(ff, encoding="latin1")	×
128	self.vprint("Starshade dV Map loaded from cache.")	1✔
129	dVMap = A["dVMap"]	1✔
130	else:
131	self.vprint('Cached Starshade dV map file not found at "%s".' % dVpath)	1✔
132	# looping over all target list and desired slew times to generate dV map
133	self.vprint("Starting dV calculations for %s stars." % TL.nStars)	1✔
134	tic = time.perf_counter()	1✔
135	for i in range(len(dt)):	1✔
136	dVMap[i, :] = self.impulsiveSlew_dV(	1✔
137	dt[i], TL, old_sInd, sInd_sorted, currentTime
138	) # sorted
139	if not i % 5:	1✔
140	self.vprint(" [%s / %s] completed." % (i, len(dt)))	1✔
141	toc = time.perf_counter()	1✔
142	B = {"dVMap": dVMap}	1✔
143	with open(dVpath, "wb") as ff:	1✔
144	pickle.dump(B, ff)	1✔
145	self.vprint("dV map computation completed in %s seconds." % (toc - tic))	1✔
146	self.vprint("dV Map array stored in %r" % dVpath)	1✔
147
148	return dVMap, angles, dt	1✔
149
150	def boundary_conditions(self, rA, rB):	1✔
151	"""Creates boundary conditions for solving a boundary value problem
152
153	This method returns the boundary conditions for the starshade transfer
154	trajectory between the lines of sight of two different stars. Point A
155	corresponds to the starshade alignment with star A; Point B, with star B.
156
157	Args:
158	rA (float 1x3 ndarray):
159	Starshade position vector aligned with current star of interest
160	rB (float 1x3 ndarray):
161	Starshade position vector aligned with next star of interest
162
163	Returns:
164	float 1x6 ndarray:
165	Star position vector in rotating frame in units of AU
166	"""
167
168	BC1 = rA[0] - self.rA[0]	1✔
169	BC2 = rA[1] - self.rA[1]	1✔
170	BC3 = rA[2] - self.rA[2]	1✔
171
172	BC4 = rB[0] - self.rB[0]	1✔
173	BC5 = rB[1] - self.rB[1]	1✔
174	BC6 = rB[2] - self.rB[2]	1✔
175
176	BC = np.array([BC1, BC2, BC3, BC4, BC5, BC6])	1✔
177
178	return BC	1✔
179
180	def send_it(self, TL, nA, nB, tA, tB):	1✔
181	"""Solves boundary value problem between starshade star alignments
182
183	This method solves the boundary value problem for starshade star alignments
184	with two given stars at times tA and tB. It uses scipy's solve_bvp method.
185
186	Args:
187	TL (float 1x3 ndarray):
188	TargetList class object
189	nA (integer):
190	Integer index of the current star of interest
191	nB (integer):
192	Integer index of the next star of interest
193	tA (astropy Time array):
194	Current absolute mission time in MJD
195	tB (astropy Time array):
196	Absolute mission time for next star alignment in MJD
197
198	Returns:
199	float nx6 ndarray:
200	State vectors in rotating frame in normalized units
201	"""
202
203	angle, uA, uB, r_tscp = self.lookVectors(TL, nA, nB, tA, tB)	1✔
204
205	vA = self.haloVelocity(tA)[0].value / (2 * np.pi)	1✔
206	vB = self.haloVelocity(tB)[0].value / (2 * np.pi)	1✔
207
208	# position vector of occulter in heliocentric frame
209	self.rA = uA * self.occulterSep.to("au").value + r_tscp[0]	1✔
210	self.rB = uB * self.occulterSep.to("au").value + r_tscp[-1]	1✔
211
212	a = ((np.mod(tA.value, self.equinox[0].value) * u.d)).to("yr").value * (	1✔
213	2 * np.pi
214	)
215	b = ((np.mod(tB.value, self.equinox[0].value) * u.d)).to("yr").value * (	1✔
216	2 * np.pi
217	)
218
219	# running shooting algorithm
220	t = np.linspace(a, b, 2)	1✔
221
222	sG = np.array(	1✔
223	[
224	[self.rA[0], self.rB[0]],
225	[self.rA[1], self.rB[1]],
226	[self.rA[2], self.rB[2]],
227	[vA[0], vB[0]],
228	[vA[1], vB[1]],
229	[vA[2], vB[2]],
230	]
231	)
232
233	sol = solve_bvp(	1✔
234	self.equationsOfMotion_CRTBP, self.boundary_conditions, t, sG, tol=1e-10
235	)
236
237	s = sol.y.T	1✔
238	t_s = sol.x	1✔
239
240	return s, t_s	1✔
241
242	def calculate_dV(self, TL, old_sInd, sInds, sd, slewTimes, tmpCurrentTimeAbs):	1✔
243	"""Finds the change in velocity needed to transfer to a new star line of sight
244
245	This method sums the total delta-V needed to transfer from one star
246	line of sight to another. It determines the change in velocity to move from
247	one station-keeping orbit to a transfer orbit at the current time, then from
248	the transfer orbit to the next station-keeping orbit at currentTime + dt.
249	Station-keeping orbits are modeled as discrete boundary value problems.
250	This method can handle multiple indeces for the next target stars and calculates
251	the dVs of each trajectory from the same starting star.
252
253	Args:
254	TL (:ref:`TargetList`):
255	TargetList class object
256	old_sInd (int):
257	Index of the current star
258	sInds (~numpy.ndarray(int)):
259	Integer index of the next star(s) of interest
260	sd (~astropy.units.Quantity(~numpy.ndarray(float))):
261	Angular separation between stars in rad
262	slewTimes (~astropy.time.Time(~numpy.ndarray)):
263	Slew times.
264	tmpCurrentTimeAbs (~astropy.time.Time):
265	Current absolute mission time in MJD
266
267	Returns:
268	~astropy.units.Quantity(~numpy.ndarray(float)):
269	Delta-V values in units of length/time
270	"""
271
272	if old_sInd is None:	×
273	dV = np.zeros(slewTimes.shape)	×
274	else:
275	dV = np.zeros(slewTimes.shape)	×
276	badSlews_i, badSlew_j = np.where(slewTimes.value < self.occ_dtmin.value)	×
277	for i in range(len(sInds)):	×
278	for t in range(len(slewTimes.T)):	×
279	dV[i, t] = self.dV_interp(slewTimes[i, t], sd[i].to("deg"))	×
NEW 280	dV[badSlews_i, badSlew_j] = np.inf	×
281	return dV * u.m / u.s	×
282
283	def impulsiveSlew_dV(self, dt, TL, nA, N, tA):	1✔
284	"""Finds the change in velocity needed to transfer to a new star line of sight
285
286	This method sums the total delta-V needed to transfer from one star
287	line of sight to another. It determines the change in velocity to move from
288	one station-keeping orbit to a transfer orbit at the current time, then from
289	the transfer orbit to the next station-keeping orbit at currentTime + dt.
290	Station-keeping orbits are modeled as discrete boundary value problems.
291	This method can handle multiple indeces for the next target stars and calculates
292	the dVs of each trajectory from the same starting star.
293
294	Args:
295	dt (float 1x1 ndarray):
296	Number of days corresponding to starshade slew time
297	TL (float 1x3 ndarray):
298	TargetList class object
299	nA (integer):
300	Integer index of the current star of interest
301	N (integer):
302	Integer index of the next star(s) of interest
303	tA (astropy Time array):
304	Current absolute mission time in MJD
305
306	Returns:
307	float nx6 ndarray:
308	State vectors in rotating frame in normalized units
309	"""
310
311	if dt.shape:	1✔
312	dt = dt[0]	×
313
314	if nA is None:	1✔
315	dV = np.zeros(len(N))	×
316	else:
317	# if only calculating one trajectory, this allows loop to run
318	if N.size == 1:	1✔
319	N = np.array([N])	×
320
321	# time to reach star B's line of sight
322	tB = tA + dt * u.d	1✔
323
324	# initializing arrays for BVP state solutions
325	sol_slew = np.zeros([2, len(N), 6])	1✔
326	t_sol = np.zeros([2, len(N)])	1✔
327	for x in range(len(N)):	1✔
328	# simulating slew trajectory from star A at tA to star B at tB
329	sol, t = self.send_it(TL, nA, N[x], tA, tB)	1✔
330	sol_slew[:, x, :] = np.array([sol[0], sol[-1]])	1✔
331	t_sol[:, x] = np.array([t[0], t[-1]])	1✔
332
333	# starshade velocities at both endpoints of the slew trajectory
334	r_slewA = sol_slew[0, :, 0:3]	1✔
335	r_slewB = sol_slew[-1, :, 0:3]	1✔
336	v_slewA = sol_slew[0, :, 3:6]	1✔
337	v_slewB = sol_slew[-1, :, 3:6]	1✔
338
339	if len(N) == 1:	1✔
340	t_slewA = t_sol[0]	×
341	t_slewB = t_sol[1]	×
342	else:
343	t_slewA = t_sol[0, 0]	1✔
344	t_slewB = t_sol[1, 1]	1✔
345
346	r_haloA = (self.haloPosition(tA) + self.L2_dist * np.array([1, 0, 0]))[	1✔
347	0
348	] / u.AU
349	r_haloB = (self.haloPosition(tB) + self.L2_dist * np.array([1, 0, 0]))[	1✔
350	0
351	] / u.AU
352
353	v_haloA = self.haloVelocity(tA)[0] / u.AU * u.year / (2 * np.pi)	1✔
354	v_haloB = self.haloVelocity(tB)[0] / u.AU * u.year / (2 * np.pi)	1✔
355
356	dvA = self.rot2inertV(r_slewA, v_slewA, t_slewA) - self.rot2inertV(	1✔
357	r_haloA.value, v_haloA.value, t_slewA
358	)
359	dvB = self.rot2inertV(r_slewB, v_slewB, t_slewB) - self.rot2inertV(	1✔
360	r_haloB.value, v_haloB.value, t_slewB
361	)
362
363	if len(dvA) == 1:	1✔
364	dV = np.linalg.norm(dvA) * u.AU / u.year * (2 * np.pi) + np.linalg.norm(	×
365	dvB
366	) * u.AU / u.year * (2 * np.pi)
367	else:
368	dV = np.linalg.norm(dvA, axis=1) * u.AU / u.year * (	1✔
369	2 * np.pi
370	) + np.linalg.norm(dvB, axis=1) * u.AU / u.year * (2 * np.pi)
371
372	return dV.to("m/s")	1✔
373
374	def minimize_slewTimes(self, TL, nA, nB, tA):	1✔
375	"""Minimizes the slew time for a starshade transferring to a new star
376	line of sight
377
378	This method uses scipy's optimization module to minimize the slew time for
379	a starshade transferring between one star's line of sight to another's under
380	the constraint that the total change in velocity cannot exceed more than a
381	certain percentage of the total fuel on board the starshade.
382
383	Args:
384	TL (float 1x3 ndarray):
385	TargetList class object
386	nA (integer):
387	Integer index of the current star of interest
388	nB (integer):
389	Integer index of the next star of interest
390	tA (astropy Time array):
391	Current absolute mission time in MJD
392
393	Returns:
394	tuple:
395	opt_slewTime (float):
396	Optimal slew time in days for starshade transfer to a new
397	line of sight
398	opt_dV (float):
399	Optimal total change in velocity in m/s for starshade
400	line of sight transfer
401
402	"""
403
404	def slewTime_objFun(dt):	×
405	if dt.shape:	×
406	dt = dt[0]	×
407
408	return dt	×
409
410	def slewTime_constraints(dt, TL, nA, nB, tA):	×
411	dV = self.calculate_dV(dt, TL, nA, nB, tA)	×
412	dV_max = self.dVmax	×
413
414	return (dV_max - dV).value, dt - 1	×
415
416	dt_guess = 20	×
417	Tol = 1e-3	×
418
419	t0 = [dt_guess]	×
420
421	res = optimize.minimize(	×
422	slewTime_objFun,
423	t0,
424	method="COBYLA",
425	constraints={
426	"type": "ineq",
427	"fun": slewTime_constraints,
428	"args": ([TL, nA, nB, tA]),
429	},
430	tol=Tol,
431	options={"disp": False},
432	)
433
434	opt_slewTime = res.x	×
435	opt_dV = self.calculate_dV(opt_slewTime, TL, nA, nB, tA)	×
436
437	return opt_slewTime, opt_dV.value	×
438
439	def minimize_fuelUsage(self, TL, nA, nB, tA):	1✔
440	"""Minimizes the fuel usage of a starshade transferring to a new star
441	line of sight
442
443	This method uses scipy's optimization module to minimize the fuel usage for
444	a starshade transferring between one star's line of sight to another's. The
445	total slew time for the transfer is bounded with some dt_min and dt_max.
446
447	Args:
448	TL (float 1x3 ndarray):
449	TargetList class object
450	nA (integer):
451	Integer index of the current star of interest
452	nB (integer):
453	Integer index of the next star of interest
454	tA (astropy Time array):
455	Current absolute mission time in MJD
456
457	Returns:
458	tuple:
459	opt_slewTime (float):
460	Optimal slew time in days for starshade transfer to a
461	new line of sight
462	opt_dV (float):
463	Optimal total change in velocity in m/s for starshade
464	line of sight transfer
465
466	"""
467
468	def fuelUsage_objFun(dt, TL, nA, N, tA):	×
469	dV = self.calculate_dV(dt, TL, nA, N, tA)	×
470	return dV.value	×
471
472	def fuelUsage_constraints(dt, dt_min, dt_max):	×
473	return dt_max - dt, dt - dt_min	×
474
475	dt_guess = 20	×
476	dt_min = 1	×
477	dt_max = 45	×
478	Tol = 1e-5	×
479
480	t0 = [dt_guess]	×
481
482	res = optimize.minimize(	×
483	fuelUsage_objFun,
484	t0,
485	method="COBYLA",
486	args=(TL, nA, nB, tA),
487	constraints={
488	"type": "ineq",
489	"fun": fuelUsage_constraints,
490	"args": ([dt_min, dt_max]),
491	},
492	tol=Tol,
493	options={"disp": False},
494	)
495	opt_slewTime = res.x	×
496	opt_dV = res.fun	×
497
498	return opt_slewTime, opt_dV	×
499
500	def calculate_slewTimes(self, TL, old_sInd, sInds, sd, obsTimes, tmpCurrentTimeAbs):	1✔
501	"""Finds slew times and separation angles between target stars
502
503	This method determines the slew times of an occulter spacecraft needed
504	to transfer from one star's line of sight to all others in a given
505	target list.
506
507	Args:
508	TL (:ref:`TargetList`):
509	TargetList class object
510	old_sInd (int):
511	Integer index of the most recently observed star
512	sInds (~numpy.ndarray(int)):
513	Integer indices of the star of interest
514	sd (~astropy.units.Quantity):
515	Angular separation between stars in rad
516	obsTimes (~astropy.time.Time(~numpy.ndarray)):
517	Observation times for targets.
518	tmpCurrentTimeAbs (~astropy.time.Time(~numpy.ndarray)):
519	Current absolute mission time in MJD
520
521	Returns:
522	~astropy.units.Quantity:
523	Time to transfer to new star line of sight in units of days
524	"""
525
526	if old_sInd is None:	×
527	slewTimes = np.zeros(len(sInds)) * u.d	×
528	else:
529	obsTimeRangeNorm = (obsTimes - tmpCurrentTimeAbs).value	×
530	slewTimes = obsTimeRangeNorm[0, :] * u.d	×
531
532	return slewTimes	×
533
534	def log_occulterResults(self, DRM, slewTimes, sInd, sd, dV):	1✔
535	"""Updates the given DRM to include occulter values and results
536
537	Args:
538	DRM (dict):
539	Design Reference Mission, contains the results of one complete
540	observation (detection and characterization)
541	slewTimes (astropy Quantity):
542	Time to transfer to new star line of sight in units of days
543	sInd (integer):
544	Integer index of the star of interest
545	sd (astropy Quantity):
546	Angular separation between stars in rad
547	dV (astropy Quantity):
548	Delta-V used to transfer to new star line of sight in units of m/s
549
550	Returns:
551	dict:
552	Design Reference Mission dicitonary, contains the results of one
553	complete observation (detection and characterization)
554	"""
555
556	DRM["slew_time"] = slewTimes.to("day")	1✔
557	DRM["slew_angle"] = sd.to("deg")	1✔
558
559	dV = dV.to("m/s")	1✔
560	slew_mass_used = self.scMass * (	1✔
561	1 - np.exp(-dV.value / (self.slewIsp.value * const.g0.value))
562	)
563	DRM["slew_dV"] = dV	1✔
564	DRM["slew_mass_used"] = slew_mass_used.to("kg")	1✔
565	self.scMass = self.scMass - slew_mass_used	1✔
566	DRM["scMass"] = self.scMass.to("kg")	1✔
567	if self.twotanks:	1✔
568	self.slewMass = self.slewMass - slew_mass_used	1✔
569	DRM["slewMass"] = self.slewMass.to("kg")	1✔
570
571	return DRM	1✔

dsavransky / EXOSIMS / 11294204788

Source File Press 'n' to go to next uncovered line, 'b' for previous

Source File
Press 'n' to go to next uncovered line, 'b' for previous