On Tue, Jan 14, 2014 at 1:40 PM, Grace Roberts <xgrace_rober...@hotmail.co.uk> wrote: > > -For certain initial conditions the programme displays impossible orbits, > showing the satellite making immediate sharp turns or jumping up and down in > velocity. The exact same code can also end up displaying different graphs > when run multiple times. For example when initial conditions are set at: > xx0=[1000.,10000.,1000.,10000.].
The trajectory is decaying into the planet. In real life it hits the surface. In the model it continues until the radius is 0. As the radius approaches zero, odeint will have problems converging. I ran your calculation for Phobos, assuming a perfectly circular orbit. I started with `Rsx == Rxy` and `Vsx = -Vsy`. It steps through a cycle in the expected time. https://en.wikipedia.org/wiki/Phobos_%28moon%29 Output Plot: http://i.imgur.com/Jh0sbnT.png import numpy as np from scipy import integrate import matplotlib.pyplot as plt G = 6.674e-11 Mm = 6.4185e23 # Mars mass Rs = 9.376e6 # Phobos orbit semi-major axis Vs = 2.138e3 # Phobos orbit average velocity Ts = 2.7554e4 # Phobos orbit period def f(y, t): rx, ry, vx, vy = y r = np.hypot(rx, ry) ax = -G * Mm * rx / r ** 3 ay = -G * Mm * ry / r ** 3 return vx, vy, ax, ay Rsx = Rsy = (Rs ** 2 / 2) ** 0.5 Vsy = (Vs ** 2 / 2) ** 0.5 Vsx = -Vsy y0 = [Rsx, Rsy, Vsx, Vsy] t = np.linspace(0, Ts, Ts * 10) y, info = integrate.odeint(f, y0, t, full_output=1) rx, ry, vx, vy = y.T plt.subplot(211) plt.title('Position') plt.plot(t, rx, t, ry) plt.grid(); plt.axis('tight') plt.subplot(212) plt.title('Velocity') plt.plot(t, vx, t, vy) plt.grid(); plt.axis('tight') plt.show() _______________________________________________ Tutor maillist - Tutor@python.org To unsubscribe or change subscription options: https://mail.python.org/mailman/listinfo/tutor
