program horseshoe
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
!     Project5 - Physics 355 Chaos and dynamics
!     Chaos in the Lorenz model
!     follows set of 1000 trajectories, writes at the end (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)                   :: dx0                            ! distance between initial trajectories  

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

      integer(i4b)                :: n,nx,ny,nz                     ! index of trajectory


      character(len=9),parameter  :: infile='input.dat'             ! input: parameter values 
      character(len=9),parameter   :: outfile='horse.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
      read(1,*) cdummy, dx0
      close(1)
      print *,'Done reading.'


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

      print *,'Start integration...'
      
! Initialize motion
      DO nx=0,9
      DO ny=0,9
      DO nz=0,9
         n=nz+10*ny+100*nx
         x=x0+nx*dx0
         y=y0+ny*dx0
         z=z0+nz*dx0

         t=0.D0
         do while (t<tmax)
            call rkstep(x,y,z,t)  
            t=t+dt
         enddo
         write(3,'(I6,128(1x,e12.6))') n,x,y,z


      ENDDO
      ENDDO
      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