Earth-Mars Free Return

Earth-Mars Free ReturnLet’s try using the Fortran Astrodynamics Toolkit and Pikaia to solve a real-world orbital mechanics problem. In this case, computing the optimal Earth-Mars free return trajectory in the year 2018. This is a trajectory that departs Earth, and then with no subsequent maneuvers, flies by Mars and returns to Earth. The building blocks we need to do this are:

  1. An ephemeris of the Earth and Mars [DE405]
  2. A Lambert solver [1]
  3. Some equations for hyperbolic orbits [given below]
  4. An optimizer [Pikaia]

The \(\mathbf{v}_{\infty}\) vector of a hyperbolic planetary flyby (using patched conic assumptions) is simply the difference between the spacecraft’s heliocentric velocity and the planet’s velocity:

  • \(\mathbf{v}_{\infty} = \mathbf{v}_{heliocentric} – \mathbf{v}_{planet}\)

The hyperbolic turning angle \(\delta\) is the angle between \(\mathbf{v}_{\infty}^{-}\) (the \(\mathbf{v}_{\infty}\) before the flyby) and the \(\mathbf{v}_{\infty}^{+}\) (the \(\mathbf{v}_{\infty}\) vector after the flyby). The turning angle can be used to compute the flyby periapsis radius \(r_p\) using the equations:

  • \(e = \frac{1}{\sin(\delta/2)}\)
  • \(r_p = \frac{\mu}{v_{\infty}^2}(e-1)\)

So, from the Lambert solver, we can obtain Earth to Mars, and Mars to Earth trajectories (and the heliocentric velocities at each end). These are used to compute the \(\mathbf{v}_{\infty}\) vectors at: Earth departure, Mars flyby (both incoming and outgoing vectors), and Earth return.

There are three optimization variables:

  • The epoch of Earth departure [modified Julian date]
  • The outbound flight time from Earth to Mars [days]
  • The return flyby time from Mars to Earth [days]

There are also two constraints on the heliocentric velocity vectors before and after the flyby at Mars:

  • The Mars flyby must be unpowered (i.e., no propulsive maneuver is performed). This means that the magnitude of the incoming (\(\mathbf{v}_{\infty}^-\)) and outgoing (\(\mathbf{v}_{\infty}^+\)) vectors must be equal.
  • The periapsis flyby radius (\(r_p\)) must be greater than some minimum value (say, 160 km above the Martian surface).

For the objective function (the quantity that is being minimized), we will use the sum of the Earth departure and return \(\mathbf{v}_{\infty}\) vector magnitudes. To solve this problem using Pikaia, we need to create a fitness function, which is given below:

 subroutine obj_func(me,x,f)

 !Pikaia fitness function for an Earth-Mars free return

 use fortran_astrodynamics_toolkit

 implicit none

    class(pikaia_class),intent(inout) :: me !pikaia class
    real(wp),dimension(:),intent(in)  :: x  !optimization variable vector
    real(wp),intent(out)              :: f  !fitness value

    real(wp),dimension(6) :: rv_earth_0, rv_mars, rv_earth_f
    real(wp),dimension(3) :: v1,v2,v3,v4,vinf_0,vinf_f,vinfm,vinfp
    real(wp) :: jd_earth_0,jd_mars,jd_earth_f,&

    real(wp),parameter :: mu_sun      = 132712440018.0_wp  !sun grav param [km3/s2]
    real(wp),parameter :: mu_mars     = 42828.0_wp         !mars grav param [km3/s2]
    real(wp),parameter :: rp_mars_min = 3390.0_wp+160.0_wp !min flyby radius at mars [km]
    real(wp),parameter :: vinf_weight = 1000.0_wp          !weights for the
    real(wp),parameter :: rp_weight   = 10.0_wp            ! penalty functions

    !get times:
    jd_earth_0      = mjd_to_jd(x(1))           !julian date of earth departure
    dt_earth_out    = x(2)                      !outbound flyby time [days]
    dt_earth_return = x(3)                      !return flyby time [days]
    jd_mars         = jd_earth_0 + dt_earth_out !julian date of mars flyby
    jd_earth_f      = jd_mars + dt_earth_return !julian date of earth arrival

    !get earth/mars locations (wrt Sun):
    call get_state(jd_earth_0,3,11,rv_earth_0) !earth at departure
    call get_state(jd_mars,   4,11,rv_mars)    !mars at flyby
    call get_state(jd_earth_f,3,11,rv_earth_f) !earth at arrival

    !compute lambert maneuvers:
    call lambert(rv_earth_0,rv_mars,   dt_earth_out*day2sec,   mu_sun,v1,v2) !outbound
    call lambert(rv_mars,   rv_earth_f,dt_earth_return*day2sec,mu_sun,v3,v4) !return

    !compute v-inf vectors:
    vinf_0 = v1 - rv_earth_0(4:6) !earth departure
    vinfm  = v2 - rv_mars(4:6)    !mars flyby (incoming)
    vinfp  = v3 - rv_mars(4:6)    !mars flyby (outgoing)
    vinf_f = v4 - rv_earth_f(4:6) !earth return

    !the turning angle is the angle between vinf- and vinf+
    delta   = angle_between_vectors(vinfm,vinfp) !turning angle [rad]
    vinfmag = norm2(vinfm)                       !vinf vector mag (incoming) [km/s]
    e       = one/sin(delta/two)                 !eccentricity
    rp      = mu_mars/(vinfmag*vinfmag)*(e-one)  !periapsis radius [km]

    !the constraints are added to the fitness function as penalties:
    if (rp>=rp_mars_min) then
        rp_penalty = zero    !not active
        rp_penalty = rp_mars_min - rp
    end if
    vinf_penalty = abs(norm2(vinfm) - norm2(vinfp))

    !fitness function (to maximize):
    f = - ( norm2(vinf_0) +            &
            norm2(vinf_f) +            &
            vinf_weight*vinf_penalty + &
            rp_weight*rp_penalty       )

 end subroutine obj_func

