program project2
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
!     Project2 - Physics 5413 - Chaos and Dynamics
!     Biochemical Switch
!     
!     bioswitch2:    9 Feb 2022
!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Parameters
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

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

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

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

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

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

      real(r8b)                   :: smin,smax,s, ds                   ! 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 *,'Finds stable fixed and follows it with increasing s' 
     print *,'Degradation r'
     read*, r
     print *,'minimum stimulus concentration s_min'
     read *, smin
     print *,'maximum stimulus concentration s_max'
     read *, smax
	 print *,'stimulus step ds'
	 read *,ds
     print *,'Initial value of x'
     read *,g0 
      
     open(3,file=datfile,status='replace')
     rewind(3)
     write(3,*) '# Project 2 - Bioswitch'
     write(3,*) '# x0= ', g0
	 write(3,*) '# s_min =', smin
     write(3,*) '# s_max= ', smax
     write(3,*) '# tmax= ',tmax 	 
     write(3,*) '# =================================================================='
     write(3,*) '# r  s  x*'
     close(3)

     open(3,file=datfile,position='append') 
      
! stimulus loop
    s=smin
    do while (s.le.smax)

! Initialize motion
     g(0)=g0
     gtemp=g0
     t=tmin
     do it=1, Nt
        call rkstep(gtemp,t)
        g(it)=gtemp        
        t=t+dt
     enddo

     write(3,'(128(1x,e12.6))') r, s, g(Nt)
     s=s+ds
	 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