PyKEP to TransX conversion

Keithth G

New member
Joined
Nov 20, 2014
Messages
272
Reaction score
0
Points
0
I haven't had much time form Orbiter of late, but I've been musing on how to make use of PyKEP / PyGMO output within Orbiter. PyKEP / PyGMO is a free mission planning tool from the European Space Agency. However, it knows nothing about Orbiter and produces optimal trajectory plans in ways that are, most decidedly, Orbiter-unfriendly. To make practical use of PyKEP / PyGMO within Orbiter, then, one needs a way of converting PyKEP optimal trajectory plans into something that can be implemented in Orbiter using widely used MFDS - such as TransX and IMFD.

The purpose of this note, then, is to capture a few thoughts on these transformations. This is as much as an 'aide de memoire' for me as anything else since on more than one occasion, I have had recourse to rifling through my own posts to remind myself of useful equations and conversions.

Why convert from PyKEP-speak to Orbiter-speak?
PyKEP works by optimising a series of 'lambert arcs' that connect a series of planetary encounters. For example, for a simple Hohmann transfer from Earth to Jupiter, say, the two planetary encounters are: a departure from Earth; and an arrival at Jupiter. When PyKEP has completed its optimisation task, usually by minimising total [MATH]\Delta V[/MATH] requirements, it reports back the optimal times of the planetary encounters - e.g., the departure date from Earth; and the arrival date at Jupiter. However, PyKEP doesn't provide any useful information about how to set up the Earth 'escape plan'. Now, one can use TransX / IMFD and, by a process of trial-and-error, one can back-solve the encounter dates to work out what the 'escape plan' has to be to arrive at Jupiter on the planned arrival date. For more complicated trajectory plans involving a long sequence of slingshots and DSMs, using this trial-and-error technique to re-construct PyKEP's optimal trajectory solution in Orbiter can be challenging, if not impossible. Surely, there must be a better way than trial-and-error.

What output does PyKEP / PyGMO produce?
Let's start with the information that PyKEP does produce. First and foremost, it yields the times of the encounters with the planetary bodies. But PyKEP also has an ephemeris engine so it knows the position and velocity (the 'state vectors') of the planetary bodies. Now, PyKEP uses JPL's Low Precision ephemeris whereas Orbiter uses the VSOP87 ephemeris solution. Although these ephemerides are not exactly the same, for practical purposes the two ephemerides can be treated as equivalent.

Given that PyKEP knows 'where' the planetary bodies are, there is a function in PyKEP that allows one to re-construct the lambert arc connecting each encounter in the sequence. The lambert-arc is generally just a ballistic elliptical trajectory that connects the two planetary encounters and so it corresponds to the trajectory that a spacecraft needs to take in order to go from one encounter to the next in the specified amount of time. Part of the reconstruction of the Lambert arc is information about the velocity that the spacecraft must have when it leaves one planetary encounter; and the velocity that the spacecraft will have when it arrives at the next planetary encounter.

So, PyKEP provides the following information:

* the position of the spacecraft (and planetary body) at each encounter;
* the velocity of the planetary body at each encounter;
* the velocity of the spacecraft as it arrives at a planetary encounter; and
* the velocity of the spacecraft as it leaves a planetary encounter.

From this information we need to construct information that can be used in Orbiter's mission tools - e.g. TransX and IMFD. (In this note, I'm going to focus on TransX, since that tool is better suited to setting up slingshots than IMFD).

As an example from a PyKEP optimisation of a VEGA-style trajectory (Venus-Earth Gravity Assist), after extracting the relevant information from the lambert arcs, PyKEP tells me that at the Venus encounter (MJD 58002.0765):

* the position of Venus and the spacecraft (in a heliocentric co-ordinate system) is:
[MATH]r_v=(5070159035.023689, 107556835822.169, 1182577911.1819444)[/MATH]
* the velocity of Venus is:
[MATH] v_v=(-35101.3332441319, 1466.4390915935328, 2045.7339135134034)[/MATH]
* the spacecraft's approach velocity vector (heliocentric coordinates again) is:
[MATH] v_{in}=(-36817.31294777, -5469.38165681, -195.38174684)[/MATH]
* the spacecraft's departure velocity vector (heliocentric coordinates) is:
[MATH] v_{out} = (-41823.78381252, -1170.11804522, 63.30610391)[/MATH]

(In case you are wondering, the positions are the x-y-z coordinates in metres; and the velocities are the x-y-z coordinates in metres per second.)

From this information, we can accurately and reliably construct the 'escape plan' (if we are on Venus and leaving it); and the slingshot (if we are arriving from somewhere else and want to head of to another planetary destination). How do we do that? Well, let's begin by talking about TransX's coordinate system.

TransX's coordinate system - escape plans
Rather than use the x-y-z coordinate system of the underlying ephemeris, TransX works in another coordinate system consisting of 'prograde', 'outward' and 'plane'. These three quantities can be thought of as a new x-y-z coordinate system - in that the coordinates are all at right-angles to each other and form a dextral coordinate system in the same way that the underlying heliocentric x-y-z coordinate system does. To set up an escape plan, we need first to construct the unit vectors corresponding to TransX's 'prograde', 'outward' and 'plane' triplet of directions. Let's call these unit vectors [MATH]\hat{f}[/MATH], [MATH]\hat{o}[/MATH] and [MATH]\hat{p}[/MATH]. Here 'f' (or 'forward') stands for prograde to avoid confusion with using 'p' for 'plane'. These vectors are constructed as follows:

