; 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.
npart = 200
tab = dblarr(npart,5) ; results

for k=0,npart-1 do begin
;-----initial conditions
	ystart=12.
; separate calculations zooming in on interesting regions
; xmin is denoted as x0(left) in plots; the left-hand value of x0
; xinterv  is the range of x0 plotted on horizontal axis
	xmin = 1.7 & xinterv = 0.5           ; plots hill3-zoom1*png, large-scale 
	xmin = 2.148 & xinterv = 0.003 
	xmin = 2.1488 & xinterv = 0.0002  
	xmin = 2.14884d0 & xinterv = 3d-6 
	xmin = 2.1488405d0 & xinterv = 1d-7  
	xmin = 2.14884057d0 & xinterv = 1d-8   ; plots hill3-zoom18*png, micro-scale 

	x = xmin + xinterv*k/npart
	y = ystart
	vx = 0.d0  & vy = -3.d0/2.*x
	r = sqrt(x^2+y^2)  &  r0 = r
	dt0 = 0.0001d0   &      dt = dt0
	Tfinal = 60.   &        t = 0d0
	i = -1L
	q = 0.5  ; first time step needs dt/2 for v update
	N = floor(Tfinal/dt0*4) 	; max # of steps
	epsilon = 0.5d-5
;-------declare a table of results
	res = dblarr(N*1L,6)
;______________________________________________________________________________
while ((t LT Tfinal) AND (abs(y) LE ystart) AND (i LT N-1)) do begin ;time loop
	i = i + 1
	t = t + dt
;------find r and -3/r^3  (r3 := r^3)
	r2 = (x*x +y*y) > epsilon ; softening at r< ~1d-5  to prevent huge f
	r = sqrt(r2)
	r3 = r2*r
;------compute forces
	f = -3d0/r3
	fx = f*x +3.d0*x +2.d0*vy
	fy = f*y -2.d0*vx -adot/2.d0
;------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
;-----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>epsilon) -1.5*x^2 +adot*y/2. 
	q = 1.0d0      ; q de-synchronizes v's and x's in 1st step 
;-----recommend dt for the next step
	dt = dt0*r*sqrt(r) 	;  ~r^(3/2), like e.g. the orbital time 
	dt = (dt > 1d-7) < 0.003d0  ; minimum and maximum enforced

endwhile
;__________________________________________________________________________
next:
;---record final x
	tab(k,0) = res(0,1) ; x0
	tab(k,1) = x
    tab(k,2) = y
	tab(k,3) = t
	tab(k,4) = max(abs(res(*,5)))
;----truncate array for plotting
	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/4*4 EQ k) then begin
	scope = (3*xinterv < 3.) > 1.25
	if(k EQ 0)then  plot,res(*,1),res(*,2), $
		xrange=[-scope,scope],yrange=[-scope,scope], $
		title='Hill problem trajectories dt0 = '+string(dt0),$
		charsize=2
	if(k GT 0)then oplot,res(*,1),res(*,2)
	;oplot,res(*,1),res(*,2),psym=4,symsize=0.25
	oplot,[-10.,10.],[0.,0.],color=1234 
	oplot,[0.,0.],[-10.,10.],color=1234
   endif
;------warn about insufficient simul. time
	if (t GT 0.99*Tfinal) then print,' k=',k,' t(end)=',t

endfor ; end of loop over initial conditions

stop 

plot,tab(*,0)-tab(0,0),tab(*,1),charsize=2,thick=1,$
	title='final vs. initial x.  x0(left)='+string(tab(0,0)),xtitle='x0-x0(left)',$
    ytitle='x(end)
oplot,tab(*,0)-tab(0,0),tab(*,1)*0,line=1
end