      PROGRAM MAIN
      IMPLICIT NONE
C     .. Parameters ..
      INTEGER NZMAX,NFMAX,IMESH,IMLOC,NEV
      PARAMETER (NZMAX=1000,NFMAX=NZMAX,IMESH=10000,IMLOC=10,NEV=10)
C     ..
C     .. Local Scalars ..
      DOUBLE COMPLEX ANS,ANSBC,ANSLC,ANSRC,ANSTC,
     &               ANSWR1,ANSWR2,DMISS1,DMISS2,EBL,EBR,ELAM,
     &               ELAM1,ELAM2,ETL,ETR,F
      DOUBLE PRECISION A,B,CONST,CONST1,CONST2,RELAM,TOL
      INTEGER I,IFAIL,ILOC,IMATCH,IREVCM,NLOC,NMESH,NZACT1,NZACT2,
     &        NZBC,NZLC,NZRC,NZTC,NEIGS
C     ..
C     .. Local Arrays ..
      DOUBLE COMPLEX FARRY1(0:NFMAX,1:2),FARRY2(0:NFMAX,1:2),
     &       FB(0:NFMAX,1:2,NEV),FBC1(0:NFMAX,1:2),FL(0:NFMAX,1:2,NEV),
     &       FLC1(0:NFMAX,1:2),FR(0:NFMAX,1:2,NEV),FRC1(0:NFMAX,1:2),
     &       FT(0:NFMAX,1:2,NEV),FTC1(0:NFMAX,1:2),ZARRY1(0:NZMAX),
     &       ZARRY2(0:NZMAX),ZB(0:NZMAX,NEV),ZBC1(0:NZMAX),
     &       ZL(0:NZMAX,NEV),ZLC1(0:NZMAX),ZR(0:NZMAX,NEV),
     &       ZRC1(0:NZMAX),ZT(0:NZMAX,NEV),ZTC1(0:NZMAX),
     &       FBC2(0:NFMAX,1:2),FRC2(0:NFMAX,1:2),FTC2(0:NFMAX,1:2),
     &       FLC2(0:NFMAX,1:2),ZBC2(0:NZMAX),ZLC2(0:NZMAX),
     &       ZRC2(0:NZMAX),ZTC2(0:NZMAX)
      DOUBLE COMPLEX ANSL(NEV),ANSR(NEV),ANST(NEV),ANSB(NEV),EIGVAL(NEV)
      INTEGER NZL(NEV),NZR(NEV),NZT(NEV),NZB(NEV),IEIGS(NEV)
      DOUBLE PRECISION XLOC(0:IMLOC),XMESH(0:IMESH)
C     ..
C     .. Scalars in Common ..
      DOUBLE PRECISION BETA
      INTEGER IPROB
C     ..
C     .. External Subroutines ..
      EXTERNAL BCHILD,BXSIDE,LCHILD,LRCHLD,MKBOX,RCHILD,RFNBOX,SL4BCS,
     &         SL4COF,SLMISS,SLSMBC,SLSMPC,TCHILD,TBCHLD,ALLEVS,
     &         GETMSH
C     ..
C     .. Intrinsic Functions ..
      INTRINSIC ABS,ATAN,CMPLX,DBLE,DIMAG,NINT
C     ..
C     .. Common blocks ..
      COMMON BETA,IPROB
