; Leapfrog integration method
; used for dynamical equations of classical mechanics: accelerations = f(positions)
; 
; here, we study one particle in harmonic oscillator potential f = -(k/m)*x, with 
; omega^2  = k/m = 1, in time from t=0 to t=20.
; In leapfrog, at the end of every time step we have velocities
; lagging the positions by dt/2. (Factor q = 0.5 in the first step only, 
; later v is out of sync with x, and q remains equal 1. 
;
; The offset by dt/2 and the correct order of first updating v and then
; using it to update x, makes the method so good. It's 2nd order acurate for x(t), 
; despite only ONE force evaluation per timestep. Compare that with 
; method in inte1.pro, where 2nd order accuracy is reached with 
; TWO r.h.s evaluations. In real world, r.h.s. evaluations can be very 
; time-consuming for high-dimesional systems, like N gravitating bodies.
; For instance, exact computation of gravity requires N*(N-1)/2  = O(N^2) 
; pairwise force evaluations, i.e. on the order of 1e12 distance and force 
; calaculations for a million-body system. That's why we'd be happy 
; to have as few force evaluations per timestep as possible, while keeping 
; (at least) 2nd order of accuracy. There is a symplectic algorithm that
; I would recommend, which achieves 4th order accuracy with only 3 r.h.s.
; evaluations (https://en.wikipedia.org/wiki/Symplectic_integrator)
;
; P. Art Jan 2014-2017
;
t = 0.d0
x = 1.d0
v = 0.d0
dt = 0.001d0
Tfinal = 20.d0
q = 0.5d0  ; first time step needs dt/2 for v update
N = Tfinal/dt
; declare a table of results
res = 1d0*fltarr(N,4)

for i=0,N-1 do begin

	t = t+dt
	f = -x       ; harmonic oscillator, frequency = 1
	v = v + f*dt *q
	x = x + v*dt

;	store results
	res(i,0) = t
	res(i,1) = x
	res(i,2) = v
;	compute and store energy 
	res(i,3) = 0.5d0*(v+dt/2d0*f)^2 + 0.5d0*x^2  ; mechanical energy, notice 
	     ; that we needed to synchronize v to be that at the end of time step
	if (i/10*10 EQ i) then print,i,t,x,v,res(i,3)/res(0,3)-1.
	q = 1.0	   ; programming trick to avoid a slower if(i EQ 1) then..
endfor

; plot all the results, energy (relative) error multiplied by 10^5

plot,res(*,0),res(*,1), xtitle =' time', ytitle='x, v, energy error * 1e5', $
   charsize=2, title=' Math pendulum equation dx/dt = -x with dt=0.001'
oplot,res(*,0),res(*,2) 
oplot,res(*,0),(res(*,3)/res(0,3)-1.d0)*1d5
oplot,res(*,0),0*res(*,0)

end