; 
; Lorenz system
  sigma = 8d0 ; 10d0
  b = 3d0 ; 8d0/3d0
  r0 = 28d0
; r = variable parameter
;
; 
; 2016	P. Art.

read,' r0 = ',r0

;-----loop over init. cond.
for k=1,40 do begin
	r = r0 *(1+0.02d0*k)
;-----initial conditions
	x = 1. & y = 1. & z = 1.
	dt0 = 0.001d0   &      dt = dt0
	Tfinal = 200.   &        t = 0d0
	i = -1L
	N = Tfinal/dt +1	; max # of steps
;-------declare a table of results
	res = fltarr(N,6)

;______________________________________________________________________________
while ((t LT Tfinal) AND (i LT N-1)) do begin    ; main time loop
	i = i + 1
	t = t + dt
;------evaluate prediction
	dxdt0 = sigma*(y - x)
	dydt0 = r*x - y - x*z
	dzdt0 = x*y - b*z
	xp = x + dxdt0*dt
	yp = y + dydt0*dt
	zp = z + dzdt0*dt
;------evaluate predicted r.h.s. at the end of timestep
	dxdt1 = sigma*(yp - xp)
	dydt1 = r*xp - yp - x*zp
	dzdt1 = xp*yp - b*zp
;------perform update averaging 0 and 1 r.h.s prediction
	dxdt =  0.5d0*(dxdt0 + dxdt1)
	dydt =  0.5d0*(dydt0 + dydt1)
	dzdt =  0.5d0*(dzdt0 + dzdt1)
; 
	x = x + dxdt*dt
	y = y + dydt*dt
	z = z + dzdt*dt
;-----store results
	res(i,0) = t
	res(i,1) = x
	res(i,2) = y
	res(i,3) = z
	res(i,4) = x-xp
;----print
;	if (i/500*500 EQ i) then begin
;	print,t,x,y,z
;	endif
endwhile
;_____________________________________________________________________________


;----truncate arrays 
	res = res(0:i-2,*)
 
;----plot trajectory on (x,y) plane
	ti = '(y,z) Lorentz system, r='+string(r)
	plot,res(*,2),res(*,3), title=ti,$
;		xrange=[-20,20],yrange=[-20,20], title='(y,z) Lorentz system',$
		charsize=2
	oplot,[-300.,300.],[0.,0.],color=1234 
	oplot,[0.,0.],[-300.,300.],color=1234
	wait,1.

;----Poincare map x=0
	m=0L
	for i=0L,N-6 do begin
	if (res(i,2)*res(i+1,2) LT 0.) then m=[m,i]  ; y=0 crossing
	endfor
	m(0)=m(1)
	;print,m

;----plot Poincare map on (x,z) plane
	poincx = res(m,1) & poincy = res(m,3)
	ti = '(x,z) Poincare map, r='+string(r)
	plot,poincx,poincy,title=ti,charsize=2,psym=4,symsize=2
	wait,3.


endfor ; end of loop over initial conditions
end