C     ..

      WRITE (6,FMT=*) 'Specify problem:'
      READ (5,FMT=*) IPROB
      IF (IPROB.LE.0.OR.IPROB.GT.3) STOP
      IF (IPROB.EQ.1) THEN
       BETA = 10.D0
       B = 2.D0*ATAN(1.D0)
       A = -B
      ELSE IF (IPROB.EQ.2) THEN
       A = 0.D0
       B = 10.D0
      ELSE 
       A = 0.D0
       B = 1.D0
      END IF
      WRITE (6,FMT=*) 'Specify meshing tolerance:'
      READ (5,FMT=*) TOL
      NMESH = IMESH
      ELAM = CMPLX(1.D0,1.D0)
      CALL GETMSH(A,B,XMESH,NMESH,IMATCH,ELAM,SL4COF,SL4BCS,
     &            TOL,IFAIL)

      DO 20 I = 0,NMESH
          WRITE (6,FMT=*) 'I,XMESH(I):',I,XMESH(I)
   20 CONTINUE

      NLOC = 10
      DO 30 I = 0,NLOC
          XLOC(I) = (DBLE(NLOC-I)*A+DBLE(I)*B)/DBLE(NLOC)
   30 CONTINUE
      ILOC = NINT(DBLE(NLOC)/2.D0)
      IMATCH = NINT(DBLE(NMESH)/2.D0)

      WRITE (6,FMT=*) 'NMESH,IMATCH in driver:',NMESH,IMATCH
      WRITE (6,FMT=*) 'Specify RELAM:'
      READ (5,FMT=*) RELAM

      ELAM1 = CMPLX(RELAM,0.D0)
      ELAM2 = CMPLX(RELAM,1.D8)
      IFAIL = 0
      IREVCM = 1

   40 CALL BXSIDE(ELAM1,ELAM2,ELAM,F,CONST,IREVCM,ZARRY1,NZMAX,NZACT1,
     &            FARRY1,NFMAX,ANSWR1,IFAIL)

      IF (IREVCM.NE.0) THEN
C Write out mesh to check

          CALL SLMISS(XMESH,NMESH,IMATCH,ELAM,SL4COF,SL4BCS,DMISS1,
     &                CONST1,IFAIL)
          IF (IFAIL.NE.0) THEN
              WRITE (6,FMT=*) 'Failure in SLMISS, IFAIL=',IFAIL
              STOP
          END IF
          CALL SLMISS(XLOC,NLOC,ILOC,ELAM,SLSMPC,SLSMBC,DMISS2,CONST2,
     &                IFAIL)
          IF (IFAIL.NE.0) THEN
              WRITE (6,FMT=*) 'Failure in SLMISS, IFAIL=',IFAIL
              STOP
          END IF
          F = DMISS1/DMISS2
          CONST = CONST1 - CONST2
C          WRITE (19,FMT=*) 'ELAM,FUN,FUNSMP:',ELAM,DMISS1,DMISS2
C          WRITE (7,FMT=*) 'ELAM,F,CONST:',ELAM,F,CONST
          GO TO 40
      END IF

      IF (IREVCM.EQ.0 .AND. IFAIL.NE.0) THEN
          WRITE (6,FMT=*) 'Failure in BXSIDE, IFAIL=',IFAIL
          STOP
      END IF

C If we are here then we have done the first integral; now
C do the second.

      ELAM1 = CMPLX(RELAM,0.D0)
      ELAM2 = CMPLX(RELAM,-1.D8)
      IFAIL = 0
      IREVCM = 1

   50 CALL BXSIDE(ELAM1,ELAM2,ELAM,F,CONST,IREVCM,ZARRY2,NZMAX,NZACT2,
     &            FARRY2,NFMAX,ANSWR2,IFAIL)

      IF (IREVCM.NE.0) THEN
          CALL SLMISS(XMESH,NMESH,IMATCH,ELAM,SL4COF,SL4BCS,DMISS1,
     &                CONST1,IFAIL)
          IF (IFAIL.NE.0) THEN
              WRITE (6,FMT=*) 'Failure in SLMISS, IFAIL=',IFAIL
              STOP
          END IF
          CALL SLMISS(XLOC,NLOC,ILOC,ELAM,SLSMPC,SLSMBC,DMISS2,CONST2,
     &                IFAIL)
          IF (IFAIL.NE.0) THEN
              WRITE (6,FMT=*) 'Failure in SLMISS, IFAIL=',IFAIL
              STOP
          END IF
          F = DMISS1/DMISS2
          CONST = CONST1 - CONST2
