c This program solves for the population in the excited state
c   for the simple form of the Quantum Zeno effect
c
c pz = rho-tilde_z
c
      program Quantum_Zeno_Simp
      implicit real*8 (a-h,o-z)
      parameter (ntyp = 11)
      integer nmeas(0:ntyp)
      real*8 tstp(0:ntyp),pz(0:ntyp),tstr(0:ntyp),p2(0:ntyp)

c name the output files
      open(unit=20,file='quant_zeno_simp.txt')
      open(unit=21,file='quant_zeno_simpf.txt')
      
c define pi
      pi = 2.d0*asin(1.d0)

c define Omega
      rabomeg = 2.d0*pi/1.0d-6
c define dt
      dt = 0.25d0*pi/rabomeg/(2**ntyp)
c set tfinal
      tfin = pi/rabomeg
c set the number of time steps
      numtim = tfin/dt
c initialize the Zeno probabilities and times
c pz = probability in initial state
c tstp = the time between measurements
c tstr = the time since last measurement
c nmeas = number of measurements
       do 1 k = 0, ntyp
       pz(k) = 0.5d0
       tstp(k) = 2.d0*tfin/(dfloat((2**(k+1))))
       tstr(k) = 0.0
       nmeas(k) = 2**k
 1     continue

c the 10 loop is the time steps
      do 10 j = 0, numtim
c define the time for step j
      tim =  dfloat(j)*dt
c compute the probability at time t to be in the excited state for the different measurements
       do 12 k = 0, ntyp
       p2(k) = 0.5d0-pz(k)*cos(rabomeg*(tim-tstr(k)))
 12     continue
c write out the probabilities
      write(20,900)tim*1.e6,(p2(k),k=0,ntyp)
c test whether a measurement occurs
       do 11 k = 0, ntyp
        if(tim .ge. (tstr(k)+tstp(k))) then
        pz(k) = pz(k)*cos(rabomeg*(tim-tstr(k)))
        tstr(k) = tim
        end if
 11    continue
 10   continue
c end of time loop

c output the final probabilities
      do 13 k = 0,ntyp
      j = 2*k+2
      l = k+1
      frac = pi/nmeas(k)
c the exact value
      dum = 0.5*dcos(frac)**(nmeas(k))
      write(21,900)0.5d0*(2**(k+1)),p2(k),0.5d0-dum
 13   continue
      
c define the format of output
 900  format(18(1pe15.7))
      stop 'end of program'
      end
