Tyear   = 3600.*365.*24.
nyears  = 1
nsteps  = long(10000)
dt =       nyears*tyear/(nsteps)

M_sun =   double(1.999e30)
m_earth = double(6e24)
G =       double(6.67e-11)
R_earth = double(1.5e11)

fact = .9

V_earth =fact* 2.*!pi*r_earth / tyear

L= double(m_earth * v_earth * R_earth)

r =        dblarr(nsteps)
rdot =     dblarr(nsteps)
theta =    dblarr(nsteps)
thetadot = dblarr(nsteps)

t =  findgen(nsteps)*dt


r(0) = r_earth & theta(0) = 0.
rdot(0)= 0.0
thetadot(0) =  L/m_earth/r(0)^2  

   ;for leapfrog, velocity is at half time step

for i = 1,nsteps -1 do begin

  r(i) = r(i-1) + rdot(i-1)*dt 
  rdot(i) = rdot(i-1) + dt * (  L^2 / m_earth^2 / r(i)^3  - G * M_sun / r(i)^2  )

  ;thetadot is a half time step 

  thetadot(i) = L/m_earth/r(i)^2  
  theta(i) = theta(i-1) + dt* thetadot(i)

endfor

!p.charsize = 1
!p.multi = [0,2,2]

plot,r*cos(theta)/r_earth  , r*sin(theta)/r_earth,xtitle = 'X (au)',ytitle = 'Y (au)',$
xrange = [-1.,1.],yrange = [-1.,1.]

;KE = 0.5* (rdot^2 + r^2*thetadot^2)
KE = 0.5*(rdot^2 + L^2 / m_earth^2 / r^2)
PE = -1.*G*m_sun/r
Etot = Ke + PE
plot,t/tyear,Etot,xtitle = 'time (years)',title = 'Energy',ytitle = 'Joules',yrange = [min(pe),max(ke)]
oplot,t/tyear,kE,line = 1
oplot,t/tyear,PE,line = 2

plot,t/tyear,(Etot-Etot(1))/abs(Etot(1)),$
          xtitle = 'time (years)', $
          ytitle = 'fraction',title = 'energy conservation'


plot,r,theta,title = 'Orbit (r,theta)',ytitle = 'rad',xtitle = 'm'

print,'major axis',  (max(r) + min(r))/2., ' m'


end