plotting orbital trajectories in python
Question: How can I setup the three body problem in python? How to I define the function
to solve the ODEs?

The three equations are  
`x'' = -mu / np.sqrt(x ** 2 + y ** 2 + z ** 2) * x`,  
`y'' = -mu / np.sqrt(x ** 2 + y ** 2 + z ** 2) * y`, and  
`z'' = -mu / np.sqrt(x ** 2 + y ** 2 + z ** 2) * z`.

Written as 6 first order we have

`x' = x2`,

`y' = y2`,

`z' = z2`,

`x2' = -mu / np.sqrt(x ** 2 + y ** 2 + z ** 2) * x`,

`y2' = -mu / np.sqrt(x ** 2 + y ** 2 + z ** 2) * y`, and

`z2' = -mu / np.sqrt(x ** 2 + y ** 2 + z ** 2) * z`

I also want to add in the path Plot o Earth's orbit and Mars which we can
assume to be circular. Earth is `149.6 * 10 ** 6` km from the sun and Mars
`227.9 * 10 ** 6` km.

    
    
    #!/usr/bin/env python                                                             
    #  This program solves the 3 Body Problem numerically and plots the trajectories      
    
    import pylab
    import numpy as np
    import scipy.integrate as integrate
    import matplotlib.pyplot as plt
    from numpy import linspace
    
    mu = 132712000000  #gravitational parameter
    r0 = [-149.6 * 10 ** 6, 0.0, 0.0]
    v0 = [29.0, -5.0, 0.0]
    dt = np.linspace(0.0, 86400 * 700, 5000)  # time is seconds
    


Answer: As you've shown, you can write this as a system of six first-order ode's:

    
    
    x' = x2
    y' = y2
    z' = z2
    x2' = -mu / np.sqrt(x ** 2 + y ** 2 + z ** 2) * x
    y2' = -mu / np.sqrt(x ** 2 + y ** 2 + z ** 2) * y
    z2' = -mu / np.sqrt(x ** 2 + y ** 2 + z ** 2) * z
    

You can save this as a vector:

    
    
    u = (x, y, z, x2, y2, z2)
    

and thus create a function that returns its derivative:

    
    
    def deriv(u, t):
        n = -mu / np.sqrt(u[0]**2 + u[1]**2 + u[2]**2)
        return [u[3],      # u[0]' = u[3]
                u[4],      # u[1]' = u[4]
                u[5],      # u[2]' = u[5]
                u[0] * n,  # u[3]' = u[0] * n
                u[1] * n,  # u[4]' = u[1] * n
                u[2] * n]  # u[5]' = u[2] * n
    

Given an initial state `u0 = (x0, y0, z0, x20, y20, z20)`, and a variable for
the times `t`, this can be fed into `scipy.integrate.odeint` as such:

    
    
    u = odeint(deriv, u0, t)
    

where `u` will be the list as above. Or you can unpack `u` from the start, and
ignore the values for `x2`, `y2`, and `z2` (you must transpose the output
first with `.T`)

    
    
    x, y, z, _, _, _ = odeint(deriv, u0, t).T
    

