c
c Simulation of the Lorenz map
c        - Euler or Runge-Kutta (2nd order usual one)
c Fortran version written by H. Nakanishi, need to be compiled with "-lpepl".
c Line 89 was originally incorrect, corrected to "t1 = t(i-1)+0.5*dt"
c
	program lorenz
c
c       Declare the arrays we will need
c
	implicit real*8 (a-h,o-z)
	dimension x(500003),y(500003),z(500003),t(500003)
c
c       Use subroutines to do the work
c
	call initialize(x,y,z,t,sigma,r,b,dt,lsym,nsym,k,nstep,m,nout)
	call calculate(x,y,z,t,sigma,r,b,dt,lsym,nsym,k,nstep,m)
	if(nout.eq.1) then
		open(1,file='out.lorenz')
		rewind(1)
		do 10 i=1,nstep
			write(1,15) t(i),x(i),y(i),z(i)
10		continue
		close(1)
15		format(1x,1p,4(e12.5,1x))
	else
		call display(x,y,z,t,sigma,r,b,dt,lsym,nsym,k,nstep,m)
	endif
	stop
	end
c
	subroutine initialize(x,y,z,t,sigma,r,b,dt,lsym,nsym,k,n,m,nout)
c
c       Initialize variables
c
	implicit real*8 (a-h,o-z)
        dimension x(1),y(1),z(1),t(1)
	character ans,yes
	yes='y'
	print *,'Lorenz model: Euler (1) or Runge-Kutta [2nd ord] (2)?'
	read(5,*) m
	print *,'parameters: sigma, r, b ?'
	read(5,*) sigma,r,b
	print *,'initial x, y, z ?'
	read(5,*) x(1),y(1),z(1)
	print *,'time step, how many steps ?'
	read(5,*) dt,n
	t(1)=0
	print *,'record output to a file [out.lorenz] ?'
	read(5,14) ans
	if(ans.eq.yes) then
		nout=1
	else
		nout=2
		print *,'plot [z vs t,z vs x](1) or [x , y vs t](2)?'
		read(5,*) k
		if(k.ne.1.and.k.ne.2) then
			print *,'must select 1 or 2 ...'
			stop
		endif
		print *,'set line, symbol?'
		read(5,14) ans
14		format(a1)
		if(ans.eq.yes) then
			print *,'lsym, nsym ? -> '
			read(5,*) lsym,nsym
		else
			lsym=0
			nsym=1
		endif
	endif
	return
	end
c
	subroutine calculate(x,y,z,t,sigma,r,b,dt,lsym,nsym,k,nstep,m)
c
c       Now use Euler or Runge-Kutta (2nd order)
c
	implicit real*8 (a-h,o-z)
	dimension x(1),y(1),z(1),t(1)
	nstep=min(nstep,500000)
	if(m.eq.2) then
	   do 10 i = 2,nstep
		t(i) = t(i-1)+dt
		call dv(x(i-1),y(i-1),z(i-1),t(i-1),sigma,r,b,dx,dy,dz)
         	x1 = x(i-1)+0.5*dt*dx
         	y1 = y(i-1)+0.5*dt*dy
         	z1 = z(i-1)+0.5*dt*dz
		t1 = t(i-1)+0.5*dt
		call dv(x1,y1,z1,t1,sigma,r,b,dx,dy,dz)
		x(i)=x(i-1)+dt*dx
		y(i)=y(i-1)+dt*dy
		z(i)=z(i-1)+dt*dz
10	   continue
	else
	   do 20 i = 2,nstep
		t(i) = t(i-1)+dt
		call dv(x(i-1),y(i-1),z(i-1),t(i-1),sigma,r,b,dx,dy,dz)
         	x(i) = x(i-1)+dt*dx
         	y(i) = y(i-1)+dt*dy
         	z(i) = z(i-1)+dt*dz
20	   continue
	endif
	return
	end
c
	subroutine dv(x0,y0,z0,t0,sigma,r,b,dx,dy,dz)
	implicit real*8 (a-h,o-z)
	dx=sigma*(y0-x0)
	dy=-x0*z0+r*x0-y0
	dz=x0*y0-b*z0
	return
	end
c
	subroutine display(x,y,z,t,sigma,r,b,dt,lsym,nsym,k,n,m)
c
c       First set up title and label axes for graph. Plotting is a lot
c       of work in fortran.
c       This version displays output as well as writes to a file "graph.out".
c
	implicit real*8 (a-h,o-z)
	dimension x(1),y(1),z(1),t(1)
	real u(n+3),v(n+3),fac,pos,s1,r1,b1,dt1
	s1=sigma
	r1=r
	b1=b
	dt1=dt
	fac=0.8
	pos=10.2
	call usrmon(.true.,.false.,-0.5,9.,-0.5,12.)
	call pltlun(19,.true.,.false.)
	call pltlfn('graph.out')
	call plots
	call factor(fac)
	call plot(0.,1.,-3)
	if(m.eq.2) then
		call text(0.2,pos,0.2,'Lorenz: Runge-Kutta2',0.,20,0)
	else
		call text(0.2,pos,0.2,'Lorenz: Euler',0.,13,0)
	endif
	if(k.eq.1) then
		do 15 i=1,n
			u(i)=x(i)
			v(i)=z(i)
