c  This program compares the results of two output files and computes
c  delphi, delphi_c, and delpsi2 which can then easily be plotted.
	IMPLICIT None
	Integer nm,n,i,itt,itest,nk1,nk2,nk3,nt,nktot,it,it1,it2,junk,iprint
	Real*8 c,c1,c2,c3,delk1,delk2,delk3,k,t,t1,t2,delt
	Real*8 pi,corr,mphi,xi,lamphi,g,mu
	Real*8 fr,frp,fi,fip,tt,ttest,phi,phip,phi2p,psi2,delphi
	Real*8 B,state,ratio,delphip,delphi2p,Omega,Wm1,psi2term
	Real*8 phitwo,phitwop,phitwo2p,psi2two
	Real*8 epsilon,delpsi2,delpsi2exact,diff,ratio1
	Real*8 ratio2,ratio3,I5tot,psi2an,psi2twoan
	Real*8 delphic,g2phidelpsi2,g2phidelpsi2an
	PARAMETER(nm=16000,itt=21000)
	DIMENSION c1(nm),c2(nm),c3(nm),k(nm)
	Common/in1/c1,c2,c3,k
	Common/in2/nk1,nk2,nk3
	OPEN(14,FILE='rfs.dat')
	OPEN(15,FILE='rfs.out')
	Open(16,file='rfs.2')
	OPEN(19,FILE='compare_delphi.out')
c reading the variable junk is just a trick to allow for the labeling of the 
c data file.
c  The first set of variables are not used in this program
	READ(14,*) junk
	READ(14,*) mphi,xi,lamphi,g
	READ(14,*) junk
	READ(14,*) DELk1,nk1,delk2,nk2,delk3,nk3,DELT,NT
	pi = 2.d0*dasin(1.d0)
	Print *,'pi = ',pi
        Print *,'nt = ',nt
!$$$This is the lower limit cutoff for the integrals that go like 1/k in
!$$$the analytic approximation.  

!$$$c  We'll set the parameter mu in the integrals of the analytic approximation
!$$$c   to be equal to 1.
	mu = 1.d0
c  Now read the junk in the rfs.out file so it is ready to go
	Read(15,*) junk
	Read(15,*) junk
	Read(15,*) junk
	Read(15,*) junk
	Read(15,*) junk
	Read(15,*) junk
	Read(15,*) junk
	Read(15,*) junk
	Read(15,*) junk
	Read(16,*) junk
	Read(16,*) junk
	Read(16,*) junk
	Read(16,*) junk
	Read(16,*) junk
	Read(16,*) junk
	Read(16,*) junk
	Read(16,*) junk
	Read(16,*) junk
	Do it=1,nt+1
	  Read(15,*) t,phi,phip,phi2p,psi2,psi2an
          If(it.eq.1.and.t.gt.1.d-12) then
            Print*,' it = ',it,'t = ',t
            Stop
          Endif
          Read(16,*) t1,phitwo,phitwop,phitwo2p,psi2two,psi2twoan
          If(dabs(t-t1).gt.1.d-12) then
            Print*,'it = ',it,' t = ',t,' t1 = ',t1
            Stop
          Endif
            If(it.eq.1) then
              c = (phitwo-phi)/phi
            Endif
            delphi = phitwo - phi
            delphic = delphi - c*phi
c  I want for delpsi2 to have the source term for delphic
c  Of course this is for the exact expresion.  The one in the linear
c  response equation is different slightly.
            delpsi2 = g**2*(psi2two - psi2)*phi
            delpsi2 = psi2two - psi2
            g2phidelpsi2 = g**2*(psi2two - psi2)*phi
	    g2phidelpsi2an = g**2*(psi2twoan-psi2an)*phi
	  Write(19,27) t,phi,delphi,delphic,psi2,delpsi2,g2phidelpsi2,g2phidelpsi2an
	Enddo
  27    format(10(e21.14,1x))
	STOP
	END