1from pykep.trajopt import mga_1dsm, launchers
2from pykep.planet import jpl_lp
3from pykep import epoch_from_string
4
5import numpy as np
6from numpy.linalg import norm
7from math import log, acos, cos, sin, asin, exp
8from copy import deepcopy
9
10class _tandem_udp(mga_1dsm):
11    """
12    This class represents a rendezvous mission to Saturn modelled as an MGA-1DSM transfer. Mission parameters are
13    inspired to the TandEM mission. A launcher model (i.e. Atlas 501) is also used, so that the final mass delivered
14    at Saturn is the main objective of this optimization problem.
15
16    The problem draws inspiration from the work performed in April 2008 by the
17    European Space Agency working group on mission analysis on the mission named TandEM. TandEM is an interplanetary
18    mission aimed at reaching Titan and Enceladus (two moons of Saturn).
19
20    .. note::
21
22       The Titan and Enceladus Mission (TandEM), an ambitious scientific mission to study the Saturnian system
23       with particular emphasis on the moons Titan and Enceladus, was selected in October 2007 as a candidate mission
24       within the ESA Cosmic Vision plan. In February 2009, TandEM exited the Cosmic Vision programme when ESA
25       and NASA chose EJSM-Laplace as the L-class outer Solar System mission candidate.
26
27    .. note::
28
29       A significantly similar version of this problem was part of the no longer maintained GTOP database,
30       https://www.esa.int/gsp/ACT/projects/gtop/gtop.html. The exact definition is, though, different and results
31       cannot thus not be compared to those posted in GTOP.
32    """
33    def __init__(self, prob_id = 1, constrained = True):
34        """
35        The TandEM problem of the trajectory gym consists in 48 different instances varying in fly-by sequence and
36        the presence of a time constraint.
37
38        Args:
39            - prob_id (``int``): The problem id defines the fly-by sequence.
40            - constrained (``bool``): Activates the constraint on the time of flight
41              (fitness will thus return two numbers, the objective function and the inequality constraint violation)
42        """
43        # Redefining the planets as to change their safe radius
44        earth = jpl_lp('earth')
45        earth.safe_radius = 1.05
46        # We need the Earth eph in the fitness
47        venus = jpl_lp('venus')
48        venus.safe_radius = 1.05
49        mars = jpl_lp('mars')
50        mars.safe_radius = 1.05
51        jupiter = jpl_lp('jupiter')
52        jupiter.safe_radius = 1.7
53        saturn = jpl_lp('saturn')
54
55        # Defining the different sequences
56        seq_tandem = []
57        seq_tandem.append([earth, venus, venus, venus, saturn])
58        seq_tandem.append([earth, venus, venus, earth, saturn])
59        seq_tandem.append([earth, venus, venus, mars, saturn])
60        seq_tandem.append([earth, venus, venus, jupiter, saturn])
61
62        seq_tandem.append([earth, venus, earth, venus, saturn])
63        seq_tandem.append([earth, venus, earth, earth, saturn])
64        seq_tandem.append([earth, venus, earth, mars, saturn])
65        seq_tandem.append([earth, venus, earth, jupiter, saturn])
66
67        seq_tandem.append([earth, venus, mars, venus, saturn])
68        seq_tandem.append([earth, venus, mars, earth, saturn])
69        seq_tandem.append([earth, venus, mars, mars, saturn])
70        seq_tandem.append([earth, venus, mars, jupiter, saturn])
71
72        seq_tandem.append([earth, earth, venus, venus, saturn])
73        seq_tandem.append([earth, earth, venus, earth, saturn])
74        seq_tandem.append([earth, earth, venus, mars, saturn])
75        seq_tandem.append([earth, earth, venus, jupiter, saturn])
76
77        seq_tandem.append([earth, earth, earth, venus, saturn])
78        seq_tandem.append([earth, earth, earth, earth, saturn])
79        seq_tandem.append([earth, earth, earth, mars, saturn])
80        seq_tandem.append([earth, earth, earth, jupiter, saturn])
81
82        seq_tandem.append([earth, earth, mars, venus, saturn])
83        seq_tandem.append([earth, earth, mars, earth, saturn])
84        seq_tandem.append([earth, earth, mars, mars, saturn])
85        seq_tandem.append([earth, earth, mars, jupiter, saturn])
86
87        if prob_id > 24 or type(prob_id) != int or prob_id < 1:
88            raise ValueError("TandEM problem id must be an integer in [1, 24]")
89
90        super().__init__(
91            seq = seq_tandem[prob_id - 1],
92            t0 = [5475, 9132],
93            tof = [[20, 2500], [20, 2500], [20, 2500], [20, 2500]],
94            vinf = [2.5, 4.9],
95            add_vinf_dep = False,
96            add_vinf_arr = True,
97            tof_encoding = 'direct',
98            multi_objective = False,
99            orbit_insertion = True,
100            e_target = 0.98531407996358,
101            rp_target = 80330000,
102            eta_lb = 0.01,
103            eta_ub = 0.99,
104            rp_ub = 10
105        )
106
107        self.prob_id = prob_id
108        self.constrained = constrained
109
110    def fitness(self, x):
111        T, Vinfx, Vinfy, Vinfz = self._decode_times_and_vinf(x)
112        # We transform it (only the needed component) to an equatorial system rotating along x
113        # (this is an approximation, assuming vernal equinox is roughly x and the ecliptic plane is roughly xy)
114        earth_axis_inclination = 0.409072975
115        # This is different from the GTOP tanmEM problem, I think it was bugged there as the rotation was in the wrong direction.
116        Vinfz = - Vinfy * sin(earth_axis_inclination) + Vinfz * cos(earth_axis_inclination)
117        # And we find the vinf declination (in degrees)
118        sindelta = Vinfz / x[3]
119        declination = asin(sindelta) / np.pi * 180.
120        # We now have the initial mass of the spacecraft
121        m_initial = launchers.atlas501(x[3] / 1000., declination)
122        # And we can evaluate the final mass via Tsiolkowsky
123        Isp = 312.
124        g0 = 9.80665
125        DV = super().fitness(x)[0]
126        DV = DV + 165.  # losses for 3 swgbys + insertion
127        m_final = m_initial * exp(-DV / Isp / g0)
128        # Numerical guard for the exponential
129        if m_final == 0:
130            m_final = 1e-320
131        if self.constrained:
132            retval = [-log(m_final), sum(T) - 3652.5]
133        else:
134            retval = [-log(m_final), ]
135        return retval
136
137    def get_nic(self):
138        return int(self.constrained)
139
140    def get_name(self):
141        return "TandEM problem, id: " + str(self.prob_id)
142
143    def get_extra_info(self):
144        retval = "\t Sequence: " + \
145            [pl.name for pl in self._seq].__repr__() + "\n\t Constrained: " + \
146            str(self.constrained)
147        return retval
148
149    def pretty(self, x):
150        """
151        prob.plot(x)
152
153        - x: encoded trajectory
154
155        Prints human readable information on the trajectory represented by the decision vector x
156
157        Example::
158
159          print(prob.pretty(x))
160        """
161        super().pretty(x)
162        T, Vinfx, Vinfy, Vinfz = self._decode_times_and_vinf(x)
163        # We transform it (only the needed component) to an equatorial system rotating along x
164        # (this is an approximation, assuming vernal equinox is roughly x and the ecliptic plane is roughly xy)
165        earth_axis_inclination = 0.409072975
166        # This is different from the GTOP tanmEM problem, I think it was bugged there as the rotation was in the wrong direction.
167        Vinfz = - Vinfy * sin(earth_axis_inclination) + Vinfz * cos(earth_axis_inclination)
168        # And we find the vinf declination (in degrees)
169        sindelta = Vinfz / x[3]
170        declination = asin(sindelta) / np.pi * 180.
171        m_initial = launchers.atlas501(x[3] / 1000., declination)
172        # And we can evaluate the final mass via Tsiolkowsky
173        Isp = 312.
174        g0 = 9.80665
175        DV = super().fitness(x)[0]
176        DV = DV + 165.  # losses for 3 swgbys + insertion
177        m_final = m_initial * exp(-DV / Isp / g0)
178        print("\nInitial mass:", m_initial)
179        print("Final mass:", m_final)
180        print("Declination:", declination)
181
182    def __repr__(self):
183        return "TandEM (Trajectory Optimisation Gym P12, multiple instances)"
184
185# Problem P12: TandEM mission MGA1DSM, single objective, direct encoding, time constrained
186tandem = _tandem_udp
187