C          WRITE (7,FMT=*) 'ELAM,F,CONST:',ELAM,F,CONST
          GO TO 50
      END IF

      IF (IREVCM.EQ.0 .AND. IFAIL.NE.0) THEN
          WRITE (6,FMT=*) 'Failure in BXSIDE, IFAIL=',IFAIL
          STOP
      END IF

C Write out the final answer
C
      ANS = (ANSWR1-ANSWR2)/ (8.D0*ATAN(1.D0))
      WRITE (6,FMT=*) 'The final answer was',ANS
      WRITE (6,FMT=*) 'Numbers of function evaluations:',NZACT1,NZACT2

C Now decide how many eigenvalues are in a given box.
      WRITE (6,FMT=*) 'Specify bottom right corner of box:'
      READ (5,FMT=*) EBR
      WRITE (6,FMT=*) 'Specify top right corner of box:'
      READ (5,FMT=*) ETR
      WRITE (6,FMT=*) 'Specify top left corner of box:'
      READ (5,FMT=*) ETL
      WRITE (6,FMT=*) 'Specify bottom left corner of box:'
      READ (5,FMT=*) EBL

      IFAIL = 0

      CALL MKBOX(EBR,ETR,ETL,EBL,ZL,NZMAX,NZL,ZR,NZMAX,NZR,ZT,NZMAX,NZT,
     &           ZB,NZMAX,NZB,FL,NFMAX,FR,NFMAX,FT,NFMAX,FB,NFMAX,ANSL,
     &           ANSR,ANST,ANSB,ANS,XMESH,NMESH,IMATCH,SL4COF,SL4BCS,
     &           IFAIL)
      NEIGS = NINT(DIMAG(ANS))
      WRITE (6,FMT=*) 'Number of eigenvalues was',NEIGS
      WRITE (6,FMT=*) 'ANS was',ANS
      WRITE (6,FMT=*) 'Numbers NZR,NZT,NZL,NZB:',NZR(1),NZT(1),NZL(1),
     &                 NZB(1)

      IF (NEIGS.LE.NEV) THEN
      WRITE (6,FMT=*) 'Computing all eigenvalues:'
      WRITE (6,FMT=*) 'Specify tolerance:'
      TOL = 1.D-7
      READ (5,FMT=*) TOL
      IFAIL = 0
      WRITE (6,FMT=*) 'NMESH, IMATCH before calling ALLEVS:'
      WRITE (6,FMT=*) NMESH,IMATCH
      DO 77 I = 1,NEIGS
         EIGVAL(I) = CMPLX(0.D0,0.D0)
77    CONTINUE
      IEIGS(1) = NEIGS
      CALL ALLEVS(ZL,NZMAX,NZL,ZR,NZMAX,NZR,ZT,NZMAX,NZT,ZB,NZMAX,NZB,
     &            FL,NFMAX,FR,NFMAX,FT,NFMAX,FB,NFMAX,ANSL,ANSR,ANST,
     &            ANSB,ZLC1,NZMAX,ZRC1,NZMAX,ZTC1,NZMAX,ZBC1,NZMAX,
     &            FLC1,NFMAX,FRC1,NFMAX,FTC1,NFMAX,FBC1,NFMAX,
     &            ZLC2,NZMAX,ZRC2,NZMAX,ZTC2,NZMAX,ZBC2,NZMAX,
     &            FLC2,NFMAX,FRC2,NFMAX,FTC2,NFMAX,FBC2,NFMAX,
     &            XMESH,NMESH,IMATCH,SL4COF,SL4BCS,TOL,
     &            EIGVAL,NEIGS,IEIGS,IFAIL)

      WRITE (6,FMT=*) 'Here are the eigenvalues:'
      DO 17 I = 1,NEIGS
      WRITE (6,FMT=*) EIGVAL(I)
      WRITE (16,FMT=*) EIGVAL(I)
