Program Lorenz1
*	Compututation of Lyapanov exponents for Lorenz equations
*	using simple Euler integration
*	Program based on Wolf et al., Physica 16D, 285 (1985)
*	Method discussed in Tobochnik & Gould, Computer in Physics,
*	Nov/Dec 1989.
*	December 22, 1989
	real*8 x,y,z,dx(3),dy(3),dz(3)
	real*8 scum(3)
	real*8 t,t0,dt
*	parameters of Lorenz equations
	real*8 a,b,r
*	output variables
	integer ntot,nout,n0
	common /xyz/ x,y,z,dx,dy,dz
	common /time/ t,dt,t0
	common /output/ ntot,nout,n0
	common /spect/scum
	common/params/a,b,r
	call initial
*	iterate Lorentz eqs., linearized equations, and accumulate data
	call update
	stop
	end

	subroutine initial
	real*8 x,y,z,dx(3),dy(3),dz(3)
	real*8 scum(3)
	real*8 t,t0,dt
	real*8 a,b,r	
	integer ntot,nout,n0
	common /xyz/ x,y,z,dx,dy,dz
	common /time/ t,dt,t0
	common /output/ ntot,nout,n0
	common /spect/scum
	common/params/a,b,r
	a = 16.0d0
	b = 4.0d0
	r = 45.92d0
	write(*,*) 'a, b, r, =',a,b,r
*	initial values of x, y, z
	x = 10.0d0
	y = 1.0d0
	z = 0.0d0
*	initialize unit difference vectors
	do 20 i = 1,3
	   dx(i) = 0.0d0
	   dy(i) = 0.0d0
	   dz(i) = 0.0d0
	   scum(i) = 0.0d0
20	continue
*	initial vectors (1,0,0), (0,1,0), (0,0,1)
	dx(1) = 1.0d0
	dy(2) = 1.0d0
	dz(3) = 1.0d0
*	try dt = 0.001
	write(*,*) 'time step dt = '
	read(*,*) dt
*	try n0 = 10,000 iterations
	write(*,*) '# of iterations for transients = ?'
	read(*,*) n0
	t0 = float(n0)*dt
*	collect data after t0
*	try sample time of 100,000 iterations
	write(*,*) '# of iterations for spectrum = ?'
	read(*,*) ntot
	ntot = n0 + ntot
	write(*,*) 'number of steps between output = ?'
	read(*,*) nout
	write(*,*) '   t       lambda(1)       lambda(2)    lambda(3)'
	t = 0.0d0
	end

	subroutine update
	real*8 x,y,z,dx(3),dy(3),dz(3)
	real*8 scum(3)
	real*8 t,t0,dt
	real*8 a,b,r	
	integer ntot,nout,n0
	common /xyz/ x,y,z,dx,dy,dz
	common /time/ t,dt,t0
	common /output/ ntot,nout,n0
	common /spect/scum
	common/params/a,b,r	
*	local variables
	real*8 dot12,dot13,dot23,norm(3)
	integer i,j
	do 100 i = 1,ntot
	   call euler
	   t = t + dt
*	   normalize first difference vector
	   norm(1) = sqrt( dx(1)*dx(1) + dy(1)*dy(1) + dz(1)*dz(1) )
	   dx(1) = dx(1)/norm(1)
	   dy(1) = dy(1)/norm(1)
	   dz(1) = dz(1)/norm(1)
*	   apply Gram-Schmidt procedure
	   dot12 = dx(1)*dx(2) + dy(1)*dy(2) + dz(1)*dz(2)
	   dx(2) = dx(2) - dot12*dx(1)
	   dy(2) = dy(2) - dot12*dy(1)
	   dz(2) = dz(2) - dot12*dz(1)
	   norm(2) = sqrt( dx(2)*dx(2) + dy(2)*dy(2) + dz(2)*dz(2) )
	   dx(2) = dx(2)/norm(2)
	   dy(2) = dy(2)/norm(2)
	   dz(2) = dz(2)/norm(2)
	   dot13 = dx(1)*dx(3) + dy(1)*dy(3) + dz(1)*dz(3)
	   dot23 = dx(2)*dx(3) + dy(2)*dy(3) + dz(2)*dz(3)
	   dx(3) = dx(3) - dot13*dx(1) - dot23*dx(2)
	   dy(3) = dy(3) - dot13*dy(1) - dot23*dy(2)
	   dz(3) = dz(3) - dot13*dz(1) - dot23*dz(2)
	   norm(3) = sqrt( dx(3)*dx(3) + dy(3)*dy(3) + dz(3)*dz(3) )
	   dx(3) = dx(3)/norm(3)
	   dy(3) = dy(3)/norm(3)
	   dz(3) = dz(3)/norm(3)
*	   accumulate data after transient time t0
	   if (i .gt. n0) then
	      do 10 j = 1,3
	         scum(j) = scum(j) + log(norm(j))/log(2.0)
10	      continue
	      if (mod(i,nout) .eq. 0) then
		 write(*,'(f10.2,3(5x,f10.5))') t,(scum(j)/(t-t0),j = 1,3)
	      endif
	   end if
100	continue
	end
	
	subroutine euler
*	solve diff. eqs. using Euler method
	real*8 x,y,z,dx(3),dy(3),dz(3)
	real*8 scum(3)
	real*8 t,t0,dt
	real*8 a,b,r	
	integer ntot,nout,n0
	common /xyz/ x,y,z,dx,dy,dz
	common /time/ t,dt,t0
	common /output/ ntot,nout,n0
	common /spect/scum
	common/params/a,b,r	
*	local variables
	real*8 xp,yp,zp,dxp(3),dyp(3),dzp(3)
	integer i
*	update Lorenz equations
	xp = x + dt*a*( y - x )
	yp = y + dt*( x*(r - z) - y )
	zp = z + dt*( x*y - b*z )
*	update linearized equations for each vector
	do 10 i = 1,3
	   dxp(i) = dx(i) + a*( dy(i) - dx(i) )*dt
	   dyp(i) = dy(i) + ( dx(i)*(r - z) - x*dz(i) - dy(i) )*dt
	   dzp(i) = dz(i) + ( dx(i)*y + x*dy(i) - b*dz(i) )*dt
10	continue
	x = xp
	y = yp
	z = zp
	do 20 i = 1,3
	   dx(i) = dxp(i)
	   dy(i) = dyp(i)
	   dz(i) = dzp(i)
20	continue
	end