christopherrabotin / smd

Space Mission Design - A SPICE-enhanced continuous thrust interplanetary mission propagator and vizualizer

Space Mission Design (smd)

Space Mission Design allows one to perform an initial space mission design, around a given celestial body or between celestial bodies.

This package was written to support my thesis and my astrodynamics courses (ASEN 6008 Space Mission Design / Interplanetary Mission Design) at the University of Colorado Boulder.

Features

Note: this list may not be up to date with the latest developments.

• Propagation of an orbit around a celestial body
• Direct closed-loop optimization of continuous thrust via Naasz and Ruggiero control laws.
• VSOP87 support via the amazing https://github.com/soniakeys/meeus
• Patched conics for interplanetary missions
• Stream orbital elements as CSV for live visualization of how they change
• Export as a set of NASA Cosmographia files (cf. http://cosmoguide.org/) for really cool visualization of the overall mission
• Export mission state as CSV (cf. the examples/statOD/main.go)

Usage

If running smd and planning on changing reference frames (e.g. when doing patched conics) to attempting to include third body dynamics, you will need to define the SMD_CONFIG environment variable. This must define whether using VSOP87 or SPICE for frame transformations. An example of such a file is found in conf.toml. Important: this configuration file must be called conf.toml (but it can be placed in any directory). Note: the availability of this file will only occur in the function which gets the heliocentric orbit of a given planet. So definitely make sure this is configured before running a long simulation or it will crash when you're looking away.

smd's Issues

IMD: Homework 5

VEEGA transfer

Argument of perigee correction plane based on efficiency

Petropoulos (via Hatten) proposes two (working) laws to correct the argument of perigee: the one which is in plane does not alter the other orbital elements, but the one which is out of plane will. Both have different efficiency factors, so what would be quite interesting is switching between those laws based on which one is more efficient at a given time.

IMD: Homework 3

StatOD : HW2

Range and range rate measurements for the three following orbits:

• LEO: a = 7000 km, e = 0.001, i = 30 degrees, Ω = 80 degrees, ω = 40 degrees, and an initial true anomaly of ν = 0 degrees.
• GEO: a = 42163 km, e = 1 × 10−5 , i = 1 degrees, Ω = 70 degrees, ω = 0 degrees, and an initial true anomaly of ν = 180 degrees.
• Hyperbolic flyby: ECI initial conditions are given by r0 = [−7737.559071593195 −43881.87809094457 0.0]km, v0 = [3.347424567061589 3.828541915617483 0.0] km/s.

And the three following ground stations (lat/long (degrees)):

• Station 1 -35.398333, 148.981944
• Station 2 40.427222, 355.749444
• Station 3 35.247164, 243.205

Fix SOI edge case

When propagating (after #27), the SOI switch fails (h=NaN => \dot{a}=NaN). This must be fixed.

IMD: Lab 4 - Contour plots

Pseudocode:

Departure_Dates = vector of departure dates // In days, must be plottable
Arrival_Dates = vector of arrival dates // In days, must be plottable
Loop Departure_Dates
	Loop Arrival_Dates
		[R_Earth, V_Earth] = Earth_RV(JD_Depart)
		[R_Mars,  V_Mars]  = Mars_RV(JD_Arrive)
		[V_initial, V_final] = Lambert(R_Earth, R_Mars, dt)
		Find mission design parameters and stick them into their appropriate arrays
	End Loop
End Loop
Plot contours in a figure with C3, TOF and V_{\infty}

Pseudo SPICE integration

Possibly generate the ephemeral data needed prior to simulation, eventually simply specify whether to use VSOP or a CSV. Then, use the same format as JPL horizon. Another option is calling the script every time, but that would be incredibly slow and annoying if using third body dynamics.

Naasz control law

Around page 66 of the Hatten MSc. thesis. It's much simpler to implement that Petropoulos as it's similar to Ruggerio.

State Transition Matrix from Orbit and Perturbations

There should be an OrbitEstimation struct which can generate and propagate a state transition matrix from an initial orbit. The reason why I think this fits better here than in gokalman is because more structs of smd will be used in the future (including MissionState to export the results of the estimation). The only thing needed from gokalman is the interfaces for the kalman filters.

Coverage by coverall stalled

Not sure why but it doesn't seem to report the status to Github anymore.

SpiceyPy standalone converter

A simple to use script that converts state vectors between frames.

IMD: HW4

Gravity assist computing function. Should return radius of periapse, turning angle, B-plane target, B vector and theta. Should be a separate file which allows one to compute a number of different gravity assists based on different parameters for plotting.

Lawden's spiral

r = (r_s*(sin alpha)^6)/(1-3*(sin(alpha))^2)
theta = theta_0 - 4*alpha-3*cotan(alpha) -- cotan == cot

Convert to OSS

Once API is more stable and the coverage is greater than 90%, convert this project to open source software. This will require re-organizing the code to make it importable with go get. I'll probably also have to find a new name...

Fix out of plane maneuvers (inclination and RAAN)

Currently there is a dirty dirty hack which does the max thrust of pi/2 for both. This is not scientifically correct and needs to be fixed.

Pure helio simulation

Simulate helio trajectories from Earth to Mars, and from Mars to Earth, by only increasing the velocity and all the way until Mars's orbit (no ephemeris).

Do not use direct optimization for this. Only velocity increase. This will be a core result.

Probably:

• Use go-flags to specify the name of the thrusters and the count of each.
• Start from Earth position and go to Mars position. Then later, input the hyperbolic vectors for a more precise simulation.

Petropoulos test cases

Write the Petropoulos test cases from his 2004 paper, include both the single and multiple orbital element targeting

Support circular equatorial orbits

Cf. page 114 of Vallado. Tests already written.

Differential correction of interplanetary targeting on true anomaly

Make this a command so that scenarios can target a planet correctly. Use a dichotomy method. There is already such an example, but now I want this is a command to make it simple to use. Maybe have it as a tool of smd.

Code Lambert multirev

Also add more Lambert test cases from those provided in homework.

StatOD - Homework 1

Create a numerical simulation of a spacecraft in an orbit with a = 7000 km, e = 0.001, i = 30 degrees, Ω = 80 degrees, ω = 40 degrees, and an initial true anomaly of ν = 0 degrees. Include only μ and J2 in the dynamics for the system. Integrate the system for 24 hours. This is your reference trajectory.
Integrate a second trajectory by perturbing the initial state with a state deviation vector, δx = [1km, 0, 0, 0, 10m/s, 0].
Compare the state of the second trajectory with respect to the reference trajectory over the course of 24 hours. Use the STM computed around the reference trajectory to perform a second propagation of δx.

Add efficiency calculations

Several control laws provide an efficiency calculation so that if the efficiency isn't good enough we can stop thrusting and wait for a better time to continue the maneuver. This is important.

Update for RK4 API changes

Fixing ChristopherRabotin/ode#3 led to API breaking changes.

Side note: it may be of interest to change the update of the time in the SetState function instead of the Stop function.

Add J2 drifts

They are trivial to add (cf. Naasz page 8), but I should probably have some kind cutoff factor since they are negligible starting a given distance.

StatOD: HW4

• Ensure CKF works
• Ensure EKF works
• Add state noise compensation

PlanetTarget waypoint precision

Make PlanetTarget more precise

Spiralfdr should run on multiple threads

I should be able to run simulations in different go routines and with a setting for the max number of simultaneous running routines.

Pinkham's spiral

f_r = (q*sqrt(mu*p)*v_r)/(2*r**2)
f_theta = (q*sqrt(mu*p)*v_theta)/(2*r**2)
p = p_s * exp*(q0)
r = (p*(1+q**1))/(1+exp(q0)*(1+q**2)*k*cos(theta-omega))

Propagate in inertial frame (Cartesian)

This will be essential for propagating restricted circular three body problem (there is no other solution). This may also help solve #35.

Heliocentric inclination is wrong

Example: when switching an orbit from an Earth centric to a heliocentric, the inclination goes from 28 degrees to 23 degrees. This does not seem correct, but I may be wrong in this judgment. More investigation is needed. Clearly, some tilt is needed, as without any tilt, we get the following 0 degree inclination orbit.

Petropoulos control Q-law

Implement the Q-law as per the 2004 paper. May need to rewrite part of the Orbit Target code to allow for both Ruggerio and Petropoulos, which both rely on the same basic control laws and defer only on the summation.

IMD: Mars B Plane targeting

Investigate with changing the semi major axis also changes the eccentricity

That makes no sense at all. There is no out of plane thrust, so the eccentricity should not change at all, and I triple checked the equations in Mission.Func()... WTF?!

case OptiΔaCL:
		ctrl = func(o Orbit) []float64 {
			sinν, cosν := math.Sincos(o.ν)
			return unitΔvFromAngles(math.Atan2(o.e*sinν, 1+o.e*cosν), 0.0)
		}
		break

Targeted spiral

It should be possible to target a spiral in such a way as to reach SOI aligned with the planet, regardless of thrust.

Add estimation Phi tests and J2 tests

• J2 tests: assignment: http://ccar.colorado.edu/ASEN5070/assignments/Fall2015/ASEN5070_HW02.pdf ; solutions: http://ccar.colorado.edu/ASEN5070/solutions/Fall2015/ASEN5070_HW02_IntegCheck.pdf
• Phi tests: http://ccar.colorado.edu/ASEN5070/solutions/Fall2015/ASEN5070_HW04_IntegCheck.pdf (includes assignment)

Arbitrary export of orbital element functions

Possibly allow for the CSV exporter to call any function of the Orbit. Specifically, I need to be able to plot the specific energy of an orbit. This will be useful when/if migrating to a Cartesian propagation (which will be needed for the circular restricted three body problem).

StatOD: HW3

• Orbit simulations: generate noisy measurements for the three ground stations for the LEO orbit from HW2, using J2 and J3
• Run estimators with the following:
  • CKF
  • EKF
  • Batch KF
• Plot errors of position and velocity with 3 sigma bounds (i.e. six plots per KF)
• Plot post fit residuals for all three, should look like white noise
• Answer questions.

Mars Orbit Insertion

Grand finale: perform injection to Mars from interplanetary trajectory. Once this is done, move onto final project.

Orbit struct to store R, V and c instead of OE

OEs break when doing hyperbolic orbits. This may be the main issue TBH, and it definitely is an issue when trying to visualize these kind of orbits in Cosmographia.

TOF based targeting

Use TOF to compute best time to start the launch for a given thrust. Effectively, run with an initial guess from a launch date, then propagate until the target orbit is reached, get the TOF for this, and re-propagate this time targeting the destination orbit at the time of launch plus initial TOF. This may require a few iterations to get right (dichotomy?).

Several levels of precision

• When integrating a new mission, be able to specify the time step
• When generating a new orbit, be able to specify the precision (probably using an Enum). To avoid breaking things, it's probably better to have a new constructor solely for this.

Fix J3 perturbations

J3 perturbations are currently disabled because they are buggy: https://github.com/ChristopherRabotin/smd/blob/master/perturbations.go#L55 . These need to be fixed for HW3 of StatOD.

Propagate in Lagrange Equinoctial

Maybe even transform the Orbit struct to be based on these. They support all kinds of orbits, and make sense when reading the values. Lagrange VOP allow for arbitrary perturbations which is desired in #36. In addition, there are many paper which explain how to handle these for J_n and the continuous thrusting should work out of the box.

Support for circular orbits

Ugh, I need to support circular orbits somehow.

Lambert solver

Add a tools folder that has specific useful tools, including a Lambert solver (using universal variables).

Hohmann transfer

Hohmann transfers must be able to provide time of flight, delta-v, departure and arrival velocities. This must also be supported for patched conics.

Only provide an instantaneous burn initially, with zero fuel consumption. That will make things much easier for the propagation. Anything more isn't yet requested in the homework.

Spiral in/out using Cartesian and just increasing the velocity

The semi major axis optimization isn't necessary the best for this. The vehicle keeps on hitting the right spot but it cannot seem to do the transition correctly from helio to SOI.

Earth GroundStation & Exporting changes

MissionState

• Rename MissionState to State
• Add a NewPureState() which has an empty Spacecraft struct. This will allow State to also be used in the export for OrbitEstimate

ExportConfig

• Add Exporter as interface with methods Init(e ExportConfig) and Export(st State). Add a few exporters including ExportCosmo, ExportOE, ExportSC, ExportStations.
• Add an attribute of type []Exporter to have a list of exporters.

GroundStation

Add a GroundStation struct which has the following attributes and methods:

• Attributes:
  • RangeNoiseMean float64 (will use normal distribution so this is only the mean)
  • RangeRateNoiseMean float64 (idem)
  • Elevation
  • Frequency
• Methods:
  • OnVisible func(st Snapshot)

Mission

• Add computeSTM flag and effectively merge OrbitEstimate and Mission. GMAT has the capability of generating the STM on demand as well, and that makes sense given the similarities.
• Mission should be able to make station measurements too. Possibly could that be done by pushing to a channel, something like Mission.SetMeasurementChannel(->someChan (Measurement)) and then in the calling function just wait on measurements to come through.

Arbitrary perturbations

Celestial Perturbations

The best way to handle this is probably by a "Pertubations" struct which supports J_n of different bodies (probably just stick to J2 to start with) and which other bodies to consider during the orbit. I think that this should only be added in Func() of Mission. Possibly, this could support three bodies initially and determine which body to consider based on the distance to each. We would also log which body is now being considered the third body.

Arbitrary Pertubations

For Stat OD, I'll need to add in random perturbations on specific orbital elements (specifically, HW1 requires a 1km delta on the first component of the position and a 10m/s delta on the first component of the velocity. Seems like this will also require #42. Although, for this, I could also convert the orbit to an RV, then perform the ECI perturbation, and convert this back to Keplerian OE.

Python 2 warning suppression

Cf. http://stackoverflow.com/questions/14463277/how-to-disable-python-warnings