c This program solves the time dependent equations for
c   a 3 state system meant to mimic the quantum Zeno effect.
c State 1 = initial state, 2 = excited state and 3 = state used to probe P_2
c
c The idea is that at a time a laser couples 2 and 3 where 3 has a high radiative
c   decay rate. The laser is on for only a short time to mimic an instantaneous
c   measurement
c
c This program uses the 2nd order Runge-Kutta method to propagate the density
c   matrix
c
      program Quantum_Zeno_dens
      implicit real*8 (a-h,o-z)
      complex *16 rho(3,3), rhoh(3,3), drho(3,3), icom
c define the complex number i
      icom = dcmplx(0.d0,1.d0)
c define pi
      pi = 2.d0*dasin(1.d0)
c the output will be in whatever is defined as file
      open(unit=20,file='quant_zeno_dens4.txt')
c the electric field is assumed proportional to exp(-t^2/tau^2) cos(omega t)
      tau = 0.00125d-6
c define Omega_12
      om12 = 2.d0*pi/1.0d-6
c define Omega_23 prefactor
      om23pref = 3200.0*om12
c define the Decay rate for state 3
      gam3 = 0.5*om23pref
      
c set the time step (you should check convergence with respect to dt)
      dt = 0.00001d0/(100*om12)
c define the starting (tstr) and final time (tfin)
      tstr = 0.d0
      tfin = pi/om12
c define the number of steps in the numerical integration
      numtim = nint(((tfin-tstr)/dt))
c define the time where the measurement occurs
      tme0 = tfin/4.0
      tme = tme0
c initialize the density matrix
       do 1 j = 1, 3
       do 1 k = 1, 3
       rho(j,k)=cmplx(0.0,0.0)
 1     continue
      rho(1,1)=cmplx(1.0,0.0)
      
c nmod is roughly how many steps to skip so the number of output points
c    are the denominator
      nmod = numtim/2000
      
c the 10 loop is the time steps
      do 10 itim = 0, numtim
c define the time for step j
      tim = tstr + dfloat(itim)*dt
      if(tim .gt. (tme + 0.1*tme0)) tme = tme + tme0
c compute the 2-3 coupling at time t
      om23 = 0.0
      dum = ((tim-tme)/tau)**2
       if(dum .lt. 36.0) om23 = om23pref*exp(-dum)
c compute the derivative of the density matrix at t
      drho(1,1) = icom*0.5*om12*(rho(1,2)-rho(2,1))
      drho(1,2) = icom*0.5*om12*(rho(1,1)-rho(2,2))
     1           +icom*0.5*om23*rho(1,3)
      drho(2,1) = conjg(drho(1,2))
      drho(1,3) = -icom*0.5*om12*rho(2,3)+icom*0.5*om23*rho(1,2)
     1           - 0.5*gam3*rho(1,3)
      drho(3,1) = conjg(drho(1,3))
      drho(2,2) = -icom*0.5*om12*(rho(1,2)-rho(2,1))
     1          +icom*0.5*om23*(rho(2,3)-rho(3,2))+gam3*rho(3,3)
      drho(2,3) = -icom*0.5*om12*rho(1,3)
     1          +icom*0.5*om23*(rho(2,2)-rho(3,3))
     1           - 0.5*gam3*rho(2,3)
      drho(3,2) = conjg(drho(2,3))
      drho(3,3) = -icom*0.5*om23*(rho(2,3)-rho(3,2))-gam3*rho(3,3)
c do the 1/2 step
       do 11 j = 1, 3
       do 11 k = 1, 3
       rhoh(k,j) = rho(k,j) + 0.5*dt*drho(k,j)
 11    continue

c compute the derivative of the density matrix at t+dt/2
c compute the 2-3 coupling at time t+dt/2
      om23 = 0.0
      dum = ((tim+0.5*dt-tme)/tau)**2
       if(dum .lt. 36.0) om23 = om23pref*exp(-dum)
      drho(1,1) = icom*0.5*om12*(rhoh(1,2)-rhoh(2,1))
      drho(1,2) = icom*0.5*om12*(rhoh(1,1)-rhoh(2,2))
     1           +icom*0.5*om23*rhoh(1,3)
      drho(2,1) = conjg(drho(1,2))
      drho(1,3) = -icom*0.5*om12*rhoh(2,3)+icom*0.5*om23*rhoh(1,2)
     1           - 0.5*gam3*rhoh(1,3)
      drho(3,1) = conjg(drho(1,3))
      drho(2,2) = -icom*0.5*om12*(rhoh(1,2)-rhoh(2,1))
     1          +icom*0.5*om23*(rhoh(2,3)-rhoh(3,2))+gam3*rhoh(3,3)
      drho(2,3) = -icom*0.5*om12*rhoh(1,3)
     1          +icom*0.5*om23*(rhoh(2,2)-rhoh(3,3))
     1           - 0.5*gam3*rhoh(2,3)
      drho(3,2) = conjg(drho(2,3))
      drho(3,3) = -icom*0.5*om23*(rhoh(2,3)-rhoh(3,2))-gam3*rhoh(3,3)
c do the 1/2 step
       do 12 j = 1, 3
       do 12 k = 1, 3
       rho(k,j) = rho(k,j) + dt*drho(k,j)
 12    continue

c every nmod steps output the populations
      if(mod(itim,nmod) .eq. 0) then
c write the time in microsecs, P_1, P_2, rho_12
      p1 = dreal(rho(1,1))
      p2 = dreal(rho(2,2))
      p3 = dreal(rho(3,3))
       if(p3 .lt. 1.e-18) p3 = 0.0
      write(6,900)(tim+dt)*1.e6,p1,p2,rho(1,2)
c write the timein microsecs, P_1, P_2, P_1+P_2 to the output file
      write(20,900)(tim+dt)*1.e6,p1,p2,p3,rho(1,2),rho(1,3),rho(2,3)
      end if
 10   continue
c end of time loop
      write(6,*)'tau * Gamma_3',tau*gam3
c define the format of output
 900  format(18(1pe15.7))
      stop 'end of program'
      end
