c  This program reads the results of modes_efield.f and computes a
c   rate.
        Implicit None
        Integer i,j,ik,nk,n,junk,it,nt,kn
        Real*8 k,lam,t2,t,B2,pi,b
        Real*8 sum,delt,rate,Ikz,ktest,p,rateex,ratedif
        Real*8 int,Itest,inttest,Itexact,diftest
        Real*8 Ar,Ai,Br,Bi,testw,testB2,difB2
        Real*8 jz,jz1,jz2,intjz1,intjz2,w0,Ijz1,Ijz2,tmatch
        Real*8 cosw0,sinw0,tprint
        Real*16 kk,ki,delk
        Parameter(kn=1000000)
        Dimension k(kn),p(kn),Ar(kn),Br(kn),Ai(kn),Bi(kn),B2(kn)
        Open(unit = 10, file = 'rate_efield.dat')
        OPEN(UNIT = 11, FILE = 'modes_efield.dat')
        OPEN(UNIT = 12, FILE = 'B2.1')
        OPEN(UNIT = 13, FILE = 'rate.out')
c        Open(unit = 14, file = 'rate.tst')
        Open(unit=15,file='jz.integrand')
        Open(unit=16,file='jz.out')
        Open(unit=17,file='jz.integrand1')
        READ(10,*) JUNK
c  tprint is the time the second integrand is printed out
        Read(10,*) delt, nt, tprint
        READ(11,*) JUNK
        READ(11,*) ki, delk, nk
        READ(11,*) JUNK
        Read(11,*) lam,t2,b
        pi = 2.d0*dasin(1.d0)
        Ikz = 0.d0
        Ijz1 = 0.d0
        Ijz2 = 0.d0
        Itest = 0.d0
        Do ik = 1,nk
          If(ik.eq.1) then
            kk = ki
            k(ik) = kk
          Else
            kk = kk+delk 
            k(ik) = kk
          Endif
c  Note that tmatch is the same for all values of k
777       Read(12,*) t2,tmatch,ktest,Ar(ik),Ai(ik),Br(ik),Bi(ik),B2(ik)
          If(dabs(k(ik)-ktest).gt.1.d-10) then
c            Print*,' Wrong value of k.  k = ',k(ik),' ktest = ',ktest
c            Stop
             Go to 777
          Endif
          If(ik.eq.1.or.ik.eq.nk) then 
            p(ik) = delk/3.d0
            n = 1
          ElseIf(n.eq.1) then 
            p(ik) = 4.d0*delk/3.d0
            n = 0
          Else
            p(ik) = 2.d0*delk/3.d0
            n = 1
          Endif 
c  This is to test the Bogolubov coefficients
          testw = Ai(ik)**2 + Ar(ik)**2 - Bi(ik)**2 - Br(ik)**2
          testB2 = Br(ik)**2+Bi(ik)**2
          difB2 = (B2(ik)-testB2)/B2(ik)
c          Write(14,19) k, testw,B2,testB2,difB2
          int = dlog(1.d0+B2(ik))
c  I'll start just doing intzj at t = 0 so below is only for that case
          w0 = dsqrt((k(ik)+t2/2.d0)**2+lam)
          intjz1 = ((k(ik)+t2/2.d0)/w0)*B2(ik)
c          intjz2 = ((k(ik)+t2/2.d0)/w0)*2.d0*(Ar(ik)*Br(ik)+Ai(ik)*Bi(ik))
c These are the values at the time t = tmatch.  See below for other
c  times
          cosw0 = 1.d0
          sinw0 = 0.d0
          intjz2 = ((k(ik)+t2/2.d0)/w0)*2.d0*((Ar(ik)*Br(ik)+Ai(ik)*Bi(ik))*cosw0
     -              - (Ar(ik)*Bi(ik)-Ai(ik)*Br(ik))*sinw0)
          Write(15,19) k(ik),B2(ik),int,intjz1,intjz2
c  This is to test the integration routine
          inttest = k(ik)**2
          Ikz = Ikz+int*p(ik)
          Ijz1 = Ijz1+intjz1*p(ik)
          Ijz2 = Ijz2+intjz2*p(ik)
          Itest = Itest+inttest*p(ik)
        Enddo
        Itexact = (1.d0/3.d0)*(k(nk)**3-ki**3)
        diftest = (Itexact-Itest)/Itexact
c  Divide by 2*t2 to get a rate
        rate = Ikz/(2.d0*pi*t2)
        rateex = dlog(1.d0+dexp(-2.d0*pi*lam))/(2.d0*pi)
        ratedif = (rateex-rate)/rateex
        jz1 = Ijz1/(2.d0*pi)
        jz2 = Ijz2/(2.d0*pi)
        jz = jz1+jz2
        Write(13,*) 'Test integral results'
        Write(13,19) ki,k(nk),Itexact,Itest,diftest
        Write(13,*) 'Rate results with scaled variables b = ',b,' and t2 = ',t2
        Write(13,19) Ikz,rate,rateex,ratedif
        Write(13,*) ' B2 contribution to the current = ', jz1
        Write(13,*) 'Oscillating contribution to the curreent = ', jz2
        Write(13,*) 'Total current at the matching time ',jz 
c  Now start the time do loop
        t = tmatch-delt
        Do it=1,nt
          t = t+delt
          Ijz2 = 0.d0
          Do ik = 1,nk
            w0 = dsqrt((k(ik)+t2/2.d0)**2+lam)
            cosw0 = qcos(2.q0*w0*(t-tmatch))
            sinw0 = qsin(2.q0*w0*(t-tmatch))
            intjz2 = ((k(ik)+t2/2.d0)/w0)*2.d0*((Ar(ik)*Br(ik)+Ai(ik)*Bi(ik))*cosw0
     -                - (Ar(ik)*Bi(ik)-Ai(ik)*Br(ik))*sinw0)
            If(dabs(t-tprint).lt.1.d-5) then
              Write(17,19) k(ik),intjz2
            Endif
            Ijz2 = Ijz2+intjz2*p(ik)
          Enddo
          jz2 = Ijz2/(2.d0*pi)
          jz = jz1+jz2
          Write(16,19) t,jz1,jz2,jz
        Enddo
  19    Format(1x,9(d14.7,1x))
        Stop
        End