; leapfrog method
; Hill problem dot(dot(x)) = -3/r3*x +3x +2 dy/dt
; and dot(dot(y)) = -3/r3*y -2 dx/dt 
;
; in leapfrog method, at the end of every time step we have velocities
; lagging the positions by dt/2. 
; q = 0.5 in the first step only, later v is out of sync with x,
; which makes the method to be so good (2nd order in time)
; 2014	P. Art.

adot = 0.0
read,' planet migration speed da/dt = ',adot

;-----loop over init. cond.
for k=1,40 do begin
;-----initial conditions
	x = 3.75-k/15.d0  & y=12.d0
	vx = 0.d0  & vy = -3.d0/2.*x
	r = sqrt(x^2+y^2)  &  r0 = r
	dt0 = 0.001d0   &      dt = dt0
	Tfinal = 50.   &        t = 0d0
	i = -1
	q = 0.5  ; first time step needs dt/2 for v update
	N = floor(Tfinal/0.001*1.2)	; max # of steps
;-------declare a table of results
	res = fltarr(N,6)
;_________________________________________________________________________
while ((t LT Tfinal) AND (r LE r0) AND (i LT N-1)) do begin    ; main time loop
	i = i + 1
	t = t + dt
;------find r and -3/r^3  (r3 := r^3)
	r2 = (x*x +y*y) >0.9999999d-4 ; softening below r~1d-4 enforced to prevent huge f
	r = sqrt(r2)
	r3 = r2*r
;------compute forces
	f = -3./r3
	fx = f*x +3.*x +2.*vy
	fy = f*y -2.*vx -adot/2.
;------do one leapfrog step
	vx = vx + fx*dt *q
	vy = vy + fy*dt *q
	x = x + vx*dt
	y = y + vy*dt
;-----store results
	res(i,0) = t
	res(i,1) = x
	res(i,2) = y
	res(i,3) = vx
	res(i,4) = vy
;-----if very close to the planet, remove far from it so the program stops
	if (r LT 0.01) then x=100.
;-----compute and store energy constant
	vxsynch = vx+dt/2*fx
	vysynch = vy+dt/2*fy ; v's synch'ed with x,y
;-----Jacobi const
	res(i,5) = (vxsynch^2+vysynch^2)/2. -3./r -1.5*x^2 +adot*y/2. 
	q = 1.0d0      ; q de-synchronizes v's and x's in 1st step 
;-----recommend dt for next step
	dt = dt0*r*sqrt(r) 	;  ~r^(3/2), like e.g. the orbital time 
	dt = (dt > 0.000001d0) < 0.01d0  ; minimum and maximum enforced

endwhile
;_____________________________________________________________________________


;----truncate arrays 
	res = res(0:i-2,*)
;----calculate error of Jacobi const.
	energy_err = res(*,5)/res(0,5)-1.
;----plot trajectory on (x,y) plane
	if(k EQ 1)then  plot,res(*,1),res(*,2), $
		xrange=[-3,3],yrange=[-3,3], title='Hill problem trajectories',$
		charsize=2
	if(k GT 1)then oplot,res(*,1),res(*,2)
	;oplot,res(*,1),res(*,2),psym=4,symsize=0.5
	oplot,[-10.,10.],[0.,0.],color=1234 
	oplot,[0.,0.],[-10.,10.],color=1234

endfor ; end of loop over initial conditions
end