program project2
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
!     Project2 - Physics 5413 - Chaos and Dynamics
!     Biochemical Switch
!     
!     project2:    28 Jan 2008    first version
!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Parameters
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

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

       
      real(r8b), parameter        :: tmin=0.D0,tmax=100.D0           ! start and stop time 
      integer(i4b), parameter     :: Nt=1000                         ! number of time steps
      real(r8b), parameter        :: dt=(tmax-tmin)/Nt

      character(len=12),parameter :: datfile='trajecto.dat'           ! output file 
     

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Variables
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

      real(r8b)                   :: t                              ! time
      real(r8b)                   :: g(0:Nt)                        ! gene product
      real(r8b)                   :: gtemp

      real(r8b)                   :: s=0.D0                              ! stimulus
      real(r8b)                   :: r=1.D0                              ! degradation
      real(r8b)                   :: g0=0.5D0                            ! initial value of g
      

      integer(i4b)                :: it



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

     print *,'Project 2 - biochemical switch'
	 print *,'Computes individual trajectory'
     print *,'Degradation r'
     read*, r
     print *,'Stimulus concentration s'
     read *, s
     print *,'Initial value of x'
     read *,g0 
      
     open(3,file=datfile,status='replace')
     rewind(3)
     write(3,*) '# Project 2 - Bioswitch'
     write(3,*) '# s=   ',s,' r= ',r
     write(3,*) '# tmin= ',tmin, ' tmax= ',tmax 
     write(3,*) '# =================================================================='
     write(3,*) '# it  t   x'
     close(3)

! Initialize motion
     g(0)=g0
     gtemp=g0
     t=tmin

     do it=1, Nt
        call rkstep(gtemp,t)
        g(it)=gtemp        
        t=t+dt
     enddo


     open(3,file=datfile,position='append') 
     do it=0,Nt
        t=tmin+it*dt
        write(3,'(I6,1x,F8.4,128(1x,e12.6))') it,t,g(it)
     enddo
     close(3)



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

      subroutine rkstep(z1,tt)
      real(r8b)    :: z1,tt
      real(r8b)    :: k11,k12,k13,k14


      k11=dt*g1(z1,tt)

      k12=dt*g1(z1+0.5*k11,tt+0.5*dt)

      k13=dt*g1(z1+0.5*k12,tt+0.5*dt)

      k14=dt*g1(z1+k13,tt+dt)

      z1=z1+(k11+2*k12+2*k13+k14)/6.D0

      return
      end subroutine rkstep

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
      function g1(z1,tt)
      real(r8b)    :: z1,z2,tt,g1
      g1=s-r*z1+(z1*z1)/(1+z1*z1)  
      return
      end function g1



      end