If you want to be really accurate, you should use a programming language to simulate the trajectories. From the examples, the eccentricities are off in the solutions above. Here is Python code I use all the time to create this kind of stuff. Note: I work in orbital mechanics.
Here are some notes I am working on that may help you if you want to code it up. This notes discuss orbital mechanics. They aren't done yet but with the link you can always check to see when it is complete.
Orbital Mechanics Notes
Here is the Earth and Mars in astronomical units. You can add_patches for the spheres for the Sun, Earth, and Moon. The code just needs to be changed to fit your problem.
#!/usr/bin/env ipython
# This program solves the 3 Body Problem numerically and plots the
# trajectories
import numpy as np
from scipy.integrate import odeint
from mpl_toolkits.mplot3d import Axes3D
import pylab
mu = 1.0
# r0 = [-149.6 * 10 ** 6, 0.0, 0.0] # Initial position
# v0 = [-5.04769, -29.9652, 0.0] # Initial velocity
# u0 = [-149.6 * 10 ** 6, 0.0, 0.0, 29.9652, -5.04769, 0.0]
u0 = [-1.0, 0.0, 0.0, 0.169474, -1.0067, 0.0]
e0 = [-1.0, 0.0, 0.0, 0.0, -1.0, 0.0]
m0 = [1.53, 0.0, 0.0, 0.0, 1.23152, 0.0]
def deriv2(e, dt):
n = -mu / np.sqrt(e[0] ** 2 + e[1] ** 2 + e[2] ** 2)
return [e[3], # dotu[0] = u[3]'
e[4], # dotu[1] = u[4]'
e[5], # dotu[2] = u[5]'
e[0] * n, # dotu[3] = u[0] * n
e[1] * n, # dotu[4] = u[1] * n
e[2] * n] # dotu[5] = u[2] * n
def deriv(u, b):
n = -mu / np.sqrt(u[0] ** 2 + u[1] ** 2 + u[2] ** 2)
return [u[3], # dotu[0] = u[3]'
u[4], # dotu[1] = u[4]'
u[5], # dotu[2] = u[5]'
u[0] * n, # dotu[3] = u[0] * n
u[1] * n, # dotu[4] = u[1] * n
u[2] * n] # dotu[5] = u[2] * n
def deriv3(m, t):
n = -mu / np.sqrt(m[0] ** 2 + m[1] ** 2 + m[2] ** 2)
return [m[3], # dotu[0] = u[3]'
m[4], # dotu[1] = u[4]'
m[5], # dotu[2] = u[5]'
m[0] * n, # dotu[3] = u[0] * n
m[1] * n, # dotu[4] = u[1] * n
m[2] * n] # dotu[5] = u[2] * n
b = np.arange(0.0, 2 * np.pi, .01)
dt = np.arange(0.0, 2 * np.pi, .01)
t = np.arange(0.0, 2.5 * np.pi, .01)
u = odeint(deriv, u0, b)
e = odeint(deriv2, e0, dt)
m = odeint(deriv3, m0, t)
x, y, z, x2, y2, z2 = u.T
x3, y3, z3, x4, y4, z5 = e.T
x6, y6, z6, x7, y7, z7 = m.T
fig = pylab.figure()
ax = fig.add_subplot(111, projection='3d')
ax.plot(x3, y3, z3)
ax.plot(x, y, z)
ax.plot(x6, y6, z6)
pylab.axis((-2, 2, -2, 2))
pylab.show()
Ignore the u trajectory. The part that will aid you is the e for earth and m for mars. But as I said, you will need to make the appropriate adjustments.
Ignore the green path.
Here is the image this creates:

Here is a plot with a 3d sphere for you to look at too:
