c  This program solves the equations for phi and delta phi using an analytic
c  approximation for <phi^2> - actually two such approximations one with a
c  positive and one with a negative sign.
c  

	IMPLICIT NONE
        Integer it,nt,JUNK, nbad,nok
	Real*8 t,ti,Y,phiplus,phiminus,phi0,delphiplus,delphiminus,delphiexplus
        Real*8 delphiexminus,errphiclass
        Real*8 delphi0,delphicplus,delphicminus,delphicexplus,delphicexminus 
	Real*8 eps,h1,hmin
	Real*8 delphiex,delphicex,delt
        Real*8 phiclass,phi1plus,phi1minus,difplusdelphi,difplusdelphic
        Real*8 difminusdelphi,difminusdelphic
	DIMENSION Y(14)
	External DERIVS, BSSTEP

	OPEN(UNIT = 11, FILE = 'delphian.dat')
        Open(unit=12,file='delphiplus.out')
        Open(unit=13,file='delphiminus.out')
	READ(11,*) JUNK
	READ(11,*) delt,nt,phi0,delphi0
	READ(11,*) JUNK
	Read(11,*) eps,h1,hmin
c  We will start at t = 0
c  Note that the equation for delphian is linear so it really doesn't matter
c   what the initial value is.
c
c  I will set it up so that the first thing I do is integrate the harmonic
c   oscillator as a check:  so that is Y(1) and Y(2).  Then I'll do the
c   equation for phi with 1+phi^2  and then delphi.  Then I'll do the equation
c   for phi with 1-phi^2 and then delphi.
	   Y(1)= 1.d0
	   Y(2)= 0.d0
	   Y(3) = phi0
	   Y(4) = 0.d0
	   Y(5)= delphi0
	   Y(6)= 0.d0
	   Y(7) = phi0
	   Y(8) = 0.d0
	   Y(9)= delphi0
	   Y(10)= 0.d0
c  Now do the exact ones for a larger phi0
           Y(11) = phi0+delphi0
           Y(12) = 0.d0
           Y(13) = phi0+delphi0
           Y(14) = 0.d0
           phiclass = Y(1)
           phiplus = Y(3)
           delphiplus = Y(5)
           phiminus = Y(7)
           delphiplus = Y(9)
           errphiclass = dcos(t)-Y(1)
           phi1plus = Y(11)
           phi1minus = Y(13)
           delphicplus = delphiplus-(delphi0/phi0)*phiplus
           delphicminus = delphiminus-(delphi0/phi0)*phiminus
           delphiexplus = phi1plus - phiplus
           delphiexminus = phi1minus - phiminus
           delphicexplus = delphiexplus-(delphi0/phi0)*phiplus
           difplusdelphi = delphiexplus-delphiplus
           difplusdelphic = delphicexplus-delphicplus
           difminusdelphi = delphiexminus-delphiminus
           difminusdelphic = delphicexminus-delphicminus
           t = 0.d0
	   WRITE(12,35) t,errphiclass,phiplus,difplusdelphi,delphiplus,
     1           difplusdelphic,delphicplus 
	   WRITE(13,35) t,errphiclass,phiminus,difminusdelphi,delphiminus,
     1           difminusdelphic,delphicminus 
c I set t = 0 above
	   DO IT = 1,nt
             ti = t
	     t = t+delt
             CALL ODEINT(Y,14,ti,t,EPS,H1,HMIN,NOK,NBAD,DERIVS,BSSTEP)
             phiclass = Y(1)
             phiplus = Y(3)
             delphiplus = Y(5)
             phiminus = Y(7)
             delphiminus = Y(9)
             errphiclass = dcos(t)-Y(1)
             phi1plus = Y(11)
             phi1minus = Y(13)
             delphicplus = delphiplus-(delphi0/phi0)*phiplus
             delphicminus = delphiminus-(delphi0/phi0)*phiminus
             delphiexplus = phi1plus - phiplus
             delphiexminus = phi1minus - phiminus
             delphicexplus = delphiexplus-(delphi0/phi0)*phiplus
             delphicexminus = delphiexminus-(delphi0/phi0)*phiminus
             difplusdelphi = delphiexplus-delphiplus
             difplusdelphic = delphicexplus-delphicplus
             difminusdelphi = delphiexminus-delphiminus
             difminusdelphic = delphicexminus-delphicminus
	     WRITE(12,35) t,errphiclass,phiplus,difplusdelphi,delphiplus,
     1           difplusdelphic,delphicplus 
	     WRITE(13,35) t,errphiclass,phiminus,difminusdelphi,delphiminus,
     1           difminusdelphic,delphicminus 
           Enddo
  35    Format(1x,10(d14.7,1x))
	STOP
	END

!	--------------------------------------------------------------
!	-------------- SUBROUTINES -----------------------------------
!	--------------------------------------------------------------

	SUBROUTINE DERIVS(t,Y,DYDT,NVAR)
        Implicit None 
	INTEGER nvar
	REAL*8 t, Y(nvar), DYDT(nvar) 
        Real*8 phiplus,delphiplus,phiminus,delphiminus,phi1plus,phi1minus
             phiplus = Y(3)
             delphiplus = Y(5)
             phiminus = Y(7)
             delphiminus = Y(9)
             phi1plus = Y(11)
             phi1minus = Y(13)
	DYDT(1) = Y(2)
	DYDT(2) = -Y(1)
	DYDT(3) = Y(4)
	DYDT(4) = -(1.d0+phiplus**2)*phiplus
	DYDT(5) = Y(6)
	DYDT(6) = -(1.d0+3.d0*phiplus**2)*delphiplus
	DYDT(7) = Y(8)
	DYDT(8) = -(1.d0-phiminus**2)*phiminus
	DYDT(9) = Y(10)
	DYDT(10) = -(1.d0-3.d0*phiminus**2)*delphiminus
	DYDT(11) = Y(12)
	DYDT(12) = -(1.d0+phi1plus**2)*phi1plus
	DYDT(13) = Y(14)
	DYDT(14) = -(1.d0-phi1minus**2)*phi1minus
	RETURN
	END