17    CONTINUE      
      WRITE (6,FMT=*) 'IFAIL was',IFAIL
      END IF

      STOP
      END

C ------------------------------------------------------------------------------



      SUBROUTINE SL4COF(X,ELAM,P,Q3,Q2,Q1,Q0)
C     .. Scalar Arguments ..
      DOUBLE COMPLEX ELAM,Q0,Q1,Q2,Q3
      DOUBLE PRECISION P,X
C     ..
C     .. Local Scalars ..
      DOUBLE PRECISION Q,S
C     ..
C     .. Intrinsic Functions ..
      INTRINSIC CMPLX,COS,SIN
C     ..

C     .. Scalars in Common ..
      DOUBLE PRECISION BETA
      INTEGER IPROB
C     ..
C     .. Common blocks ..
      COMMON BETA,IPROB
C     ..
      GOTO (10,20,30) IPROB

10    P = 1.D0
      Q3 = CMPLX(0.D0,0.D0)
      S = 2.D0* ((BETA*SIN(2.D0*X))**2) - 4.D0*BETA*COS(2.D0*X)
      Q2 = CMPLX(-S,0.D0)
      Q1 = CMPLX(0.D0,0.D0)
      Q = ((BETA**2)* (SIN(2.D0*X)**2)-2.D0*BETA*COS(2.D0*X))**2 -
     &    8.D0* (BETA**2)* ((COS(2.D0*X)**2)- (SIN(2.D0*X)**2)) -
     &    8.D0*BETA*COS(2.D0*X)
      Q0 = CMPLX(Q,1.D0) - ELAM
      RETURN

20    P = 1.D0
      Q3 = CMPLX(0.D0,0.D0)
      Q2 = CMPLX(-2.D0*X*X,-1.D0)
      Q1 = CMPLX(0.D0,0.D0)
      Q0 = CMPLX(X**4-2.D0,X**2) - ELAM
      RETURN

30    P = 1.D0
      Q3 = CMPLX(0.D0,0.D0)
      Q2 = CMPLX(0.D0,0.D0)
      Q1 = CMPLX(0.D0,0.D0)
      Q0 = CMPLX(0.D0,0.D0) - ELAM
      RETURN

      END


C -----------------------------------------------------------


      SUBROUTINE SL4BCS(IEND,X,ELAM,U,IU,V,IV)

C     .. Scalar Arguments ..
      DOUBLE COMPLEX ELAM
      DOUBLE PRECISION X
      INTEGER IEND,IU,IV
C     ..
C     .. Array Arguments ..
      DOUBLE COMPLEX U(IU,2),V(IV,2)
C     ..
C     .. Scalars in Common ..
      DOUBLE PRECISION BETA
      INTEGER IPROB
C     ..
C     .. Common blocks ..
      COMMON BETA,IPROB
C     ..
C     .. Local Scalars ..
      DOUBLE PRECISION THETA
      DOUBLE COMPLEX IOMEG1,IOMEG2,E1,E2

      GOTO (10,20,30) IPROB

10    U(1,1) = 0.D0
      U(1,2) = 0.D0
      U(2,1) = 1.D0
      U(2,2) = 0.D0
      V(1,1) = 1.D0
      V(1,2) = 1.D0
      V(2,1) = 0.D0
      V(2,2) = 0.D0
      RETURN

20    U(1,1) = 0.D0
      U(1,2) = 0.D0
      U(2,1) = 1.D0
      U(2,2) = 0.D0
      V(1,1) = 1.D0
      V(1,2) = 1.D0
      V(2,1) = 0.D0
      V(2,2) = 0.D0
      RETURN