[MATH] \hat{f} = \frac{v_v}{\left\| v_v\right\| }[/MATH]
[MATH] \hat{p} = \frac{v_v \times r_v}{\left\| v_v \times r_v \right\|}[/MATH]
[MATH] \hat{o} = \hat{p} \times \hat{f}[/MATH]
So, using the values given above:

[MATH] \hat{f} = (-0.99743888, 0.04167031, 0.05813154) [/MATH]
[MATH] \hat{o} = ( 0.04239308, 0.99903761, 0.01125561) [/MATH]
[MATH] \hat{p} = (-0.05760657, 0.01369116, -0.99824548) [/MATH]
Having calculated these unit vectors, we can now calculate the information needed to form the escape plan. First, we need to work out the departure speed of the spacecraft relative to Venus:

[MATH] v_{out}' = v_{out} - v_v [/MATH]
In the above example, we calculate the spacecraft's departure velocity - relative to Venus - to be:

[MATH] v_{out}' = (-6722.45056839 -2636.55713681 -1982.4278096) [/MATH]
Now, we can always write $v_{out}'$ as:

[MATH] v_{out}' = v_f\,\hat{f} + v_o\,\hat{o} + v_p\,\hat{p} [/MATH]
where the values [MATH]v_f[/MATH], [MATH]v_o[/MATH] and [MATH]v_p[/MATH] are just the components of our Venus escape plan. We calculate these quantities by taking the 'dot product' of [MATH]v_{out}'[/MATH] with each of the TransX unit vectors as follows:

[MATH] v_f = v_{out}'.\hat{f} [/MATH]
[MATH] v_o = v_{out}'.\hat{o} [/MATH]
[MATH] v_p = v_{out}'.\hat{p} [/MATH]
If we apply this to our example, we calculate (in m/s) the following TransX escape plan:

[MATH] v_f = 6480.126 [/MATH]
[MATH] v_o = -2941.319 [/MATH]
[MATH] v_p = 2330.109 [/MATH]
And, to be sure, if starting from Venus one sets up this escape plan for MJD 58002.0765 we should find that TransX will predict a close encounter with Earth at the next scheduled encounter date of MJD 58572.5770.



The calculations here are quite general and can be applied to any PyKEP trajectory solution.


TransX's coordinate system - slingshots

Escape plans are generally used at the start of a sequence of encounters in TransX. In TransX, other planetary encounters (prior to arrival at the destination) are called 'slingshots' - aka 'gravity assists' and 'flybys'. Rather than using the cartesian coordinate system consisting of 'prograde', 'outward' and 'plane', TransX uses a spherical coordinate system for slingshots. Again, three components are needed: [MATH]v_\infty[/MATH] ('hyperbolic escape velocity'), [MATH]\theta[/MATH] ('outward angle') and [MATH]\phi[/MATH] ('inclination'). The 'hyperbolic escape velocity' is the speed of the spacecraft (relative to the planet) once it has escaped the planet's gravity well; the 'outward angle' is the angle between the prograde direction and the direction of departure from the planet (in the orbital plane of planet); and 'inclination' is the angle that the direction of departure makes with the planet's orbital plane. This sounds more complicated than it is in reality.

How do we calculate these quantities? Using the same notation as above:

[MATH] v_\infty = \left\| v_{out}' \right\|[/MATH]
[MATH]\theta = \tan ^{-1}\left(v_x,v_y\right)[/MATH]
[MATH]\phi = \tan^{-1}\left(\sqrt{v_x^2+v_y^2},v_z\right)[/MATH]
(Here, the arctan functions are of the 'atan2' variety - requiring both an 'x' and a 'y' input.) Normally in TransX, the slingshot 'inherits' the departure hyperbolic excess velocity from the arrival hyperbolic excess velocity. However, and in particular for powered flybys, one can specify a hyperbolic excess velocity that is faster or slower than the arrival hyperbolic excess velocity.

In our example, if we work though the numbers, we find that:

[MATH] v_\infty = 7488.177 \quad \text{m/s}[/MATH]
[MATH]\theta = -24.413 \quad \text{deg}[/MATH]
[MATH]\phi = 18.130 \quad \text{deg}[/MATH]
And, upon approach to Venus, if one inserts these numbers into TransX's slingshot plan (again for Venus arrival at MJD 58002.0765. The departure trajectory should again coincide with Earth arrival around MJD 58572.5770. However, this isn't as easy to set up *de novo* than the escape plan.

Setting up slingshots in this way is especially important if there is a DSM mid-way between one encounter and another. Without an ability to construct a planet departure trajectory in this fashion is becomes almost impossible (using TransX and IMFD) to ensure that the spacecraft is heading towards the DSM point.


PyKEP's full optimal VEGA trajectory for Earth departure MJD 57634.5001

For those interested, I've translated PyKEP's optimised VEGA-style trajectory plan for Earth departure MJD 57634.5001. This is based on a list of Jupiter trajectories provided in an earlier post by 'dgatsoulis'.

The plan is as follows:

Code:
Date of Earth departue:   2016-Sep-03 12:00:08.640002
Date of Venus encounter:  2017-Sep-06 01:50:05.290876
Date of Earth encounter:  2019-Mar-30 13:50:54.766147
Date of Jupiter arrival:  2021-Sep-24 07:39:52.717700

Transfer time from Earth to Venus:    367.58  days
Transfer time from Venus to Earth:    570.5  days
Transfer time from Earth to Jupiter:  908.74  days
Total mission duration:               1846.82  days


