Program Converge04

      IMPLICIT NONE

      INTEGER   i, j, k, l, Nx1, Np1, Np2, Np4, Nmax1, Nmax2, Nmax4, 
     /          Ntime, Nbound, Nmovie, Ncloak, Nuse1, Nuse2, Nuse4,
     /          Npick
      PARAMETER ( Np1 = 200, Np2 = 400, Np4 = 800, Nmax1 = 401, 
     /            Nmax2 = 801, Nmax4 = 1601 )
      REAL*8    DT1, DX1, DT2, DX2, DT4, DX4, thetax, thetat, beta, 
     /          Delta, Time, factor, Work(Nmax4),
     /          Q1(-1:1,-Np4:Np4), Mat1(Nmax4,Nmax4), Vec1(Nmax4),
     /          Q2(-1:1,-Np2:Np2), Mat2(Nmax2,Nmax2), Vec2(Nmax2),
     /          Q4(-1:1,-Np1:Np1), Mat4(Nmax1,Nmax1), Vec4(Nmax1),
     /          Cloak1(-Np4:Np4), Cloak2(-Np2:Np2), Cloak4(-Np1:Np1)
      INTEGER   Or1(Nmax4), Or2(Nmax2), Or4(Nmax1)

      OPEN(7, FILE = 'Converge.dat' , FORM = 'formatted'  )
         READ(7, 200, ERR = 1000 ) DT4, DX4, thetax, thetat, beta, 
     /                             Delta, Time
 200     FORMAT(6(11X,E11.4 / ),11X,E11.4)
         READ(7, 210, ERR = 1000 ) Nbound, Nmovie, Npick, Ncloak
 210     FORMAT(3(11X,I5 / ),11X,I5)
      CLOSE(7)

      Nx1 = NINT(1.0D0/DX4)
      Nuse4 = 8*Nx1 + Nbound
      Nuse2 = 4*Nx1 + Nbound
      Nuse1 = 2*Nx1 + Nbound
      DT2 = DT4 / 2.0D0
      DX2 = DX4 / 2.0D0
      DT1 = DT4 / 4.0D0
      DX1 = DX4 / 4.0D0

      IF ( Nx1 .GT. Np1 ) THEN
         PRINT *, 'The coarsest grid is too fine.'
         STOP
      ENDIF
      Ntime = NINT(Time/DT4)

C Set up your initial conditions

      CALL INITIAL(DT1,DX1,4*Nx1,Q1,Delta,beta,Cloak1)
      CALL INITIAL(DT2,DX2,2*Nx1,Q2,Delta,beta,Cloak2)
      CALL INITIAL(DT4,DX4,Nx1,Q4,Delta,beta,Cloak4)

      CALL MATRIX(DT1,DX1,thetax,thetat,beta,Nbound,Nuse4,Mat1,Work,Or1)
      CALL MATRIX(DT2,DX2,thetax,thetat,beta,Nbound,Nuse2,Mat2,Work,Or2)
      CALL MATRIX(DT4,DX4,thetax,thetat,beta,Nbound,Nuse1,Mat4,Work,Or4)

C Now evolve the solutions to calculate the norm

      OPEN( UNIT = 13, file = 'Factor.dat' )

      DO i = 1, Ntime
         DO k = 2*i-1, 2*i
            DO l = 2*k-1, 2*k
               CALL EVOL(DT1,DX1,thetax,thetat,beta,Nbound,4*Nx1,Q1,
     /                   Nuse4,Mat1,Vec1,Or1,Cloak1,Ncloak)
            ENDDO
            CALL EVOL(DT2,DX2,thetax,thetat,beta,Nbound,2*Nx1,Q2,
     /                Nuse2,Mat2,Vec2,Or2,Cloak2,Ncloak)
         ENDDO
         CALL EVOL(DT4,DX4,thetax,thetat,beta,Nbound,Nx1,Q4,Nuse1,
     /             Mat4,Vec4,Or4,Cloak4,Ncloak)

         CALL FACTORING(i,DX1,4*Nx1,Cloak1,2*Nx1,Cloak2,Nx1,Cloak4,
     /                  factor,Nmovie,Npick)
         PRINT 500, DT4*i, factor
         WRITE(13,500) DT4*i, factor
 500     FORMAT(1X,F9.6,1P,2X,E13.6)
 510     FORMAT(1X,F9.6,1P,2(2X,E13.6))
      ENDDO

      CLOSE(13)
      STOP

 1000 PRINT *, 'Error in reading file'
      STOP
      END

