Orbiter-Forum PyKEP to TransX conversion
 Register Blogs Orbinauts List Social Groups FAQ Projects Mark Forums Read

 Math & Physics Mathematical and physical problems of space flight and astronomy.

 04-06-2016, 08:42 AM #1 Keithth G Orbinaut PyKEP to TransX conversion 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 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: * the velocity of Venus is: * the spacecraft's approach velocity vector (heliocentric coordinates again) is: * the spacecraft's departure velocity vector (heliocentric coordinates) is: (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 , and . Here 'f' (or 'forward') stands for prograde to avoid confusion with using 'p' for 'plane'. These vectors are constructed as follows: So, using the values given above: 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: In the above example, we calculate the spacecraft's departure velocity - relative to Venus - to be: Now, we can always write $v_{out}'$ as: where the values , and are just the components of our Venus escape plan. We calculate these quantities by taking the 'dot product' of with each of the TransX unit vectors as follows: If we apply this to our example, we calculate (in m/s) the following TransX escape plan: 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: ('hyperbolic escape velocity'), ('outward angle') and ('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: (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: 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 by Keithth G; 03-05-2017 at 11:54 PM.
 Thanked by:
 04-08-2016, 04:50 AM #2 Keithth G Orbinaut 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?lis...CRLz_T50-rtXr9 Last edited by Keithth G; 03-05-2017 at 11:54 PM.
 Thanked by:
 04-08-2016, 05:57 AM #3 Enjo Mostly harmless 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 by Enjo; 04-08-2016 at 08:15 AM.
 04-08-2016, 06:19 AM #4 Keithth G Orbinaut Quote: 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 by Keithth G; 03-05-2017 at 11:53 PM.
 02-25-2019, 09:58 PM #6 Ajaja Orbinaut PyKEP to TransX conversion 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 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 Attached Thumbnails     Last edited by Ajaja; 02-25-2019 at 10:07 PM.
 Thanked by:
 02-26-2019, 12:35 AM #7 MontBlanc2012 Orbinaut Thanks for the tool! (Aside: Quote: 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!)

 Posting Rules BB code is On Smilies are On [IMG] code is On HTML code is Off You may not post new threads You may not post replies You may not post attachments You may not edit your posts
 Forum Jump User Control Panel Private Messages Subscriptions Who's Online Search Forums Forums Home Orbiter-Forum.com     Announcements     Meets & Greets Orbiter Space Flight Simulator     Orbiter Web Forum         OFMM         Orbiter Forum Space Station         Simpit Forum     General Questions & Help     MFD Questions & Help     Hardware & Software Help     Tutorials & Challenges     Orbiter SDK     Orbiter Visualization Project     Orbiter Beta » Orbiter Project Orbiter Addons     OrbitHangar Addons & Comments     Addons     Addon Development     Addon Requests     Addon Support & Bugs         Addon Developer Forums             Project Apollo - NASSP     Orbiter Lua Scripting Far Side of the Moon     Spaceflight News     Math & Physics     Astronomy & the Night Sky     Backyard Rocketry     Brighton Lounge     International Forum

All times are GMT. The time now is 12:46 AM.