program lorenz
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
!     Project5 - Physics 355 Chaos and dynamics
!     Chaos in the Lorenz model
!     Writes one full trajectory for times between 0 and tmax
!     
!     4 Mar 2004:    first version
!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Parameters and variables
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

      implicit none
      integer,parameter           :: r8b= kind(1.D200)
      integer,parameter           :: i4b= kind(2147483647)
      real(r8b), parameter        :: pi=3.14159265358979323D0

      real(r8b)                   :: p,b,r                          ! parameters in Lorenz equation
      real(r8b)                   :: x0,y0,z0                       ! initial values
      real(r8b)                   :: tmax,dt                        ! integration time, time step

      real(r8b)                   :: x,y,z                          ! dynamical variables
      real(r8b)                   :: t                              ! time


      character(len=9),parameter  :: infile='input.dat'             ! input: parameter values 
      character(len=10),parameter :: outfile='trajec.dat'           ! output file for trajectories
      character(len=128)          :: cdummy                         ! dummy for reading parameter file 
     

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Start main program
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

      print *,'Reading input data from input.dat...'
      open(1,file=infile)
      read(1,*) cdummy
      read(1,*) cdummy, p              
      read(1,*) cdummy, b
      read(1,*) cdummy, r
      read(1,*) cdummy, x0              
      read(1,*) cdummy, y0
      read(1,*) cdummy, z0
      read(1,*) cdummy, tmax
      read(1,*) cdummy, dt
      close(1)
      print *,'Done reading.'

      open(3,file=outfile,status='replace')
      write(3,*) '# project 5 - chaos in the Lorenz model'
      write(3,'(a,4(1x,e12.6))') '# p,b,r = ', p,b,r
      write(3,*) '# =================================================================='
      write(3,*) '# t   x   y    z'     


! Initialize motion
      x=x0
      y=y0
      z=z0
      t=0.D0

      print *,'Start integration...'
      do while (t<tmax)
         call rkstep(x,y,z,t)  
         write(3,'(128(1x,e12.6))') t,x,y,z
         t=t+dt
      enddo
      print *,'Done.'

      close(3)




!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Now the internal subroutines
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
      contains

      subroutine rkstep(z1,z2,z3,tt)
      real(r8b)    :: z1,z2,z3,tt
      real(r8b)    :: k11,k12,k13,k14
      real(r8b)    :: k21,k22,k23,k24
      real(r8b)    :: k31,k32,k33,k34

      k11=dt*g1(z1,z2,z3,tt)
      k21=dt*g2(z1,z2,z3,tt)
      k31=dt*g3(z1,z2,z3,tt)

      k12=dt*g1(z1+0.5*k11,z2+0.5*k21,z3+0.5*k31,tt+0.5*dt)
      k22=dt*g2(z1+0.5*k11,z2+0.5*k21,z3+0.5*k31,tt+0.5*dt)
      k32=dt*g3(z1+0.5*k11,z2+0.5*k21,z3+0.5*k31,tt+0.5*dt)

      k13=dt*g1(z1+0.5*k12,z2+0.5*k22,z3+0.5*k32,tt+0.5*dt)
      k23=dt*g2(z1+0.5*k12,z2+0.5*k22,z3+0.5*k32,tt+0.5*dt)
      k33=dt*g3(z1+0.5*k12,z2+0.5*k22,z3+0.5*k32,tt+0.5*dt)

      k14=dt*g1(z1+k13,z2+k23,z3+k33,tt+dt)
      k24=dt*g2(z1+k13,z2+k23,z3+k33,tt+dt)
      k34=dt*g3(z1+k13,z2+k23,z3+k33,tt+dt)

      z1=z1+(k11+2*k12+2*k13+k14)/6.D0
      z2=z2+(k21+2*k22+2*k23+k24)/6.D0
      z3=z3+(k31+2*k32+2*k33+k34)/6.D0

      return
      end subroutine rkstep

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
      function g1(z1,z2,z3,tt)
      real(r8b)    :: z1,z2,z3,tt,g1
      g1=p*(z2-z1)  
      return
      end function g1


!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
      function g2(z1,z2,z3,tt)
      real(r8b)    :: z1,z2,z3,tt,g2
      g2=-z1*z3 + r*z1 - z2
      return
      end function g2

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
      function g3(z1,z2,z3,tt)
      real(r8b)    :: z1,z2,z3,tt,g3
      g3=z1*z2 - b*z3
      return
      end function g3

      end