      PROGRAM MAIN

C     .. Parameters ..
      INTEGER NZMAX,NFMAX,IMESH,IMLOC
      PARAMETER (NZMAX=10000,NFMAX=NZMAX,IMESH=10000,IMLOC=10)
C     ..
C     .. Local Scalars ..
      DOUBLE COMPLEX ANS,ANSB,ANSBC,ANSL,ANSLC,ANSR,ANSRC,ANST,ANSTC,
     &               ANSWR1,ANSWR2,DMISS1,DMISS2,EBL,EBR,EIGVAL,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,NZB,NZBC,NZL,NZLC,NZR,NZRC,NZT,NZTC
C     ..
C     .. Local Arrays ..
      DOUBLE COMPLEX FARRY1(0:NFMAX,1:2),FARRY2(0:NFMAX,1:2),
     &               FB(0:NFMAX,1:2),FBC(0:NFMAX,1:2),FL(0:NFMAX,1:2),
     &               FLC(0:NFMAX,1:2),FR(0:NFMAX,1:2),FRC(0:NFMAX,1:2),
     &               FT(0:NFMAX,1:2),FTC(0:NFMAX,1:2),ZARRY1(0:NZMAX),
     &               ZARRY2(0:NZMAX),ZB(0:NZMAX),ZBC(0:NZMAX),
     &               ZL(0:NZMAX),ZLC(0:NZMAX),ZR(0:NZMAX),ZRC(0:NZMAX),
     &               ZT(0:NZMAX),ZTC(0:NZMAX)
      DOUBLE PRECISION XLOC(0:IMLOC),XMESH(0:IMESH)
C     ..
C     .. Scalars in Common ..
      DOUBLE PRECISION BETA
C     ..
C     .. External Subroutines ..
      EXTERNAL BCHILD,BXSIDE,LCHILD,MKBOX,RCHILD,RFNBOX,SL6BCS,SL6COF,
     &         SLMISS,SLSMBC,SLSMPC,TCHILD
C     ..
C     .. Intrinsic Functions ..
      INTRINSIC ABS,ATAN,CMPLX,DBLE,DIMAG,NINT
C     ..
C     .. Common blocks ..
      COMMON BETA
C     ..
      BETA = 10.D0
      WRITE (6,FMT=*) 'Specify BETA:'
      READ (5,FMT=*) BETA
      B = 2.D0*ATAN(1.D0)
      A = -B
      WRITE (6,FMT=*) 'Specify endpoints A and B:'
      READ (5,FMT=*) A,B
      WRITE (6,FMT=*) 'Specify no. of meshpoints:'
      READ (5,FMT=*) NMESH
      XMESH(0) = A/ (1.D0+ABS(A))
      XMESH(NMESH) = B/ (1.D0+ABS(B))
      DO 10 I = 1,NMESH - 1
          XMESH(I) = (DBLE(NMESH-I)*XMESH(0)+DBLE(I)*XMESH(NMESH))/
     &               DBLE(NMESH)
   10 CONTINUE
      DO 20 I = 0,NMESH
          XMESH(I) = XMESH(I)/ (1.D0-ABS(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=*) '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,SL6COF,SL6BCS,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,SL6COF,SL6BCS,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,SL6COF,SL6BCS,
     &           IFAIL)

      WRITE (6,FMT=*) 'Number of eigenvalues was',NINT(DIMAG(ANS))
      WRITE (6,FMT=*) 'ANS was',ANS
      WRITE (6,FMT=*) 'Numbers NZR,NZT,NZL,NZB:',NZR,NZT,NZL,NZB

C Now try to make a child box and compute the number of eigenvalues it
C contains.
      IFAIL = 0
      CALL TCHILD(ZL,NZMAX,NZL,ZR,NZMAX,NZR,ZT,NZMAX,NZT,ZB,NZMAX,NZB,
     &            FL,NFMAX,FR,NFMAX,FT,NFMAX,FB,NFMAX,ANSL,ANSR,ANST,
     &            ANSB,ZLC,NZMAX,NZLC,ZRC,NZMAX,NZRC,ZTC,NZMAX,NZTC,ZBC,
     &            NZMAX,NZBC,FLC,NFMAX,FRC,NFMAX,FTC,NFMAX,FBC,NFMAX,
     &            ANSLC,ANSRC,ANSTC,ANSBC,XMESH,NMESH,IMATCH,SL6COF,
     &            SL6BCS,IFAIL)
      ANS = (ANSLC+ANSRC+ANSTC+ANSBC)/ (8.D0*ATAN(1.D0))
      WRITE (6,FMT=*) 'Number of eigenvalues in top:',NINT(DIMAG(ANS))
      WRITE (6,FMT=*) 'ANS is now',ANS
      WRITE (6,FMT=*) 'Numbers NZRC,NZTC,NZLC,NZBC:',NZRC,NZTC,NZLC,NZBC

      CALL BCHILD(ZL,NZMAX,NZL,ZR,NZMAX,NZR,ZT,NZMAX,NZT,ZB,NZMAX,NZB,
     &            FL,NFMAX,FR,NFMAX,FT,NFMAX,FB,NFMAX,ANSL,ANSR,ANST,
     &            ANSB,ZLC,NZMAX,NZLC,ZRC,NZMAX,NZRC,ZTC,NZMAX,NZTC,ZBC,
     &            NZMAX,NZBC,FLC,NFMAX,FRC,NFMAX,FTC,NFMAX,FBC,NFMAX,
     &            ANSLC,ANSRC,ANSTC,ANSBC,XMESH,NMESH,IMATCH,SL6COF,
     &            SL6BCS,IFAIL)
      ANS = (ANSLC+ANSRC+ANSTC+ANSBC)/ (8.D0*ATAN(1.D0))
      WRITE (6,FMT=*) 'Number of eigenvalues in bottom:',
     &  NINT(DIMAG(ANS))
      WRITE (6,FMT=*) 'ANS is now',ANS
      WRITE (6,FMT=*) 'Numbers NZRC,NZTC,NZLC,NZBC:',NZRC,NZTC,NZLC,NZBC

      CALL LCHILD(ZL,NZMAX,NZL,ZR,NZMAX,NZR,ZT,NZMAX,NZT,ZB,NZMAX,NZB,
     &            FL,NFMAX,FR,NFMAX,FT,NFMAX,FB,NFMAX,ANSL,ANSR,ANST,
     &            ANSB,ZLC,NZMAX,NZLC,ZRC,NZMAX,NZRC,ZTC,NZMAX,NZTC,ZBC,
     &            NZMAX,NZBC,FLC,NFMAX,FRC,NFMAX,FTC,NFMAX,FBC,NFMAX,
     &            ANSLC,ANSRC,ANSTC,ANSBC,XMESH,NMESH,IMATCH,SL6COF,
     &            SL6BCS,IFAIL)
      ANS = (ANSLC+ANSRC+ANSTC+ANSBC)/ (8.D0*ATAN(1.D0))
      WRITE (6,FMT=*) 'Number of eigenvalues in left:',NINT(DIMAG(ANS))
      WRITE (6,FMT=*) 'ANS is now',ANS
      WRITE (6,FMT=*) 'Numbers NZRC,NZTC,NZLC,NZBC:',NZRC,NZTC,NZLC,NZBC

      CALL RCHILD(ZL,NZMAX,NZL,ZR,NZMAX,NZR,ZT,NZMAX,NZT,ZB,NZMAX,NZB,
     &            FL,NFMAX,FR,NFMAX,FT,NFMAX,FB,NFMAX,ANSL,ANSR,ANST,
     &            ANSB,ZLC,NZMAX,NZLC,ZRC,NZMAX,NZRC,ZTC,NZMAX,NZTC,ZBC,
     &            NZMAX,NZBC,FLC,NFMAX,FRC,NFMAX,FTC,NFMAX,FBC,NFMAX,
     &            ANSLC,ANSRC,ANSTC,ANSBC,XMESH,NMESH,IMATCH,SL6COF,
     &            SL6BCS,IFAIL)
      ANS = (ANSLC+ANSRC+ANSTC+ANSBC)/ (8.D0*ATAN(1.D0))
      WRITE (6,FMT=*) 'Number of eigenvalues in right:',NINT(DIMAG(ANS))
      WRITE (6,FMT=*) 'ANS is now',ANS
      WRITE (6,FMT=*) 'Numbers NZRC,NZTC,NZLC,NZBC:',NZRC,NZTC,NZLC,NZBC

      IFAIL = 0
      TOL = 1.0D-6
      WRITE (6,FMT=*) 'Specify eigenvalue tolerance'
      READ (5,FMT=*) TOL

      CALL RFNBOX(ZL,NZMAX,NZL,ZR,NZMAX,NZR,ZT,NZMAX,NZT,ZB,NZMAX,NZB,
     &            FL,NFMAX,FR,NFMAX,FT,NFMAX,FB,NFMAX,ANSL,ANSR,ANST,
     &            ANSB,ZLC,NZMAX,NZLC,ZRC,NZMAX,NZRC,ZTC,NZMAX,NZTC,ZBC,
     &            NZMAX,NZBC,FLC,NFMAX,FRC,NFMAX,FTC,NFMAX,FBC,NFMAX,
     &            ANSLC,ANSRC,ANSTC,ANSBC,XMESH,NMESH,IMATCH,SL6COF,
     &            SL6BCS,TOL,'C',EIGVAL,IFAIL)
      IF (IFAIL.NE.0) THEN
          WRITE (6,FMT=*) 'Failure exit from RFNBOX, IFAIL=',IFAIL
          STOP
      END IF
      WRITE (6,FMT=*) 'Final eigenvalue was',EIGVAL

      STOP
      END

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



      SUBROUTINE SL6COF(X,ELAM,P,Q5,Q4,Q3,Q2,Q1,Q0)
C     .. Scalar Arguments ..
      DOUBLE COMPLEX ELAM,Q0,Q1,Q2,Q3,Q4,Q5
      DOUBLE PRECISION P,X
C     ..
C     .. Local Scalars ..
      DOUBLE PRECISION BETA2,BETA4,BETA6,COS2X,Q,SIN2X,
     &                 QLOC,D2Q
C     ..
C     .. Intrinsic Functions ..
      INTRINSIC CMPLX,COS,SIN
C     ..

C     .. Scalars in Common ..
      DOUBLE PRECISION BETA
C     ..
C     .. Common blocks ..
      COMMON BETA
C     ..

      BETA2 = BETA*BETA
      BETA4 = BETA2*BETA2
      BETA6 = BETA2*BETA4
      SIN2X = SIN(2.D0*X)
      COS2X = COS(2.D0*X)
      P = 1.D0
      QLOC = (-2.D0*BETA*COS2X+ (BETA*SIN2X)**2)
      D2Q = 8.D0*BETA*(COS2X*(1.D0+BETA*COS2X)-BETA*SIN2X*SIN2X)
      Q5 = CMPLX(0.D0,0.D0)
      Q3 = CMPLX(0.D0,0.D0)
      Q1 = CMPLX(0.D0,0.D0)
      Q4 = CMPLX(3.D0*QLOC,0.D0)
      Q2 = CMPLX(4.D0*D2Q - 3.D0*QLOC*QLOC,0.D0)
      Q = QLOC*(QLOC*QLOC-D2Q) 
     &  - 32.D0*BETA2*(SIN2X**2)*(1.D0+BETA*COS2X)**2
     &  - 2.D0*D2Q*BETA*(-2.D0*COS2X+BETA*SIN2X*SIN2X)
     &  - 32.D0*BETA*COS2X 
     &  - 128.D0*BETA*BETA*COS(4.D0*X)

      P = 1.D0
C      Q5 = CMPLX(0.D0,0.D0)
C      Q4 = CMPLX(0.D0,0.D0)
C      Q3 = CMPLX(0.D0,0.D0)
C      Q2 = CMPLX(0.D0,0.D0)
C      Q1 = CMPLX(0.D0,0.D0)
C      Q = 0.D0
      Q0 = CMPLX(Q,1.D0) - ELAM

      RETURN
      END


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


      SUBROUTINE SL6BCS(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      U(1,1) = CMPLX(0.D0,0.D0)
C      U(1,2) = CMPLX(0.D0,0.D0)
C      U(2,1) = CMPLX(0.D0,0.D0)
C      U(2,2) = CMPLX(0.D0,0.D0)

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

C      U(1,1) = (0.D0,0.D0)
C      U(1,2) = (0.D0,0.D0)
C      U(1,3) = (0.D0,0.D0)
C      U(2,1) = (0.D0,0.D0)
C      U(2,2) = (1.D0,0.D0)
C      U(2,3) = (0.D0,0.D0)
C      U(3,1) = (0.D0,0.D0)
C      U(3,2) = (0.D0,0.D0)
C      U(3,3) = (0.D0,0.D0)

C      V(1,1) = (1.D0,0.D0)
C      V(1,2) = (0.D0,0.D0)
C      V(1,3) = (0.D0,0.D0)
C      V(2,1) = (0.D0,0.D0)
C      V(2,2) = (0.D0,0.D0)
C      V(2,3) = (0.D0,0.D0)
C      V(3,1) = (0.D0,0.D0)
C      V(3,2) = (0.D0,0.D0)
C      V(3,3) = (1.D0,0.D0)

      U(1,1) = (0.D0,0.D0)
      U(1,2) = (0.D0,0.D0)
      U(1,3) = (0.D0,0.D0)
      U(2,1) = (0.D0,0.D0)
      U(2,2) = (1.D0,0.D0)
      U(2,3) = (0.D0,0.D0)
      U(3,1) = (0.D0,0.D0)
      U(3,2) = (0.D0,0.D0)
      U(3,3) = (0.D0,0.D0)

      V(1,1) = (1.D0,0.D0)
      V(1,2) = (0.D0,0.D0)
      V(1,3) = (0.D0,0.D0)
      V(2,1) = (0.D0,0.D0)
      V(2,2) = (0.D0,0.D0)
      V(2,3) = (0.D0,0.D0)
      V(3,1) = (0.D0,0.D0)
      V(3,2) = (0.D0,0.D0)
      V(3,3) = (1.D0,0.D0)

      RETURN
      END


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

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

C     .. Scalar Arguments ..
      DOUBLE COMPLEX ELAM,Q0,Q1,Q2,Q3,Q4,Q5
      DOUBLE PRECISION P,X
C     ..
C     .. Intrinsic Functions ..
      INTRINSIC CMPLX,DIMAG
C     ..
      P = 1.D0
      Q5 = CMPLX(0.D0,0.D0)
      Q4 = CMPLX(0.D0,0.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     ..


C Set even-order derivatives to zero at endpoints:
      U(1,1) = (0.D0,0.D0)
      U(1,2) = (0.D0,0.D0)
      U(1,3) = (0.D0,0.D0)
      U(2,1) = (0.D0,0.D0)
      U(2,2) = (1.D0,0.D0)
      U(2,3) = (0.D0,0.D0)
      U(3,1) = (0.D0,0.D0)
      U(3,2) = (0.D0,0.D0)
      U(3,3) = (0.D0,0.D0)

      V(1,1) = (1.D0,0.D0)
      V(1,2) = (0.D0,0.D0)
      V(1,3) = (0.D0,0.D0)
      V(2,1) = (0.D0,0.D0)
      V(2,2) = (0.D0,0.D0)
      V(2,3) = (0.D0,0.D0)
      V(3,1) = (0.D0,0.D0)
      V(3,2) = (0.D0,0.D0)
      V(3,3) = (1.D0,0.D0)


      RETURN
      END