C***********************************************************************

      SUBROUTINE FACTORING(i,DX1,Nx4,Cloak1,Nx2,Cloak2,Nx1,Cloak4,
     /                     factor,Nmovie,Npick)

      IMPLICIT NONE

      INTEGER i, j, Nx4, Nx2, Nx1, Nmovie, Npick
      REAL*8  Cloak1(-Nx4:Nx4), Cloak2(-Nx2:Nx2), Cloak4(-Nx1:Nx1), 
     /        numer, denom, numpt, denpt, factor, trunc, DX1
      CHARACTER*64 fname

      IF ( Nmovie .EQ. 0 ) GOTO 500

      IF ( I .LT. 10 ) THEN
         WRITE(fname,901) I
 901     FORMAT('slice','.',I1) 
      ELSEIF ( I .GT. 9 .AND. I .LT. 100 ) THEN
         WRITE(fname,902) I
 902     FORMAT('slice','.',I2) 
      ELSEIF ( I .GT. 99 .AND. I .LT. 1000 ) THEN
         WRITE(fname,903) I
 903     FORMAT('slice','.',I3) 
      ELSEIF ( I .GT. 999 .AND. I .LT. 10000 ) THEN
         WRITE(fname,904) I
 904     FORMAT('slice','.',I4) 
      ELSE
         PRINT *, 'I is too high.'
         STOP
      ENDIF

      OPEN( UNIT = 99, FILE = fname )
      IF ( Npick .EQ. 1 ) THEN
         DO j = -Nx4, Nx4
            IF ( MOD(j,2) .EQ. 0 ) THEN
               trunc = ( Cloak2(j/2) - Cloak1(j) ) / 4.0D0
               IF ( ABS(trunc) .LT. 1.0D-89 ) trunc = 0.0D0
               WRITE(99,320) j*DX1, trunc
            ENDIF
         ENDDO
      ELSEIF ( Npick .EQ. 2 ) THEN
         DO j = -Nx2, Nx2
            IF ( MOD(j,2) .EQ. 0 ) THEN
               trunc = ( Cloak4(j/2) - Cloak2(j) ) / 4.0D0
               IF ( ABS(trunc) .LT. 1.0D-89 ) trunc = 0.0D0
               WRITE(99,320) j*DX1, trunc
            ENDIF
         ENDDO
      ELSE
         PRINT *, 'Invalid value for Npick.'
         STOP
      ENDIF

 320  FORMAT(1X,F9.6,2X,E17.10)
      CLOSE( UNIT = 99 )

 500  numer = 0.0D0
      denom = 0.0D0
      DO j = -Nx4, Nx4
         IF ( MOD(j,2) .EQ. 0 ) THEN
            denpt = ( Cloak2(j/2) - Cloak1(j) )**2
            denom = denom + denpt
         ENDIF
         IF ( MOD(j,4) .EQ. 0 ) THEN
            numpt = ( Cloak4(j/4) - Cloak2(j/2) )**2
            numer = numer + numpt
         ENDIF
      ENDDO

      numer = SQRT( numer / ( 2.0D0*Nx1 + 1.0D0 ) )
      denom = SQRT( denom / ( 2.0D0*Nx2 + 1.0D0 ) )
      factor = numer / denom

      RETURN
      END