TransX escape plan - Earth escape
--------------------------------------
MJD:                 57634.5001 
Prograde:             -3461.680 m/s
Outward:               -291.288 m/s
Plane:                 -298.955 m/s
Hyp. excess velocity:  3486.754 m/s
Earth escape burn:     3743.016 m/s


Venus encounter
--------------------------------------
MJD:                 58002.0765 
Approach velocity:     7488.177 m/s
Departure velocity:    7488.177 m/s
Outward angle:          -24.413 deg
Inclination:             18.130 deg
Turning angle:           52.332 deg
Periapsis altitude:    1292.475 km 
dV needed:                0.000 m/s


Earth encounter
--------------------------------------
MJD:                 58572.5770 
Approach velocity:    12080.335 m/s
Departure velocity:   12425.696 m/s
Outward angle:           54.099 deg
Inclination:              4.666 deg
dV needed for flyby:    442.676 m/s


Jupiter arrival
--------------------------------------
MJD:                 59481.3194    
Hyp. excess velocity:  6193.941 m/s
Capture burn            322.023 m/s - (C3 = 0)


Total fuel cost:       4507.714 m/s

This plan is perfectly flyable in Orbiter now. Total delta-V cost of getting to Jupiter (and being captured by it) is 4,507.7 m/s (starting from a 300 km x 300 km orbit around Earth). It involves a ballistic encounter with Venus and a powered flyby of Earth. In comparison, an EGA-style Jupiter approach will typically cost around 5,300 m/s; and a Hohmann-stye transfer to Jupiter close to 7,000 m/s. Journey time is a little over five years.
 
Last edited:

Keithth G

New member
Joined
Nov 20, 2014
Messages
272
Reaction score
0
Points
0
To prove (to myself) that the PyKEP VEGA-style Earth to Jupiter trajectory plan shown above can indeed be successfully in Orbiter 2010, I ran the plan through Orbiter. Planned dV was, 4,507 m/s from a 300 x 300 Earth parking orbit to Jupiter capture. The outturn dV (inclusive of assorted mid-course corrections) was around 4,600 m/s.

I videos this flight and various instalments of it have been posted here:

https://www.youtube.com/playlist?list=PL0E82DPqx11rTj87g_gCRLz_T50-rtXr9
 
Last edited:

Enjo

Mostly harmless
Addon Developer
Tutorial Publisher
Donator
Joined
Nov 25, 2007
Messages
1,665
Reaction score
13
Points
38
Location
Germany
Website
www.enderspace.de
Preferred Pronouns
Can't you smell my T levels?
Very nice. Definitely deserves to be archived in the Tutorials section... under the Advanced tab :)
I'm also glad that you're able to present the conjunction between LaunchMFD's azimuthal guidance with specific vessels' pitch guidance. I'm not sure how to explain why the correction would mess up more than help. Possibly because of varying thrust along the way to the parking orbit. But since you launch at the exact time point that Launch MFD suggests, the correction isn't that much needed.

I've seen only parts 1 and 2, so can't comment the whole series, but one thing that comes into my mind is that I'd strongly recommend increasing the MFD refresh rate in the Parameters tab. 0.2 s would be good enough.
 
Last edited:

Keithth G

New member
Joined
Nov 20, 2014
Messages
272
Reaction score
0
Points
0
one thing that comes into my mind is that I'd strongly recommend increasing the MFD refresh rate in the Parameters tab. 0.2 s would be good enough.

Ah, yes. Good point!
 
Last edited:

boogabooga

Bug Crusher
Joined
Apr 16, 2011
Messages
2,999
Reaction score
1
Points
0
Keithth G, please don't delete your youtube videos.
 

Ajaja

Active member
Joined
Apr 20, 2008
Messages
226
Reaction score
93
Points
28
I've rewritten the code of this amazing Keithth G's tool for new Pykep/Pygmo v 2.x and Python v3.6 with some improvements such as multi-revolution solutions, checking of altitude, SPICE and GMAT suppport, trajectory plot, sequenced input.