The lambert routine called here is simply a wrapper to solve_lambert_gooding that computes both the “short way” and “long way” transfers, and returns the one with the lowest total \(\Delta v\). The two constraints are added to the objective function as penalties (to be driven to zero). Pikaia is then initialized and called by:

 program flyby

 use pikaia_module
 use fortran_astrodynamics_toolkit

 implicit none

   integer,parameter :: n = 3

   type(pikaia_class) :: p
   real(wp),dimension(n) :: x,xl,xu
   integer  :: status,iseed
   real(wp) :: f

   !set random number seed:
   iseed = 371

   xl = [jd_to_mjd(julian_date(2018,1,1,0,0,0)),  100.0_wp,100.0_wp]
   xu = [jd_to_mjd(julian_date(2018,12,31,0,0,0)),400.0_wp,400.0_wp]

   call p%init(n,xl,xu,obj_func,status,&
            ngen               = 1000,      &
            np                 = 200,       &
            nd                 = 9,         &
            ivrb               = 0,         &
            convergence_tol    = 1.0e-6_wp, &
            convergence_window = 200,       &
            initial_guess_frac = 0.0_wp,    &
            iseed              = iseed      )

   call p%solve(x,f,status)

 end program flyby

The three optimization variables are bounded: the Earth departure epoch must occur sometime in 2018, and the outbound and return flight times must each be between 100 and 400 days.

Running this program, we can get different solutions (depending on the value set for the random seed, and also the Pikaia settings). The best one I’ve managed to get with minimal tweaking is this:

  • Earth departure epoch: MJD 58120 (Jan 2, 2018)
  • Outbound flight time: 229 days
  • Return flight time: 271 days

Which is similar to the solution shown in [2].


  1. R.H. Gooding, “A procedure for the solution of Lambert’s orbital boundary-value problem“, Celestial Mechanics and Dynamical Astronomy, Vol. 48, No. 2, 1990.
  2. Dennis Tito’s 2021 Human Mars Flyby Mission Explained., February 27, 2014.
  3. J.E. Prussing, B.A. Conway, “Orbital Mechanics”, Oxford University Press, 1993.
Tagged with: , , , , ,
6 comments on “Earth-Mars Free Return
  1. mjburgess says:

    I came across this program and your fortran astrodynamics toolkit while looking for ephemeris routines to use in conjunction with PyKEP.

    The flyby program is a very nice demonstration of the facilities offered by the toolkit and the purpose of this comment is to congratulate you on the excellent work and I do hope that you will continue to expand the toolkit and receive recognition for your efforts.

    I had the toolkit up and running in about half an hour and the flyby program in about an hour. Also access to the JPL ephemeris files from Python using the toolkit and F2PY to build the required library.

    This really shows the power of modern Fortran with the source code being readily intelligible and the interoperability not far short of what could be achieved with c++.

    so thank you and 10/10!

  2. Jacob says:

    Thanks! Glad you found it useful.

  3. mjburgess says:

    I see you have been very busy in your continuing one-man crusade to revitalise Fortran and the later part of 2015 has seen no decrease in your output. I wish you all success in this venture.

    I have continued to unearth gems deep in your fortran astrodynamics toolkit and whilst working on a trajectory propagation algorithm found that all required tools were already in place. In the process of writing a trajectory propagation module (before I discovered yours) I came across many ancient ODE solvers in Fortran such as ODEPACK which I see you are now hosting on your github account.

    I also note the recent release of your updated DDEABM. If you are turning your attention back to astrodynamics I have been attempting to write a low thrust trajectory optimiser which can use bracketed root finding so perhaps you could apply your updated version to this task. Other functions that you could usefully add to the astrodynamics toolkit would be low thrust propagation, sun angle routines for calculating solar panel output and gravity assist / B-plane routines.

  4. Jacob says:

    Yes indeed, I do want to include these algorithms (I don’t know when I’ll get to it, though). If you have anything you want to contribute via GitHub, then I’ll welcome it. Low-thrust is certainly something I want to include. Note: I did just add the beginnings of a B-plane module.

    For high-accuracy trajectory propagations, I do recommend DDEABM or DOP853, and also ODEPACK (DLSODE or DVODE). All are good (the first two I have modernized, the others not…yet). Also, by all accounts, DIVA from MATH77 is pretty good, but I’ve not used it yet (it only recently was released as open source). ODEPACK, DIVA, and my updated DDEABM all have bracketed root-finding capabilities.

  5. Oscar says:

    I am working with a 3D problem and for memory requirements the parameters to be optimized are of type real*4, What would I have to change to use the library?

  6. Jacob says:

    By real*4 I presume you mean single precision (note that the real*4 syntax is not standard Fortran). If so, all you need to do is change the declaration of wp in kind_module to wp = real32.

Leave a Reply

Your email address will not be published. Required fields are marked *