c This program solves the time dependent equations for
c   a 2 state system density matrix.
c
c The calculations have the amplitude for the E-field be constant
c
c This program uses the leapfrog algorithm 2nd order Runge-Kutta algorithm
c F(t+dt/2) = F(t) + (dt/2) F'(t)
c F(t+dt) = F(t) + dt F'(t+dt/2)
c
      program densmat_2state_constE
      implicit real*8 (a-h,o-z)
      real*8 rhotil(3), drhotil(3), rhowrk(3)
      complex *16 icom, c0(2), crotw(2)
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='densmat_const_2.0_1.0_2.0.txt')
      open(unit=21,file='blochsp_const_2.0_1.0_2.0.txt')
c the electric field is assumed proportional to exp(-t^2/tau^2) cos(omega t)
      tau = 0.1d-6
c define Omega
      rabomeg = 2.0*2.d0*pi/1.0d-6
c define the detuning
      detun = 1.0/tau
c define the general Rabi frequency
      genrab = sqrt(detun*detun + rabomeg*rabomeg)
c define the decay rate in units of 10^6 s^-1
      gamma = 2.0/1.0d-6
      
c set the time step (you should check convergence with respect to dt)
      dt = 0.0001d0*2.d0*pi/genrab
c define the starting (tstr) and final time (tfin)
      tstr = 0.d0
      tfin = 15.d0*2.d0*pi/sqrt(detun*detun + rabomeg*rabomeg)
c define the number of steps in the numerical integration
      numtim = nint(((tfin-tstr)/dt))
c initialize the coefficients in the two state calculation
      c0(1) = dcmplx(1.d0,0.d0)
      c0(2) = dcmplx(0.d0,0.d0)
c initialize the coefficients in the two state calculation
c    see page 8 of the Chap7.pdf for definition in terms of coefficients
c    rhotil(1) = rhotilde_x, rhotil(2) = rhotilde_y, rhotil(3) = rhotilde_z
      rhotil(1)=0.5*(c0(1)*conjg(c0(2))+conjg(c0(1))*c0(2))
      rhotil(2)=0.5*((c0(1)*conjg(c0(2))-conjg(c0(1))*c0(2))*icom)
      rhotil(3)=0.5*(abs(c0(1))**2 - abs(c0(2))**2)
      
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 j = 0, numtim
c define the time for step j
      tim = tstr + dfloat(j)*dt
c compute the drhotilde/dt at time tim
      drhotil(1)=-detun*rhotil(2)-0.5*gamma*rhotil(1)
      drhotil(2)=detun*rhotil(1)-rabomeg*rhotil(3)-0.5*gamma*rhotil(2)
      drhotil(3)=rabomeg*rhotil(2)-gamma*(rhotil(3)-0.5)
c do the 1/2 time step and store in a work array
      rhowrk(1) = rhotil(1) + 0.5*dt*drhotil(1)
      rhowrk(2) = rhotil(2) + 0.5*dt*drhotil(2)
      rhowrk(3) = rhotil(3) + 0.5*dt*drhotil(3)
c compute the drhotilde/dt at time tim+dt/2
      drhotil(1)=-detun*rhowrk(2)-0.5*gamma*rhowrk(1)
      drhotil(2)=detun*rhowrk(1)-rabomeg*rhowrk(3)-0.5*gamma*rhowrk(2)
      drhotil(3)=rabomeg*rhowrk(2)-gamma*(rhowrk(3)-0.5)
c do the full time step
      rhotil(1) = rhotil(1) + dt*drhotil(1)
      rhotil(2) = rhotil(2) + dt*drhotil(2)
      rhotil(3) = rhotil(3) + dt*drhotil(3)

c rhotil is now at time given by tim+dt
      tim = tstr + dfloat(j+1)*dt

c every nmod steps output the populations
      if(mod(j,nmod) .eq. 0) then
c write the time in microsecs, P_1, P_2, P_1+P_2 to the screen
      p1 =0.5+rhotil(3)
      p2 = 0.5-rhotil(3)
      write(6,900)tim*1.e6,p1,p2,p1+p2
c crotw is the rotating wave approximation to the c
      dmc = dcos(0.5d0*genrab*tim)
      dms = dsin(0.5d0*genrab*tim)
      crotw(1) = c0(1)*dcmplx(dmc,-dms*detun/genrab)-
     1           c0(2)*icom*rabomeg*dms/genrab
      crotw(2) = c0(2)*dcmplx(dmc, dms*detun/genrab)-
     1           c0(1)*icom*rabomeg*dms/genrab
c write the timein microsecs, P_1, P_2, P_1+P_2, P_1RW, P_2RW to the output file
      write(20,900)tim*1.e6,p1,p2,p1+p2,abs(crotw(1))**2,
     1  abs(crotw(2))**2
      write(21,900)tim*1.e6,rhotil(1),rhotil(2),rhotil(3)
      end if
 10   continue
      close(20)
      close(21)
c end of time loop
      
c define the format of output
 900  format(8(1pe15.7))
      stop 'end of program'
      end