30    THETA = ATAN(1.D0)/2.D0
      IOMEG1 = CMPLX(COS(THETA),SIN(THETA))
      IOMEG2 = 2.D0*IOMEG1
      E1 = EXP(IOMEG1*X)
      E2 = EXP(IOMEG2*X)
      WRITE (11,FMT=*) 'OMEGA1**4:',IOMEG1**4
      WRITE (11,FMT=*) 'OMEGA2**4:',IOMEG2**4


      U(1,1) = CMPLX(1.D0,0.D0)*E1
      U(1,2) = CMPLX(1.D0,0.D0)*E2
      U(2,1) = IOMEG1*E1
      U(2,2) = IOMEG2*E2
      V(1,1) = -E1*(IOMEG1**3)
      V(1,2) = -E2*(IOMEG2**3)
      V(2,1) = E1*IOMEG1**2
      V(2,2) = E2*IOMEG2**2

      WRITE (11,FMT=*) 'X, LAMBDA:',X,ELAM
      WRITE (11,FMT=*) 'U:'
      WRITE (11,FMT=*) U(1,1),U(1,2)
      WRITE (11,FMT=*) U(2,1),U(2,2)
      WRITE (11,FMT=*) 'V:'
      WRITE (11,FMT=*) V(1,1),V(1,2)
      WRITE (11,FMT=*) V(2,1),V(2,2)
      WRITE (12,FMT=*) 'Initial x, detU, detV:',
     &      X,U(1,1)*U(2,2)-U(1,2)*U(2,1),
     &      V(1,1)*V(2,2)-V(1,2)*V(2,1)


C      U(1,1) = CMPLX(0.D0,0.D0)
C      U(1,2) = CMPLX(1.D0,0.D0)
C      U(2,1) = CMPLX(-1.D0,0.D0)
C      U(2,2) = CMPLX(0.D0,0.D0)
CC Theta = Pi/8
C      THETA = 0.3926990816987238D0
C
C      V(2,1) = CMPLX(-3.D0*COS(5.D0*THETA),-3.D0*SIN(5.D0*THETA))
C      V(2,2) = CMPLX(SQRT(2.D0),SQRT(2.D0))
C      V(1,2) = CMPLX(-6.D0*COS(7.D0*THETA),-6.D0*SIN(7.D0*THETA))
C      V(1,1) = CMPLX(-3.5D0,0.D0)*V(2,2)

      RETURN

      END


C ------------------------------------------------------------------

      SUBROUTINE SLSMPC(X,ELAM,P,Q3,Q2,Q1,Q0)

C     .. Scalar Arguments ..
      DOUBLE COMPLEX ELAM,Q0,Q1,Q2,Q3
      DOUBLE PRECISION P,X
C     ..
C     .. Intrinsic Functions ..
      INTRINSIC CMPLX,DIMAG
C     ..
      P = 1.D0
      Q3 = CMPLX(0.D0,0.D0)
      Q2 = CMPLX(0.D0,0.D0)
      Q1 = CMPLX(0.D0,0.D0)
      Q0 = CMPLX(1.D0,-DIMAG(ELAM))

      RETURN
      END

C ------------------------------------------------------------------


      SUBROUTINE SLSMBC(ISING,X,ELAM,U,IU,V,IV)

C     .. Scalar Arguments ..
      DOUBLE COMPLEX ELAM
      DOUBLE PRECISION X
      INTEGER ISING,IU,IV
C     ..
C     .. Array Arguments ..
      DOUBLE COMPLEX U(IU,2),V(IV,2)
C     ..
C     .. Intrinsic Functions ..
      INTRINSIC CMPLX
C     ..

      U(1,1) = CMPLX(1.D0,0.D0)
      U(2,1) = CMPLX(0.D0,0.D0)
      U(1,2) = CMPLX(0.D0,0.D0)
      U(2,2) = CMPLX(1.D0,0.D0)
      V(1,1) = CMPLX(0.D0,0.D0)
      V(2,1) = CMPLX(0.D0,0.D0)
      V(1,2) = CMPLX(0.D0,0.D0)
      V(2,2) = CMPLX(0.D0,0.D0)

      RETURN
      END

C ------------------------------------------------------------------