15		continue
		call text(0.2,4.0,0.2,'z vs x',0.,6,0)
		call scalex(u,6.,n,1)
		call axctl(0.15,0.04,0.15,0.2,0.2,-1)
		call axisx(0.,0.,'x',-1,6.,0.,u(n+1),u(n+2),u(n+3),4)
		call axisx(0.,3.5,' ',1,6.,0.,u(n+1),u(n+2),u(n+3),20)
		call scalex(v,3.5,n,1)
		call axctl(0.15,0.04,0.15,0.2,0.2,2)
		call axisx(0.,0.,'z',1,3.5,90.,v(n+1),v(n+2),v(n+3),-5)
	call axisx(6.,0.,' ',-1,3.5,90.,v(n+1),v(n+2),v(n+3),-21)
		call line(u,v,n,1,lsym,nsym)
		call number(0.5,3.7,0.1,s1,0.,'''sigma = '',f4.1')
		call number(1.7,3.7,0.1,r1,0.,'''r = '',f6.2')
		call number(2.7,3.7,0.1,b1,0.,'''b = '',f10.7')
		call number(4.1,3.7,0.1,dt1,0.,'''dt = '',f7.4')
		call plot(0.,5.5,-3)
		do 25 i=1,n
			u(i)=t(i)
25		continue
		call text(0.2,4.0,0.2,'z vs Time',0.,9,0)
		call scalex(u,6.,n,1)
		call axctl(0.15,0.04,0.15,0.2,0.2,-1)
		call axisx(0.,0.,'Time (s)',-8,6.,0.,
     c  	u(n+1),u(n+2),u(n+3),4)
		call axisx(0.,3.5,' ',1,6.,0.,u(n+1),u(n+2),u(n+3),20)
		call axctl(0.15,0.04,0.15,0.2,0.2,2)
		call axisx(0.,0.,'z',1,3.5,90.,v(n+1),v(n+2),v(n+3),-5)
		call axisx(6.,0.,' ',-1,3.5,90.
     c		,v(n+1),v(n+2),v(n+3),-21)
		call line(u,v,n,1,lsym,nsym)
		call number(0.5,3.7,0.1,s1,0.,'''sigma = '',f4.1')
		call number(1.7,3.7,0.1,r1,0.,'''r = '',f6.2')
		call number(2.7,3.7,0.1,b1,0.,'''b = '',f10.7')
		call number(4.1,3.7,0.1,dt1,0.,'''dt = '',f7.4')
	else
		do 35 i=1,n
			u(i)=t(i)
			v(i)=y(i)
35		continue
		call text(0.2,3.9,0.2,'y vs Time',0.,9,0)
		call scalex(u,6.,n,1)
		call axctl(0.15,0.04,0.15,0.2,0.2,-1)
		call axisx(0.,0.,'Time (s)',-8,6.,0.,u(n+1),u(n+2),u(n+3),4)
		call axisx(0.,3.5,' ',1,6.,0.,u(n+1),u(n+2),u(n+3),20)
		call scalex(v,3.5,n,1)
		call axctl(0.15,0.04,0.15,0.2,0.2,2)
		call axisx(0.,0.,'y',1,3.5,90.,v(n+1),v(n+2),v(n+3),-5)
	call axisx(6.,0.,' ',-1,3.5,90.,v(n+1),v(n+2),v(n+3),-21)
		call line(u,v,n,1,lsym,nsym)
		call number(0.5,3.7,0.1,s1,0.,'''sigma = '',f4.1')
		call number(1.7,3.7,0.1,r1,0.,'''r = '',f6.2')
		call number(2.7,3.7,0.1,b1,0.,'''b = '',f10.7')
		call number(4.1,3.7,0.1,dt1,0.,'''dt = '',f7.4')
		call plot(0.,5.5,-3)
		do 45 i=1,n
			v(i)=x(i)
45		continue
		call text(0.2,3.9,0.2,'x vs Time',0.,9,0)
		call axctl(0.15,0.04,0.15,0.2,0.2,-1)
		call axisx(0.,0.,'Time (s)',-8,6.,0.,
     c  	u(n+1),u(n+2),u(n+3),4)
		call axisx(0.,3.5,' ',1,6.,0.,u(n+1),u(n+2),u(n+3),20)
		call scalex(v,3.5,n,1)
		call axctl(0.15,0.04,0.15,0.2,0.2,2)
		call axisx(0.,0.,'x',1,3.5,90.,v(n+1),v(n+2),v(n+3),-5)
		call axisx(6.,0.,' ',-1,3.5,90.
     c		,v(n+1),v(n+2),v(n+3),-21)
		call line(u,v,n,1,lsym,nsym)
		call number(0.5,3.7,0.1,s1,0.,'''sigma = '',f4.1')
		call number(1.7,3.7,0.1,r1,0.,'''r = '',f6.2')
		call number(2.7,3.7,0.1,b1,0.,'''b = '',f10.7')
		call number(4.1,3.7,0.1,dt1,0.,'''dt = '',f7.4')
	endif
	call plot(0.,0.,999)
	return
	end