It's based on the Keithth G's code (from https://www.orbiter-forum.com/showthread.php?t=36989 mostly), so I put it here without creating new thread.

The main solver flyby.py :
Code:
from numpy import *
from math import *
from pykep import epoch, lambert_problem, fb_vel, DAY2SEC, AU, MU_SUN, RAD2DEG
from pykep.planet import keplerian, jpl_lp, spice

class flyby:

    jpl_planet = {  'mercury'   : jpl_lp('mercury'),
                    'venus'     : jpl_lp('venus'),
                    'earth'     : jpl_lp('earth'),
                    'mars'      : jpl_lp('mars'),
                    'jupiter'   : jpl_lp('jupiter'),
                    'saturn'    : jpl_lp('saturn'),
                    'uranus'    : jpl_lp('uranus'),
                    'neptune'   : jpl_lp('neptune'),
                    'pluto'     : jpl_lp('pluto')
                 }


    mu_sun_spice = 132712440017.99e9
    spice_planet = {'mercury'   : spice('MERCURY', 'SUN', 'ECLIPJ2000', 'NONE', mu_sun_spice, 22032.080486418e9, 2439.7e3, 2439.7e3 * 1.1),
                    'venus'     : spice('VENUS', 'SUN', 'ECLIPJ2000', 'NONE', mu_sun_spice, 324858.59882646e9, 6051.9e3, 6051.9e3 * 1.1),
                    'earth'     : spice('EARTH', 'SUN', 'ECLIPJ2000', 'NONE', mu_sun_spice, 398600.4415e9, 6378.1363e3, 6378.1363e3 * 1.1),
                    'mars'      : spice('MARS', 'SUN', 'ECLIPJ2000', 'NONE', mu_sun_spice, 42828.314258067e9, 3397.0e3, 3397.0e3 * 1.1),
                    'jupiter'   : spice('JUPITER', 'SUN', 'ECLIPJ2000', 'NONE', mu_sun_spice, 126712767.8578e9, 71492.0e3, 71492.0e3 * 1.1),
                    'saturn'    : spice('SATURN', 'SUN', 'ECLIPJ2000', 'NONE', mu_sun_spice, 37940626.061137e9, 60268.0e3, 60268.0e3 * 1.1),
                    'uranus'    : spice('URANUS', 'SUN', 'ECLIPJ2000', 'NONE', mu_sun_spice, 5794549.0070719e9, 25559.0e3, 25559.0e3 * 1.1),
                    'neptune'   : spice('NEPTUNE', 'SUN', 'ECLIPJ2000', 'NONE', mu_sun_spice, 6836534.0638793e9, 25269.0e3, 25269.0e3 * 1.1),
                    'pluto'     : spice('PLUTO', 'SUN', 'ECLIPJ2000', 'NONE', mu_sun_spice, 981.600887707e9, 1162.0e3, 1162.0e3 * 1.1)
        }


    def __init__(self, flight_plan, travel_days, windows, ignore_last=False, orbit_alt=300000, days=1, spice=False, multi_revs=5):

        self.travel_days = travel_days.copy()
        self.windows = windows.copy()
        self.orbit_alt = orbit_alt
        self.ignore_last = ignore_last
        self.dim = len(flight_plan)
        self.flight_plan = flight_plan.copy()
        self.days = days
        self.spice = spice
        self.multi_revs = multi_revs

        self.planets = []
        if self.spice :
            self.mu_sun = self.mu_sun_spice
            for i in flight_plan:
                self.planets.append(self.spice_planet[i])
        else:
            self.mu_sun = MU_SUN
            for i in flight_plan:
                self.planets.append(self.jpl_planet[i])


        self.x = [0] * self.dim
        self.t = [0] * self.dim
        self.r = [0] * self.dim
        self.v = [0] * self.dim
        self.vo = [0] * self.dim
        self.vi = [0] * self.dim
        self.f = [0] * self.dim
        self.f_all = 0
        self.li_sol = []
        self.l = []

    def get_bounds(self):
        return ([-1 * x for x in self.windows], self.windows)
    
    def get_name(self):
        return "Flyby"

    def fitness(self, x):
        self.x = x
        # calculate the times
        self.t[0] = self.travel_days[0] + self.x[0]
        for i in range(1, self.dim):
            self.t[i] = self.days * self.t[i - 1] + self.travel_days[i] + self.x[i]

        # calculate the state vectors of planets
        for i in range(self.dim):
            self.r[i], self.v[i] = self.planets[i].eph(epoch(self.t[i], "mjd"))

        # calculate the solutions of the two Lambert transfers
        self.l = []
        n_sols = []
        for i in range(self.dim - 1):
            self.l.append(lambert_problem(self.r[i], self.r[i + 1], (self.t[i + 1] - self.t[i]) * DAY2SEC, self.mu_sun, False, self.multi_revs))
            n_sols.append(self.l[i].get_Nmax() * 2 + 1)
            
        # perform the dV calculations
        mu0 = self.planets[0].mu_self
        rad0 = self.planets[0].radius + self.orbit_alt
        mu1 = self.planets[-1].mu_self
        rad1 = self.planets[-1].radius + self.orbit_alt
        
        k = 1
        for i in range(self.dim - 1):
            k = k * n_sols[i]
        
        vot = [0] * self.dim
        vit = [0] * self.dim
        ft = [0] * self.dim
        self.f_all = 1.0e10
        
        for kk in range(k):
            d = kk
            li = []
            for j in range(self.dim - 1):
                d , b = divmod(d , n_sols[j])
                li.append(b)
                
            vot[0] = array(self.l[0].get_v1()[li[0]]) - self.v[0]
            ft[0] = sqrt(dot(vot[0], vot[0]) + 2 * mu0 / rad0) - sqrt(1 * mu0 / rad0)

            if ft[0] > self.f_all:
                continue

            for i in range(1,self.dim - 1):
                vit[i] = array(self.l[i - 1].get_v2()[li[i - 1]]) - self.v[i]
                vot[i] = array(self.l[i].get_v1()[li[i]]) - self.v[i]
                ft[i] = fb_vel(vit[i], vot[i], self.planets[i])
            
            vit[-1] = array(self.l[-1].get_v2()[li[-1]]) - self.v[-1]
            ft[-1] = sqrt(dot(vit[-1], vit[-1]) + 2 * mu1 / rad1) - sqrt(2 * mu1 / rad1)
            
            ft_all = sum(ft)
            if self.ignore_last :
                ft_all = ft_all - ft[-1]

            if ft_all < self.f_all:
                self.f_all = ft_all
                self.vi = vit.copy()
                self.vo = vot.copy()
                self.f = ft.copy()
                self.li_sol = li.copy()


        # check and set cost of negative altitude (using safe_radius)
        res = self.f_all
        for i in range(1, self.dim - 1):
            ta = acos(dot(self.vi[i], self.vo[i]) / sqrt(dot(self.vi[i], self.vi[i])) / sqrt(dot(self.vo[i], self.vo[i])))
            alt = (self.planets[i].mu_self / dot(self.vi[i], self.vi[i]) * (1 / sin(ta / 2) - 1)) - self.planets[i].safe_radius
            if alt < 0:
                res = res - alt

        # return the total fuel cost
        return [res]

    def get_extra_info(self):
        return "\tDimensions: " + str(self.dim)

    def plot_trajectory(self):
        from matplotlib.pyplot import figure, show
        from pykep.orbit_plots import plot_planet, plot_lambert

        fig = figure()
        ax = fig.gca(projection='3d', aspect='equal')
        ax.scatter([0], [0], [0], color='y')

        colors = {'mercury': '#7B7869', 
                  'venus'  : '#BB91A1',
                  'earth'  : '#0000FF',
                  'mars'   : '#E27B58',
                  'jupiter': '#C88B3A',
                  'saturn' : '#A49B72',
                  'uranus' : '#65868B',
                  'neptune': '#6081FF',
                  'pluto'  : '#333333'}

        d_max = 0
        for i in range(self.dim):
            r, v = self.planets[i].eph(epoch(self.t[i], "mjd"))
            d = dot(r,r)
            if d > d_max:
                d_max = d

            p = keplerian(epoch(self.t[i], "mjd"),
                          self.planets[i].osculating_elements(epoch(self.t[i], "mjd")),
                          self.planets[i].mu_central_body,
                          self.planets[i].mu_self,
                          self.planets[i].radius,
                          self.planets[i].safe_radius,
                          self.flight_plan[i])
            plot_planet(p, epoch(self.t[i], "mjd"),  units=AU, ax=ax, color=colors[self.flight_plan[i]])
            if i != self.dim - 1:
                plot_lambert(self.l[i], sol=self.li_sol[i], units=AU, ax=ax, color='c')

    
        d_max = 1.2 * sqrt(d_max) / AU
        ax.set_xlim(-d_max, d_max)
        ax.set_ylim(-d_max, d_max)
        ax.set_zlim(-d_max, d_max)

        show()

    def print_transx(self):

        print("Date of %8s departue : " % self.flight_plan[0], epoch(self.t[0], "mjd"))
        for i in range(1,self.dim - 1):
            print("Date of %8s encounter: " % self.flight_plan[i], epoch(self.t[i], "mjd"))
        print("Date of %8s arrival  : " % self.flight_plan[-1], epoch(self.t[-1], "mjd"))
        print("")

        for i in range(self.dim - 1):
            print("Transfer time from %8s to %8s:" % (self.flight_plan[i], self.flight_plan[i + 1]), round(self.t[i + 1] - self.t[i], 4), " days")
        
        print("Total mission duration:                 ", round(self.t[-1] - self.t[0], 4), " days")
        print("")
        print("")

        fward = [0] * self.dim
        plane = [0] * self.dim
        oward = [0] * self.dim
        for i in range(self.dim):
            fward[i] = self.v[i] / linalg.norm(self.v[i])
            plane[i] = cross(self.v[i], self.r[i])
            plane[i] = plane[i] / linalg.norm(plane[i])
            oward[i] = cross(plane[i], fward[i])

        print("TransX escape plan -  %8s escape" % self.flight_plan[0])
        print("--------------------------------------")
        print("MJD:                  %10.4f " % round(epoch(self.t[0], "mjd").mjd, 4))
        print("Prograde:             %10.4f m/s" % round(dot(fward[0], self.vo[0]), 4))
        print("Outward:              %10.4f m/s" % round(dot(oward[0], self.vo[0]), 4))
        print("Plane:                %10.4f m/s" % round(dot(plane[0], self.vo[0]), 4))
        print("Hyp. excess velocity: %10.4f m/s" % round(linalg.norm(self.vo[0]), 4))
        print("Earth escape burn:    %10.4f m/s" % round(self.f[0], 4))

        c3 = dot(self.vo[0],self.vo[0]) / 1000000
        dha = atan2(self.vo[0][2],sqrt(self.vo[0][0] * self.vo[0][0] + self.vo[0][1] * self.vo[0][1])) * RAD2DEG
        rha = atan2(self.vo[0][1],self.vo[0][0]) * RAD2DEG

        print("GMAT MJD:             ", epoch(self.t[0], "mjd").mjd - 29999.5)
        print("GMAT OutgoingC3:      ", c3)
        print("GMAT OutgoingRHA:     ", rha)
        print("GMAT OutgoingDHA:     ", dha)
        print("")

        for i in range(1, self.dim - 1):
            vx = dot(fward[i], self.vo[i])
            vy = dot(oward[i], self.vo[i])
            vz = dot(plane[i], self.vo[i])
            mu = self.planets[i].mu_self
            rad = self.planets[i].radius
            print("%8s encounter" % self.flight_plan[i])
            print("--------------------------------------")
            print("MJD:                 %10.4f " % round(epoch(self.t[i], "mjd").mjd, 4))
            print("Solution number:     %10d " % (1 + self.li_sol[i - 1]))
            print("Approach velocity:   %10.4f m/s" % round(linalg.norm(self.vi[i]), 4))
            print("Departure velocity:  %10.4f m/s" % round(linalg.norm(self.vo[i]), 4))
            print("Outward angle:       %10.4f deg" % round(atan2(vy, vx) * RAD2DEG, 4))
            print("Inclination:         %10.4f deg" % round(atan2(vz, sqrt(vx * vx + vy * vy)) * RAD2DEG, 4))

            a = - mu / dot(self.vi[i], self.vi[i])
            ta = acos(dot(self.vi[i], self.vo[i]) / linalg.norm(self.vi[i]) / linalg.norm(self.vo[i]))
            e = 1 / sin(ta / 2)
            rp = a * (1 - e) 
            alt = (rp - rad) / 1000
           
            print("Turning angle:       %10.4f deg" % round(ta * RAD2DEG, 4)) 
            print("Periapsis altitude:  %10.4f km " % round(alt, 4))
            print("dV needed:           %10.4f m/s" % round(self.f[i], 4))
            print("GMAT MJD:             ", epoch(self.t[i], "mjd").mjd - 29999.5)
            print("GMAT RadPer:          ", rp / 1000)

            c3 = dot(self.vi[i],self.vi[i]) / 1000000
            dha = atan2(-self.vi[i][2],sqrt(self.vi[i][0] * self.vi[i][0] + self.vi[i][1] * self.vi[i][1])) * RAD2DEG
            rha = atan2(-self.vi[i][1],-self.vi[i][0]) * RAD2DEG
            if rha < 0:
                rha = 360 + rha

            print("GMAT IncomingC3:      ", c3)
            print("GMAT IncomingRHA:     ", rha)
            print("GMAT IncomingDHA:     ", dha)

            e = cross([0,0,1],-self.vi[i])
            e = e / linalg.norm(e)
            n = cross(-self.vi[i],e)
            n = n / linalg.norm(n)
            h = cross(self.vi[i],self.vo[i])
            b = cross(h,-self.vi[i])
            b = b / linalg.norm(b)
            sinb = dot(b,e)
            cosb = dot(b,-n)
            bazi = atan2(sinb,cosb) * RAD2DEG
            if bazi < 0:
                bazi = bazi + 360
            print("GMAT IncomingBVAZI:   ", bazi) 

            c3 = dot(self.vo[i],self.vo[i]) / 1000000
            dha = atan2(self.vo[i][2],sqrt(self.vo[i][0] * self.vo[i][0] + self.vo[i][1] * self.vo[i][1])) * RAD2DEG
            rha = atan2(self.vo[i][1],self.vo[i][0]) * RAD2DEG
            if rha < 0:
                rha = 360 + rha
           
            print("GMAT OutgoingC3:      ", c3)
            print("GMAT OutgoingRHA:     ", rha)
            print("GMAT OutgoingDHA:     ", dha)


            e = cross([0,0,1],self.vo[i])
            e = e / linalg.norm(e)
            n = cross(self.vo[i],e)
            n = n / linalg.norm(n)
            h = cross(self.vi[i],self.vo[i])
            b = cross(h,self.vo[i])
            b = b / linalg.norm(b)
            sinb = dot(b,e)
            cosb = dot(b,-n)
            bazi = atan2(sinb,cosb) * RAD2DEG
            if bazi < 0:
                bazi = bazi + 360
            print("GMAT OutgoingBVAZI:   ", bazi) 

            print("")


        print("%8s arrival" % self.flight_plan[-1])
        print("--------------------------------------")
        print("MJD:                  %10.4f    " % round(epoch(self.t[-1], "mjd").mjd, 4))
        print("Solution number:      %10d " % (1 + self.li_sol[-1] + 1))
        print("Hyp. excess velocity: %10.4f m/s" % round(sqrt(dot(self.vi[-1], self.vi[-1])), 4))
        print("Orbit insertion burn  %10.4f m/s - C3 = 0" % round(self.f[-1], 4))

        c3 = dot(self.vi[-1],self.vi[-1]) / 1000000
        dha = atan2(self.vi[-1][2],sqrt(self.vi[-1][0] * self.vi[-1][0] + self.vi[-1][1] * self.vi[-1][1])) * RAD2DEG
        rha = atan2(self.vi[-1][1],self.vi[-1][0]) * RAD2DEG

        print("GMAT MJD:             ", epoch(self.t[-1], "mjd").mjd - 29999.5)
        print("GMAT IncomingC3:      ", c3)
        print("GMAT IncomingRHA:     ", rha)
        print("GMAT IncomingDHA:     ", dha)
        print("")


        print("--------------------------------------")
        print("Total fuel cost:     %10.4f m/s" % round(self.f_all, 4))


an example how to use it:

Code:
def run():
    from pygmo import archipelago, problem, sade
    from flyby import flyby

    p = flyby(['earth','venus', 'earth', 'mars', 'earth', 'jupiter'],
              [60094.0, 60240.0, 60555.0, 60717.0, 61370.0, 62432.0], 
              [15, 15, 15, 15, 15, 15], ignore_last=False, days=0)

    archi = archipelago(n=4, algo=sade(gen = 100), prob=problem(p), pop_size=10)
    archi.evolve(100)
    archi.wait()

    tmp = [isl for isl in archi]
    tmp.sort(key=lambda x: x.get_population().champion_f)

    p.fitness(tmp[0].get_population().champion_x)
    p.print_transx()
    p.plot_trajectory()


if __name__ == "__main__":
    run()


and result:
Code:
Date of    earth departue :  2023-May-27 17:50:04.843458
Date of    venus encounter:  2023-Oct-21 19:51:18.197115
Date of    earth encounter:  2024-Sep-01 09:15:58.847174
Date of     mars encounter:  2025-Feb-10 21:35:16.139662
Date of    earth encounter:  2026-Nov-25 19:46:45.000977
Date of  jupiter arrival  :  2029-Nov-06 23:59:58.709198

Transfer time from    earth to    venus: 147.0842  days
Transfer time from    venus to    earth: 315.5588  days
Transfer time from    earth to     mars: 162.5134  days
Transfer time from     mars to    earth: 652.9246  days
Transfer time from    earth to  jupiter: 1077.1759  days
Total mission duration:                  2355.2569  days


TransX escape plan -     earth escape
--------------------------------------
MJD:                  60091.7431 
Prograde:             -2499.4183 m/s
Outward:               -800.4121 m/s
Plane:                -1893.8489 m/s
Hyp. excess velocity:  3236.4201 m/s
Earth escape burn:     3669.4054 m/s
GMAT MJD:              30092.2431116141
GMAT OutgoingC3:       10.474415101344679
GMAT OutgoingRHA:      137.52804616165528
GMAT OutgoingDHA:      35.81279004742669

   venus encounter
--------------------------------------
MJD:                 60238.8273 
Solution number:              1 
Approach velocity:    5586.3478 m/s
Departure velocity:   5586.3478 m/s
Outward angle:         -20.3629 deg
Inclination:            24.8403 deg
Turning angle:          47.9394 deg
Periapsis altitude:   9161.9881 km 
dV needed:               0.0000 m/s
GMAT MJD:              30239.327293948096
GMAT RadPer:           15213.98811419788
GMAT IncomingC3:       31.2072819829571
GMAT IncomingRHA:      317.54996341944445
GMAT IncomingDHA:      57.844582760539
GMAT IncomingBVAZI:    59.63517420647543
GMAT OutgoingC3:       31.207281983365895
GMAT OutgoingRHA:      181.08245369321392
GMAT OutgoingDHA:      -21.556120580924805
GMAT OutgoingBVAZI:    150.41306260145242

   earth encounter
--------------------------------------
MJD:                 60554.3861 
Solution number:              1 
Approach velocity:    8725.6939 m/s
Departure velocity:   8725.6939 m/s
Outward angle:         -52.4279 deg
Inclination:           -13.8499 deg
Turning angle:          45.3102 deg
Periapsis altitude:   1978.3521 km 
dV needed:               0.0000 m/s
GMAT MJD:              30554.886097768212
GMAT RadPer:           8356.352126270855
GMAT IncomingC3:       76.13773321840134
GMAT IncomingRHA:      345.49247842565273
GMAT IncomingDHA:      1.1249877490661295
GMAT IncomingBVAZI:    290.86336800854485
GMAT OutgoingC3:       76.13773321926722
GMAT OutgoingRHA:      122.32012280403704
GMAT OutgoingDHA:      13.847161531688322
GMAT OutgoingBVAZI:    254.1996985918471

    mars encounter
--------------------------------------
MJD:                 60716.8995 
Solution number:              1 
Approach velocity:   10035.1259 m/s
Departure velocity:  10035.1259 m/s
Outward angle:          98.8373 deg
Inclination:            -0.9846 deg
Turning angle:           9.8285 deg
Periapsis altitude:   1142.2797 km 
dV needed:               0.0000 m/s
GMAT MJD:              30717.399492357203
GMAT RadPer:           4539.2796868534
GMAT IncomingC3:       100.70375201811501
GMAT IncomingRHA:      301.6450576071275
GMAT IncomingDHA:      5.296782669850119
GMAT IncomingBVAZI:    324.0995173504984
GMAT OutgoingC3:       100.70375201711299
GMAT OutgoingRHA:      115.89414740991516
GMAT OutgoingDHA:      2.6779688839493043
GMAT OutgoingBVAZI:    215.76863635275086

   earth encounter
--------------------------------------
MJD:                 61369.8241 
Solution number:              1 
Approach velocity:   11456.1786 m/s
Departure velocity:  11456.1786 m/s
Outward angle:          45.5618 deg
Inclination:            -9.0839 deg
Turning angle:          26.9791 deg
Periapsis altitude:   3604.6660 km 
dV needed:               0.0000 m/s
GMAT MJD:              31370.324131955756
GMAT RadPer:           9982.665965148264
GMAT IncomingC3:       131.24402832682992
GMAT IncomingRHA:      261.1904853458861
GMAT IncomingDHA:      -5.915312967108089
GMAT IncomingBVAZI:    81.59217261929847
GMAT OutgoingC3:       131.2440283277017
GMAT OutgoingRHA:      108.22232768208467
GMAT OutgoingDHA:      9.080597356809768
GMAT OutgoingBVAZI:    94.81330768532506

 jupiter arrival
--------------------------------------
MJD:                  62447.0000    
Solution number:               2 
Hyp. excess velocity:  5572.4852 m/s
Orbit insertion burn    260.7795 m/s - C3 = 0
GMAT MJD:              32447.49998506016
GMAT IncomingC3:       31.05259174483021
GMAT IncomingRHA:      123.57110111953926
GMAT IncomingDHA:      -1.796759856968099

--------------------------------------
Total fuel cost:      3930.1850 m/s

attachment.php


An example of a multi-revolution trajectory and usage of SPICE kernels:
Code:
from pykep.util import load_spice_kernel
load_spice_kernel('E:/Utils/GMAT/data/planetary_ephem/spk/DE424AllPlanets.bsp')

def run():
    from pygmo import archipelago, problem, sade
    from flyby import flyby
    p = flyby(['earth','venus','earth','earth', 'jupiter'],
              [62445, 146, 318, 827, 754], 
              [1, 10, 20, 30, 40 ] ,  days = 1, spice = True, multi_revs=3)

    archi = archipelago(n=4, algo=sade(gen = 1000), prob=problem(p), pop_size=30)
    archi.evolve(100)
    archi.wait()

    tmp = [isl for isl in archi]
    tmp.sort(key=lambda x: x.get_population().champion_f)

    p.fitness(tmp[0].get_population().champion_x)
    p.print_transx()
    p.plot_trajectory()


if __name__ == "__main__":
    run()

Code:
Date of    earth departue :  2029-Nov-04 00:00:00.000041
Date of    venus encounter:  2030-Mar-30 09:30:39.106430
Date of    earth encounter:  2031-Feb-13 05:05:36.473347
Date of    earth encounter:  2033-May-20 14:31:29.815307
Date of  jupiter arrival  :  2035-Jul-23 14:31:29.815280

Transfer time from    earth to    venus: 146.3963  days
Transfer time from    venus to    earth: 319.8159  days
Transfer time from    earth to    earth: 827.393  days
Transfer time from    earth to  jupiter: 794.0  days
Total mission duration:                  2087.6052  days


TransX escape plan -     earth escape
--------------------------------------
MJD:                  62444.0000 
Prograde:             -2732.5248 m/s
Outward:              -1020.3482 m/s
Plane:                 1494.3723 m/s
Hyp. excess velocity:  3277.3390 m/s
Earth escape burn:     3681.0665 m/s
GMAT MJD:              32444.50000000048
GMAT OutgoingC3:       10.740950997778658
GMAT OutgoingRHA:      -68.217681978829
GMAT OutgoingDHA:      -27.123393131390685

   venus encounter
--------------------------------------
MJD:                 62590.3963 
Solution number:              1 
Approach velocity:    5415.4688 m/s
Departure velocity:   5415.4688 m/s
Outward angle:          -7.7117 deg
Inclination:           -23.8836 deg
Turning angle:          51.9385 deg
Periapsis altitude:   8167.5004 km 
dV needed:               0.0000 m/s
GMAT MJD:              32590.896285954055
GMAT RadPer:           14219.400366937063
GMAT IncomingC3:       29.327302369028473
GMAT IncomingRHA:      101.60425523038312
GMAT IncomingDHA:      -47.16187793154653
GMAT IncomingBVAZI:    100.90619977938115
GMAT OutgoingC3:       29.32730236902684
GMAT OutgoingRHA:      337.2525014900944
GMAT OutgoingDHA:      20.53524717915198
GMAT OutgoingBVAZI:    45.47561229116398

   earth encounter
--------------------------------------
MJD:                 62910.2122 
Solution number:              1 
Approach velocity:    9556.9446 m/s
Departure velocity:   9556.9446 m/s
Outward angle:         -64.4535 deg
Inclination:            -0.0074 deg
Turning angle:          31.4670 deg
Periapsis altitude:   5351.9058 km 
dV needed:               0.0000 m/s
GMAT MJD:              32910.71222770077
GMAT RadPer:           11730.042065142294
GMAT IncomingC3:       91.33518934304068
GMAT IncomingRHA:      149.00208869332303
GMAT IncomingDHA:      -3.040602403492167
GMAT IncomingBVAZI:    265.0398889254878
GMAT OutgoingC3:       91.33518934303952
GMAT OutgoingRHA:      297.6663412312752
GMAT OutgoingDHA:      0.009898972591192219
GMAT OutgoingBVAZI:    275.81590537927144

   earth encounter
--------------------------------------
MJD:                 63737.6052 
Solution number:              3 
Approach velocity:    9575.5225 m/s
Departure velocity:   9575.5225 m/s
Outward angle:          24.4237 deg
Inclination:            11.5481 deg
Turning angle:          40.8381 deg
Periapsis altitude:   1735.0483 km 
dV needed:               0.0000 m/s
GMAT MJD:              33738.10520619569
GMAT RadPer:           8113.184568383894
GMAT IncomingC3:       91.69063119273551
GMAT IncomingRHA:      84.86678298539435
GMAT IncomingDHA:      -0.0005654178930187623
GMAT IncomingBVAZI:    107.82072562225858
GMAT OutgoingC3:       91.69063119273544
GMAT OutgoingRHA:      304.31647297041286
GMAT OutgoingDHA:      -11.543920197638524
GMAT OutgoingBVAZI:    76.33028582219168

 jupiter arrival
--------------------------------------
MJD:                  64531.6052    
Solution number:               2 
Hyp. excess velocity:  6236.3982 m/s
Orbit insertion burn    326.4068 m/s - C3 = 0
GMAT MJD:              34532.10520619537
GMAT IncomingC3:       38.892662492158415
GMAT IncomingRHA:      -36.9371557553744
GMAT IncomingDHA:      1.7124666576344567

--------------------------------------
Total fuel cost:      4007.4733 m/s
attachment.php
 

Attachments

  • Figure_1.png
    Figure_1.png
    144 KB · Views: 59
  • Figure_2.png
    Figure_2.png
    110 KB · Views: 59
Last edited:

MontBlanc2012

Active member
Joined
Aug 13, 2017
Messages
138
Reaction score
26
Points
28
Thanks for the tool!


(Aside:
from pykep.util import load_spice_kernel
load_spice_kernel('E:/Utils/GMAT/data/planetary_ephem/spk/DE424AllPlanets.bsp')

I had wondered how one loaded the spice data into pykep. Now I know!)
 
Top