# bad N-body integrator

from math import sqrt
import copy

G=6.673e-11
AU=1.49597e11
masses=[1.989e30, 5.972e24]
massy=[0,0,0, 0,-(29800*masses[1]/masses[0]),0,  AU,0,0, 0,29800,0]
DRO=1.2
LEO=6371000+800000 # low Earth orbit is above Earth's surface not Earth's core
zeromass=[DRO*AU,0,0,0,24320,0,  AU,LEO,0, 7455,29800,0]

universe=[masses, massy, zeromass]

dt=60.0

def derivs(masses, massy, zeromass):
  n=len(massy)/6
  v=[0 for i in range(6*n)]
  for i in range(n):
    # x_i dot is xdot_i
    v[6*i]=massy[6*i+3]
    v[6*i+1]=massy[6*i+4]
    v[6*i+2]=massy[6*i+5]
    # xdot_i dot is GMr/r^3
    ddx=0; ddy=0; ddz=0
    for j in range(n):
     if (i!=j):
      rx=massy[6*j]-massy[6*i]
      ry=massy[6*j+1]-massy[6*i+1]
      rz=massy[6*j+2]-massy[6*i+2]
      r=sqrt(rx*rx+ry*ry+rz*rz)
      ddx=ddx+G*masses[j]*rx/(r**3)
      ddy=ddy+G*masses[j]*ry/(r**3)
      ddz=ddz+G*masses[j]*rz/(r**3)
    v[6*i+3]=ddx
    v[6*i+4]=ddy
    v[6*i+5]=ddz
  nz=len(zeromass)/6
  vz=[0 for i in range(6*nz)]
  for i in range(nz):
    for ee in [0,1,2]:
      vz[6*i+ee]=zeromass[6*i+ee+3]
    ddx=0; ddy=0; ddz=0;
    for j in range(n):
      rx=massy[6*j]-zeromass[6*i]
      ry=massy[6*j+1]-zeromass[6*i+1]
      rz=massy[6*j+2]-zeromass[6*i+2]
      r=sqrt(rx*rx+ry*ry+rz*rz)
      ddx=ddx+G*masses[j]*rx/(r**3)
      ddy=ddy+G*masses[j]*ry/(r**3)
      ddz=ddz+G*masses[j]*rz/(r**3)
    vz[6*i+3]=ddx
    vz[6*i+4]=ddy
    vz[6*i+5]=ddz
   
  return [v,vz]

f=open("pingle-RK","w")
T = 0
for q in range(1000000):
 olduniverse=copy.deepcopy(universe)
 # distance Earth->Sun
 ug = sqrt((universe[1][6]-universe[1][0])**2+(universe[1][7]-universe[1][1])**2+(universe[1][8]-universe[1][2])**2)
 # two velocities
 v1 = sqrt(sum(t**2 for t in universe[1][3:5]))
 v2 = sqrt(sum(t**2 for t in universe[1][9:11]))
 towrite=[T,ug,v1,v2,universe[1][9],universe[1][10],universe[1][6],universe[1][7],universe[2][0],universe[2][1],universe[2][6],universe[2][7]]
 towrite=[T,sqrt((universe[2][6]-universe[1][6])**2+(universe[2][7]-universe[1][7])**2)]
 for u in towrite:
  f.write(str(u))
  f.write(" ")
 f.write("\n")

 # runge-kutta
 midpoint=copy.deepcopy(universe)
 # k1 based on f(y0)
 Ds1 = derivs(universe[0],universe[1],universe[2])

 D=Ds1[0]; Dz=Ds1[1]
 for i in range(len(D)):
  midpoint[1][i] = universe[1][i] + dt/2 * D[i]
 for i in range(len(Dz)):
  midpoint[2][i] = universe[2][i] + dt/2 * Dz[i]

 # k2 based on f(y0 + k1/2)
 Ds2 = derivs(midpoint[0],midpoint[1],midpoint[2])
 D=Ds2[0]; Dz=Ds2[1]
 for i in range(len(D)):
  midpoint[1][i] = universe[1][i] + dt/2 * D[i]
 for i in range(len(Dz)):
  midpoint[2][i] = universe[2][i] + dt/2 * Dz[i]

 # k3 based on f(y0 + k2/2)
 Ds3 = derivs(midpoint[0],midpoint[1],midpoint[2])
 D=Ds3[0]; Dz=Ds3[1]
 for i in range(len(D)):
  midpoint[1][i] = universe[1][i] + dt * D[i]
 for i in range(len(Dz)):
  midpoint[2][i] = universe[2][i] + dt * Dz[i]

  # k4 based on f(y0+k3)
 Ds4 = derivs(midpoint[0],midpoint[1],midpoint[2])

# print "Ds1=%s\nDs2=%s\nDs3=%s\nDs4=%s\n\n"%(Ds1, Ds2, Ds3, Ds4)

 for slice in [0,1]:
  for i in range(len(Ds1[slice])):
   universe[slice+1][i] = universe[slice+1][i] + dt * (Ds1[slice][i]+2*Ds2[slice][i]+2*Ds3[slice][i]+Ds4[slice][i])/6.0
 
 if (universe[1][3]>0 and olduniverse[1][3]<0):
  print "X velocity of planet changes sign at ",T
 if (universe[2][3]>0 and olduniverse[2][3]<0):
  print "X velocity of satellite changes sign at ",T
 T=T+dt