SUBROUTINE STRSCH(XL,XR,NUMEIG,TOL,XMESH,NMESH,IMATCH,SL4COF, & SL4BCS,ZL,IZL,NZL,ZR,IZR,NZR,ZT,IZT,NZT,ZB, & IZB,NZB,FL,IFL,FR,IFR,FT,IFT,FB,IFB,ANSL, & ANSR,ANST,ANSB,ZLC1,ZRC1,ZTC1,ZBC1,ZLC2, & ZRC2,ZBC2,ZTC2,FLC1,FRC1,FTC1,FBC1,FLC2, & FRC2,FBC2,FTC2,NEV,EIGS,IEIGS,IFAIL) C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C This routine STRSCH (Strip Search) computes all eigenvalues in a C vertical strip specified by the x-coordinates XL, XR. C C Input Variables: XL,XR,NUMEIG C C NUMEIG is the number of eigenvalues in the strip. This is C calculated using NZERO, and is input into STRSCH. C C Output Variables: EIGS, IFAIL. C C EIGS is an array containing the eigenvalues in the strip. C C EBR,ETR,ETL, and EBL are the corners of boxes used in STRSCH. C These are passed to MKBOX. NUMLOC is the number of eigenvalues C found in the box by MKBOX. This number is passed to ALLEVS. C NUMFND is the cumulative number of eigenvalues found in the C strip. STRSCH exits when either NUMFND = NUMEIG, or the next C box would be higher than HMAX. C C All other variables in the parameter list are passed to and from C MKBOX and ALLEVS, and are used in those subroutines. C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C C .. Local Scalars .. INTEGER ICOUNT,IDIRECT,IFAIL,IMATCH,IFB,IFL,IFR,IFT, & IZB,IZL,IZR,IZT,NMESH,NUMEIG,NUMLOC,NUMFND DOUBLE PRECISION XL,XR,YB,YBOLD,YT,HEIGHT,WIDTH,HMAX,TOL DOUBLE COMPLEX EBL,EBR,ETL,ETR,ANS C .. Local Arrays .. DOUBLE COMPLEX FB(0:IFB,1:2,NEV),FBC1(0:IFB,1:2), & FL(0:IFL,1:2,NEV),FLC1(0:IFL,1:2), & FR(0:IFR,1:2,NEV),FRC1(0:IFR,1:2), & FT(0:IFT,1:2,NEV),FTC1(0:IFT,1:2), & ZB(0:IZB,NEV),ZBC1(0:IZB), & ZL(0:IZL,NEV),ZLC1(0:IZL),ZR(0:IZR,NEV), & ZRC1(0:IZR),ZT(0:IZT,NEV),ZTC1(0:IZT), & FBC2(0:IFB,1:2),FRC2(0:IFR,1:2),FTC2(0:IFT,1:2), & FLC2(0:IFL,1:2),ZBC2(0:IZB),ZLC2(0:IZL), & ZRC2(0:IZR),ZTC2(0:IZT),EIGS(NEV) DOUBLE COMPLEX ANSL(NEV),ANSR(NEV),ANST(NEV),ANSB(NEV) INTEGER NZL(NEV),NZR(NEV),NZT(NEV),NZB(NEV),IEIGS(NEV) DOUBLE PRECISION XMESH(0:NMESH) C C C23456789112345678921234567893123456789412345678951234567896123456789712 C C EXTERNAL ALLEVS,MKBOX,SL4BCS,SL4COF INTRINSIC CMPLX,AIMAG,NINT IF (NUMEIG.EQ.0) THEN WRITE(6,FMT=*) 'NUMEIG =' ,0 STOP ENDIF HMAX = 1.D8 WIDTH = XR - XL HEIGHT = DMAX1(2.D0*WIDTH,1.D0) YB = 0.D0 YT = HEIGHT J = 1 IDIRECT = 1 ICOUNT = 0 ANS = CMPLX(0.D0,0.D0) NUMFND = 0 30 EBR = CMPLX(XR,YB) ETR = CMPLX(XR,YT) ETL = CMPLX(XL,YT) EBL = CMPLX(XL,YB) CALL MKBOX(EBR,ETR,ETL,EBL,ZL,IZL,NZL(1),ZR,IZR,NZR(1),ZT,IZT, & NZT(1),ZB,IZB,NZB(1),FL,IFL,FR,IFR,FT,IFT,FB,IFB, & ANSL(1),ANSR(1),ANST(1),ANSB(1),ANS,XMESH,NMESH, & IMATCH,SL4COF,SL4BCS, & IFAIL) IF (IFAIL.NE.0) THEN WRITE(6,FMT=*) 'For MKBOX' WRITE(6,FMT=*) 'IFAIL =' ,IFAIL STOP ENDIF ICOUNT = ICOUNT + 1 IF (IDIRECT.EQ.1) THEN YBOLD = YB YB = -YT YT = -YBOLD IDIRECT = 2 ELSE YBOLD = YB YB = -YB YT = YB + 2.D0*(YT - YBOLD) IDIRECT = 1 ENDIF NUMLOC = NINT(AIMAG(ANS)) WRITE(6,FMT=*) 'NUMLOC = ' ,NUMLOC IF (NINT(AIMAG(ANS)).EQ.0) GOTO 30 IFAIL = 0 IEIGS(1) = NUMLOC DO 35 I = 1, NUMLOC EIGS(NUMFND+I) = CMPLX(0.D0,0.D0) 35 CONTINUE CALL ALLEVS(ZL,IZL,NZL,ZR,IZR,NZR,ZT,IZT,NZT,ZB,IZB,NZB, & FL,IFL,FR,IFR,FT,IFT,FB,IFB,ANSL,ANSR,ANST, & ANSB,ZLC1,IZL,ZRC1,IZR,ZTC1,IZT,ZBC1,IZB, & FLC1,IFL,FRC1,IFR,FTC1,IFT,FBC1,IFB, & ZLC2,IZL,ZRC2,IZR,ZTC2,IZT,ZBC2,IZB, & FLC2,IFL,FRC2,IFR,FTC2,IFT,FBC2,IFB, & XMESH,NMESH,IMATCH,SL4COF,SL4BCS,TOL, & EIGS(NUMFND+1),NEV,IEIGS,IFAIL) IF (IFAIL.NE.0) THEN WRITE(6,FMT=*) 'For ALLEVS' WRITE(6,FMT=*) 'IFAIL =' ,IFAIL STOP ENDIF WRITE(6,FMT=*) 'Eigenvalues found in current box:' DO 40 I = 1, NUMLOC WRITE(6,FMT=*) EIGS(NUMFND+I) C WRITE(100,FMT=*) EIGS(NUMFND+I) 40 CONTINUE NUMFND = NUMFND + NUMLOC IF (NUMFND.EQ.NUMEIG) THEN RETURN ELSEIF (YT.GE.HMAX) THEN WRITE(6,FMT=*) 'Height exceeds' , HMAX WRITE(6,FMT=*) 'Number of eigenvalues in the strip is' ,NUMEIG WRITE(6,FMT=*) 'Number of eigenvalues found is' ,NUMFND WRITE(6,FMT=*) 'Here are the eigenvalues found' DO 60 I = 1, NUMFND WRITE(6,FMT=*) EIGS(I) 60 CONTINUE STOP ELSE GOTO 30 ENDIF END C ------------------------------------------------------------------------- RECURSIVE SUBROUTINE ALLEVS(ZL,IZL,NZL,ZR,IZR,NZR,ZT,IZT,NZT, & ZB,IZB,NZB,FL, & IFL,FR,IFR,FT,IFT,FB,IFB,ANSL,ANSR,ANST,ANSB, & ZLC1,IZLC1,ZRC1,IZRC1,ZTC1,IZTC1,ZBC1,IZBC1, & FLC1,IFLC1,FRC1,IFRC1,FTC1,IFTC1,FBC1,IFBC1, & ZLC2,IZLC2,ZRC2,IZRC2,ZTC2,IZTC2,ZBC2,IZBC2, & FLC2,IFLC2,FRC2,IFRC2,FTC2,IFTC2,FBC2,IFBC2, & XMESH,NMESH,IMATCH,SL4COF,SL4BCS,TOL, & EIGVAL,NEIGS,IEIGS,IFAIL) C This routine computes all the eigenvalues in a box specified by the C arrays ZL, ZR, ZT, ZB, FL, FR, FT, FB. C C The arguments are as follows. C C ZL(0:IZL,1:NEIGS): double complex array such that the elements ZL(0:NZL(1),1) C specify the points on the left hand side of the box, starting C from the top. These elements must be assigned before calling C ALLEVS, for example by a call to MKBOX. The elements C ZL(0:IZL,2:NEIGS) are for workspace and need not be assigned on C entry. C C On exit, the whole of the array ZL is liable to have been C overwritten. C C IZL: integer, the first dimension of ZL as declared in the calling (sub)program. C Unchanged on exit. C C NZL(1:NEIGS): integer array. The nodes on the left hand side of the box are in C ZL(0:NZL(1),1). C C NZL is overwritten on exit. C C ZR(0:IZR,1:NEIGS): double complex array such that the elements ZR(0:NZR(1),1) C specify the points on the right hand side of the box, starting C from the bottom. These elements must be assigned before calling C ALLEVS, for example by a call to MKBOX. The elements C ZR(0:IZR,2:NEIGS) are for workspace and need not be assigned on C entry. C C On exit, the whole of the array ZR is liable to have been C overwritten. C C IZR: integer, the first dimension of ZR as declared in the calling (sub)program. C Unchanged on exit. C C NZR(1:NEIGS): integer array. The nodes on the right hand side of the box are in C ZR(0:NZR(1),1). C C NZR is overwritten on exit. C C ZT(0:IZT,1:NEIGS): double complex array such that the elements ZT(0:NZT(1),1) C specify the points on the top side of the box, starting C from the right. These elements must be assigned before calling C ALLEVS, for example by a call to MKBOX. The elements C ZT(0:IZT,2:NEIGS) are for workspace and need not be assigned on C entry. C C On exit, the whole of the array ZT is liable to have been C overwritten. C C IZT: integer, the first dimension of ZT as declared in the calling (sub)program. C Unchanged on exit. C C NZT(1:NEIGS): integer array. The nodes on the top side of the box are in C ZT(0:NZT(1),1). C C NZT is overwritten on exit. C C C ZB(0:IZB,1:NEIGS): double complex array such that the elements ZB(0:NZB(1),1) C specify the points on the bottom side of the box, starting C from the left. These elements must be assigned before calling C ALLEVS, for example by a call to MKBOX. The elements C ZB(0:IZB,2:NEIGS) are for workspace and need not be assigned on C entry. C C On exit, the whole of the array ZB is liable to have been C overwritten. C C IZB: integer, the first dimension of ZB as declared in the calling (sub)program. C Unchanged on exit. C C NZB(1:NEIGS): integer array. The nodes on the left hand side of the box are in C ZB(0:NZB(1),1). C C NZB is overwritten on exit. C C FL(0:IFL,1:2,1:NEIGS): double complex array such that the elements FL(0:NZL(1),1:2,1) C specify the function values at points on the left hand side of the box, starting C from the top. These elements must be assigned before calling C ALLEVS, for example by a call to MKBOX. The elements C FL(0:IFL,1:2,2:NEIGS) are for workspace and need not be assigned on C entry. C C On exit, the whole of the array FL is liable to have been C overwritten. C C IFL: integer, the first dimension of FL as declared in the calling (sub)program. C Unchanged on exit. C C C FR(0:IFR,1:2,1:NEIGS): double complex array such that the elements FR(0:NZR(1),1:2,1) C specify the function values at points on the right hand side of the box, starting C from the bottom. These elements must be assigned before calling C ALLEVS, for example by a call to MKBOX. The elements C FR(0:IFR,1:2,2:NEIGS) are for workspace and need not be assigned on C entry. C C On exit, the whole of the array FR is liable to have been C overwritten. C C IFR: integer, the first dimension of FR as declared in the calling (sub)program. C Unchanged on exit. C C C FT(0:IFT,1:2,1:NEIGS): double complex array such that the elements FT(0:NZT(1),1:2,1) C specify the function values at points on the top side of the box, starting C from the right. These elements must be assigned before calling C ALLEVS, for example by a call to MKBOX. The elements C FT(0:IFT,1:2,2:NEIGS) are for workspace and need not be assigned on C entry. C C On exit, the whole of the array FT is liable to have been C overwritten. C C IFT: integer, the first dimension of FT as declared in the calling (sub)program. C Unchanged on exit. C C C FB(0:IFB,1:2,1:NEIGS): double complex array such that the elements FB(0:NZB(1),1:2,1) C specify the function values at points on the bottom side of the box, starting C from the left. These elements must be assigned before calling C ALLEVS, for example by a call to MKBOX. The elements C FT(0:IFB,1:2,2:NEIGS) are for workspace and need not be assigned on C entry. C C On exit, the whole of the array FB is liable to have been C overwritten. C C IFB: integer, the first dimension of FB as declared in the calling (sub)program. C Unchanged on exit. C C ANSL(NEIGS): double complex array, such that ANSL(1) has been assigned to C the appropriate value corresponding to ZL and FL by the preceeding C call to MKBOX. The other elements need not have been assigned before C entry. C C All elements of ANSL may be overwritten by ALLEVS. C C ANSR(NEIGS): double complex array, such that ANSL(1) has been assigned to C the appropriate value corresponding to ZR and FR61q by the preceeding C call to MKBOX. The other elements need not have been assigned before C entry. C C All elements of ANSR may be overwritten by ALLEVS. C C ANST(NEIGS): double complex array, such that ANST(1) has been assigned to C the appropriate value corresponding to ZT and FT by the preceeding C call to MKBOX. The other elements need not have been assigned before C entry. C C All elements of ANST may be overwritten by ALLEVS. C C ANSB(NEIGS): double complex array, such that ANSB(1) has been assigned to C the appropriate value corresponding to ZB and FB by the preceeding C call to MKBOX. The other elements need not have been assigned before C entry. C C All elements of ANSB may be overwritten by ALLEVS. C C ZLC1(0:IZLC1), FLC1(0:IFLC1,1:2), ZRC1(0:IZRC1), FRC1(0:IFRC1,1:2), C ZTC1(0:IZTC1), FTC1(0:IFTC1,1:2), ZBC1(0:IZBC1), FBC1(0:IFBC1,1:2), C ZLC2(0:IZLC2), FLC2(0:IFLC2,1:2), ZRC2(0:IZRC2), FRC2(0:IFRC2,1:2), C ZTC2(0:IZTC2), FTC2(0:IFTC2,1:2), ZBC2(0:IZBC2), FBC2(0:IFBC2,1:2): C double complex workspace arrays of the dimensions shown. These C need not be assigned on entry. C C IZLC1, IZRC1, IZTC1, IZBC1, IFLC1, IFRC1, IFBC1, IFTC1, C IZLC2, IZRC2, IZTC2, IZBC2, IFLC2, IFRC2, IFBC2, IFTC2: C integers specifying the array dimensions. In general these should C be large enough to ensure that C C IZLC1 > = NZL(1), IZLC2 > = NZL(1), C IZRC1 > = NZR(1), IZRC2 > = NZR(1), C IZTC1 > = NZT(1), IZTC2 > = NZT(1), C IZBC1 > = NZB(1), IZBC2 > = NZB(1), C IFLC1 > = NZL(1), IFLC2 > = NZL(1), C IFRC1 > = NZR(1), IFRC2 > = NZR(1), C IFTC1 > = NZT(1), IFTC2 > = NZT(1), C IFBC1 > = NZB(1), IFBC2 > = NZB(1). C C XMESH(0:NMESH): double precision real; the mesh used for the X-integration. C Unchanged on exit. C C NMESH: integer, number of mesh intervals for X-integration. C C IMATCH: integer, index of shooting matchpoint XMESH(IMATCH). C C SL4COF, SL4BCS: the subroutines which define the coefficients C and boundary conditions for the problem. C C TOL: double precision real. TOL > 0 must be set on input and C speficies the accuracy for the rootfinding. C C C EIGVAL(NEIGS): double complex array, *output*. On exit this array C contains the eigenvalues found. WARNING: if a box contains several C eigenvalues and the routine does not manage to find all of them, C then several elements of EIGVAL may be undefined on exit. C C NEIGS: integer. On input, NEIGS should specify the last dimensions C of various arrays (see above). NEIGS should be greater than IEIGS(1) C (see below). NEIGS is unchanged on exit. C C IEIGS(NEIGS): integer array. On entry, IEIGS(1) must specify the number C of eigenvalues in the box, while the other elements need not be C defined. On exit, the whole of the array IEIGS is liable to have C been overwritten. C C IFAIL: integer, error flag. To be set to 0 on entry. A non-zero C value on exit denotes an error. C C IMPLICIT NONE C .. Scalar Arguments .. INTEGER IZL,IZR,IZT,IZB,IFL,IFR,IFT,IFB,NMESH,IMATCH, & NEIGS,IFAIL,IZLC1,IZRC1,IZTC1,IZBC1,IZLC2,IZRC2,IZTC2,IZBC2, & IFLC1,IFLC2,IFRC1,IFRC2,IFTC1,IFTC2,IFBC1,IFBC2 DOUBLE PRECISION TOL C .. C .. Array Arguments .. INTEGER NZL(NEIGS),NZR(NEIGS),NZT(NEIGS),NZB(NEIGS),IEIGS(*) DOUBLE COMPLEX ANSL(NEIGS),ANSR(NEIGS),ANST(NEIGS),ANSB(NEIGS), & EIGVAL(NEIGS),ZL(0:IZL,NEIGS),ZR(0:IZR,NEIGS),ZT(0:IZT,NEIGS), & ZB(0:IZB,NEIGS),FL(0:IFL,2,NEIGS),FR(0:IFR,2,NEIGS), & FT(0:IFT,2,NEIGS),FB(0:IFB,2,NEIGS) DOUBLE COMPLEX ZLC1(0:IZLC1),ZRC1(0:IZRC1),ZTC1(0:IZTC1), & ZBC1(0:IZBC1),ZLC2(0:IZLC2),ZRC2(0:IZRC2),ZTC2(0:IZTC2), & ZBC2(0:IZBC2) DOUBLE COMPLEX FLC1(0:IFLC1,2),FRC1(0:IFRC1,2),FTC1(0:IFTC1,2), & FBC1(0:IFBC1,2),FLC2(0:IFLC2,2),FRC2(0:IFRC2,2),FTC2(0:IFTC2,2), & FBC2(0:IFBC2,2) DOUBLE PRECISION XMESH(0:NMESH) C .. C .. Subroutine Arguments .. EXTERNAL SL4COF,SL4BCS C .. C .. Local Scalars .. DOUBLE COMPLEX ELOC,ANSLC1,ANSLC2,ANSRC1, & ANSRC2,ANSTC1,ANSTC2,ANSBC1,ANSBC2 INTEGER I,NZLC1,NZLC2,NZRC1,NZRC2,NZTC1,NZTC2,NZBC1, & NZBC2,IFLAG1,IFLAG2,IEIGS1,IEIGS2,INDEX DOUBLE PRECISION RLENT1,RLENT2 C .. C WRITE (9,FMT=*) '****** NMESH,IMATCH:',NMESH,IMATCH IF (IMATCH.LT.0.OR.IMATCH.GT.NMESH) THEN IMATCH = NMESH/2 END IF C WRITE (9,FMT=*) '****** NMESH,IMATCH:',NMESH,IMATCH IF (IEIGS(1).LE.0) THEN IFAIL = MAX(1,IFAIL) C WRITE (9,FMT=*) 'Return from ALLEVS' RETURN END IF ELOC = (ZL(0,1)+ZR(0,1))/2.D0 IF (ABS(ZL(0,1)-ZR(0,1)).LT.TOL*ABS(ELOC)) THEN C The diagonal of the box has length less than TOL times the C midpoint. In this case we assign a cluster of IEIGS(1) eigenvalues C all to be equal to the midpoint. DO 10 I = 1,IEIGS(1) EIGVAL(I) = ELOC 10 CONTINUE C WRITE (9,FMT=*) 'Return from ALLEVS' RETURN END IF C If we are here then the box has diameter greater than we C want and the box contains at least one eigenvalue. IF (IEIGS(1).EQ.1) THEN C We have found a box containing just one eigenvalue, which C we now proceed to locate. IFLAG1 = 0 CALL RFNBOX(ZL,IZL,NZL(1),ZR,IZR,NZR(1),ZT,IZT,NZT(1),ZB,IZB, & NZB(1),FL,IFL,FR,IFR,FT,IFT,FB,IFB,ANSL(1), & ANSR(1),ANST(1),ANSB(1),ZLC1,IZLC1,NZLC1,ZRC1, & IZRC1,NZRC1,ZTC1,IZTC1,NZTC1,ZBC1,IZBC1,NZBC1, & FLC1,IFLC1,FRC1,IFRC1,FTC1,IFTC1,FBC1,IFBC1, & ANSLC1,ANSRC1,ANSTC1,ANSBC1,XMESH,NMESH,IMATCH, & SL4COF,SL4BCS,TOL,'C',EIGVAL(1),IFLAG1) IFAIL = MAX(IFAIL,IFLAG1) IF (IFLAG1.NE.0) THEN C If the error is not too dire we return the eigenvalue approximations C anyway. IF (ABS(ZL(0,1)-ZR(0,1)).LT.SQRT(TOL)*ABS(ELOC)) THEN EIGVAL(1) = 0.5D0*(ZL(0,1)+ZR(0,1)) END IF END IF C WRITE (9,FMT=*) 'Eigenvalue located,IFLAG1, IEIGS(1)=',IFLAG1, C & IEIGS(1) C WRITE (9,FMT=*) 'Return from ALLEVS' RETURN C ELSE C We form appropriate child boxes and use recursion RLENT1 = MIN(ABS(DBLE(ZL(NZL(1),1))),ABS(DBLE(ZR(0,1)))) RLENT1 = MAX(1.D0,SQRT(RLENT1)) RLENT2 = MIN(ABS(AIMAG(ZT(0,1))),ABS(AIMAG(ZB(NZB(1),1)))) RLENT2 = MAX(1.D0,SQRT(RLENT2)) C IF (ABS(ZL(NZL(1),1)-ZR(0,1))/RLENT1.GT. C & ABS(ZT(0,1)-ZB(NZB(1),1))/RLENT2) THEN IF (ABS(ZL(NZL(1),1)-ZR(0,1)).GT.ABS(ZT(0,1)-ZB(NZB(1),1))) THEN C We have a wide, short box. Form the left and right children. C WRITE (9,FMT=*) 'IMATCH,NMESH,IEIGS(1) before LCHILD:',IMATCH, C & NMESH,IEIGS(1) IFLAG1 = 0 CALL LCHILD(ZL(0,1),IZL,NZL(1),ZR(0,1),IZR,NZR(1),ZT(0,1),IZT, & NZT(1),ZB(0,1),IZB,NZB(1),FL(0,1,1),IFL,FR(0,1,1), & IFR,FT(0,1,1),IFT,FB(0,1,1),IFB,ANSL(1),ANSR(1), & ANST(1),ANSB(1),ZLC1,IZLC1,NZLC1,ZRC1,IZRC1,NZRC1, & ZTC1,IZTC1,NZTC1,ZBC1,IZBC1,NZBC1,FLC1,IFLC1,FRC1, & IFRC1,FTC1,IFTC1,FBC1,IFBC1,ANSLC1,ANSRC1,ANSTC1, & ANSBC1,XMESH,NMESH,IMATCH, & SL4COF,SL4BCS,IFLAG1) IF (IFLAG1.NE.0) THEN IFAIL = MAX(IFAIL,IFLAG1) C WRITE (9,FMT=*) 'Return from ALLEVS' RETURN END IF C WRITE (9,FMT=*) 'IMATCH,NMESH,IEIGS(1) after LCHILD:',IMATCH, C & NMESH,IEIGS(1) IEIGS1 = NINT(AIMAG(ANSLC1+ANSRC1+ANSTC1+ANSBC1)/ & (8.D0*ATAN(1.D0))) C WRITE (9,FMT=*) 'LCHILD FOUND',IEIGS1,' EIGENVALUES' IF (IEIGS1.EQ.IEIGS(1)) THEN C All the eigenvalues are in the left box. Copy the child box C arrays onto the parent box arrays and continue. IF (NZLC1.GT.IZL.OR.NZRC1.GT.IZR.OR.NZTC1.GT.IZT. & OR.NZBC1.GT.IZB) THEN IFAIL = MAX(25,IFAIL) C WRITE (9,FMT=*) 'Return from ALLEVS' RETURN END IF DO 23 I = 0,NZLC1 ZL(I,1) = ZLC1(I) FL(I,1,1) = FLC1(I,1) FL(I,2,1) = FLC1(I,2) 23 CONTINUE ANSL(1) = ANSLC1 NZL(1) = NZLC1 DO 33 I = 0,NZRC1 ZR(I,1) = ZRC1(I) FR(I,1,1) = FRC1(I,1) FR(I,2,1) = FRC1(I,2) 33 CONTINUE ANSR(1) = ANSRC1 NZR(1) = NZRC1 DO 43 I = 0,NZTC1 ZT(I,1) = ZTC1(I) FT(I,1,1) = FTC1(I,1) FT(I,2,1) = FTC1(I,2) 43 CONTINUE ANST(1) = ANSTC1 NZT(1) = NZTC1 DO 53 I = 0,NZBC1 ZB(I,1) = ZBC1(I) FB(I,1,1) = FBC1(I,1) FB(I,2,1) = FBC1(I,2) 53 CONTINUE ANSB(1) = ANSBC1 NZB(1) = NZBC1 C WRITE (9,FMT=*) 'Calling ALLEVS with IEIGS(1)=',IEIGS(1) CALL ALLEVS(ZL(0,1),IZL,NZL(1),ZR(0,1),IZR,NZR(1),ZT(0,1),IZT, & NZT(1),ZB(0,1),IZB,NZB(1),FL(0,1,1),IFL,FR(0,1,1), & IFR,FT(0,1,1),IFT,FB(0,1,1),IFB,ANSL(1),ANSR(1), & ANST(1),ANSB(1),ZLC1,IZLC1,ZRC1,IZRC1,ZTC1,IZTC1, & ZBC1,IZBC1,FLC1,IFLC1,FRC1,IFRC1,FTC1,IFTC1,FBC1, & IFBC1,ZLC2,IZLC2,ZRC2,IZRC2,ZTC2,IZTC2,ZBC2,IZBC2, & FLC2,IFLC2,FRC2,IFRC2,FTC2,IFTC2,FBC2,IFBC2, & XMESH,NMESH,IMATCH,SL4COF,SL4BCS,TOL, & EIGVAL,NEIGS,IEIGS,IFAIL) C WRITE (9,FMT=*) 'Return from ALLEVS with IEIGS(1)=',IEIGS(1) RETURN END IF C Form right child: C WRITE (9,FMT=*) 'IEIGS(1) before RCHILD:',IEIGS(1) IFLAG1 = 0 CALL RCHILD(ZL(0,1),IZL,NZL(1),ZR(0,1),IZR,NZR(1),ZT(0,1),IZT, & NZT(1),ZB(0,1),IZB,NZB(1),FL(0,1,1),IFL,FR(0,1,1), & IFR,FT(0,1,1),IFT,FB(0,1,1),IFB,ANSL(1),ANSR(1), & ANST(1),ANSB(1),ZLC2,IZLC2,NZLC2,ZRC2,IZRC2,NZRC2, & ZTC2,IZTC2,NZTC2,ZBC2,IZBC2,NZBC2,FLC2,IFLC2,FRC2, & IFRC2,FTC2,IFTC2,FBC2,IFBC2,ANSLC2,ANSRC2,ANSTC2, & ANSBC2,XMESH,NMESH,IMATCH, & SL4COF,SL4BCS,IFLAG1) C WRITE (9,FMT=*) 'IMATCH,NMESH,IEIGS(1),IFLAG1 after RCHILD:', C & IMATCH,NMESH,IEIGS(1),IFLAG1 IF (IFLAG1.NE.0) THEN IFAIL = MAX(IFAIL,IFLAG1) C WRITE (9,FMT=*) 'Return from ALLEVS' RETURN END IF IEIGS2 = NINT(AIMAG(ANSLC2+ANSRC2+ANSTC2+ANSBC2)/ & (8.D0*ATAN(1.D0))) C WRITE (9,FMT=*) 'RCHILD FOUND',IEIGS2,' EIGENVALUES' C WRITE (9,FMT=*) '-* IEIGS1,IEIGS2,NEIGS:',IEIGS1,IEIGS2,IEIGS(1) IF (IEIGS1+IEIGS2.NE.IEIGS(1)) THEN C The eigenvalue counts do not match up. Return error flag. C Eventually we shall have to work out what to do in this C case. C WRITE (9,FMT=*) '* IEIGS1,IEIGS2,NEIGS:',IEIGS1,IEIGS2, C & IEIGS(1) IFAIL = MAX(24,IFAIL) C WRITE (9,FMT=*) 'Return from ALLEVS' RETURN END IF ELSE C We have a long, thin box. Form the top and bottom children. C Form top child: C WRITE (9,FMT=*) 'IEIGS(1) before TCHILD:',IEIGS(1) IFLAG1 = 0 CALL TCHILD(ZL(0,1),IZL,NZL(1),ZR(0,1),IZR,NZR(1),ZT(0,1),IZT, & NZT(1),ZB(0,1),IZB,NZB(1),FL(0,1,1),IFL,FR(0,1,1), & IFR,FT(0,1,1),IFT,FB(0,1,1),IFB,ANSL(1),ANSR(1), & ANST(1),ANSB(1),ZLC1,IZLC1,NZLC1,ZRC1,IZRC1,NZRC1, & ZTC1,IZTC1,NZTC1,ZBC1,IZBC1,NZBC1,FLC1,IFLC1,FRC1, & IFRC1,FTC1,IFTC1,FBC1,IFBC1,ANSLC1,ANSRC1,ANSTC1, & ANSBC1,XMESH,NMESH,IMATCH, & SL4COF,SL4BCS,IFLAG1) C WRITE (9,FMT=*) 'NMESH,IMATCH,IEIGS(1) after TCHILD:',NMESH, C & IMATCH,IEIGS(1) IF (IFLAG1.NE.0) THEN IFAIL = MAX(IFAIL,IFLAG1) C WRITE (9,FMT=*) 'Return from ALLEVS' RETURN END IF IEIGS1 = NINT(AIMAG(ANSLC1+ANSRC1+ANSTC1+ANSBC1)/ & (8.D0*ATAN(1.D0))) C WRITE (9,FMT=*) 'TCHILD FOUND',IEIGS1,' EIGENVALUES' IF (IEIGS1.EQ.IEIGS(1)) THEN C All the eigenvalues are in the top box. Copy the child box C arrays onto the parent box arrays and continue. IF (NZLC1.GT.IZL.OR.NZRC1.GT.IZR.OR.NZTC1.GT.IZT. & OR.NZBC1.GT.IZB) THEN IFAIL = MAX(25,IFAIL) C WRITE (9,FMT=*) 'Return from ALLEVS' RETURN END IF DO 20 I = 0,NZLC1 ZL(I,1) = ZLC1(I) FL(I,1,1) = FLC1(I,1) FL(I,2,1) = FLC1(I,2) 20 CONTINUE ANSL(1) = ANSLC1 NZL(1) = NZLC1 DO 30 I = 0,NZRC1 ZR(I,1) = ZRC1(I) FR(I,1,1) = FRC1(I,1) FR(I,2,1) = FRC1(I,2) 30 CONTINUE ANSR(1) = ANSRC1 NZR(1) = NZRC1 DO 40 I = 0,NZTC1 ZT(I,1) = ZTC1(I) FT(I,1,1) = FTC1(I,1) FT(I,2,1) = FTC1(I,2) 40 CONTINUE ANST(1) = ANSTC1 NZT(1) = NZTC1 DO 50 I = 0,NZBC1 ZB(I,1) = ZBC1(I) FB(I,1,1) = FBC1(I,1) FB(I,2,1) = FBC1(I,2) 50 CONTINUE ANSB(1) = ANSBC1 NZB(1) = NZBC1 C WRITE (9,FMT=*) 'Current box:' C WRITE (9,FMT=*) ZT(0,1),ZL(0,1),ZB(0,1),ZR(0,1) C WRITE (9,FMT=*) 'Mesh sizes (T,L,B,R):',NZT(1), C & NZL(1),NZB(1),NZR(1) C WRITE (9,FMT=*) '* Calling ALLEVS with IEIGS(1)=',IEIGS(1) CALL ALLEVS(ZL(0,1),IZL,NZL(1),ZR(0,1),IZR,NZR(1),ZT(0,1),IZT, & NZT(1),ZB(0,1),IZB,NZB(1),FL(0,1,1),IFL,FR(0,1,1), & IFR,FT(0,1,1),IFT,FB(0,1,1),IFB,ANSL(1),ANSR(1), & ANST(1),ANSB(1),ZLC1,IZLC1,ZRC1,IZRC1,ZTC1,IZTC1, & ZBC1,IZBC1,FLC1,IFLC1,FRC1,IFRC1,FTC1,IFTC1,FBC1, & IFBC1,ZLC2,IZLC2,ZRC2,IZRC2,ZTC2,IZTC2,ZBC2,IZBC2, & FLC2,IFLC2,FRC2,IFRC2,FTC2,IFTC2,FBC2,IFBC2, & XMESH,NMESH,IMATCH,SL4COF,SL4BCS,TOL, & EIGVAL,NEIGS,IEIGS,IFAIL) C WRITE (9,FMT=*) 'Return from ALLEVS with IEIGS(1)=',IEIGS(1) RETURN END IF C Form bottom child: C WRITE (9,FMT=*) 'Before BCHILD, IEIGS(1)=',IEIGS(1) IFLAG1 = 0 CALL BCHILD(ZL(0,1),IZL,NZL(1),ZR(0,1),IZR,NZR(1),ZT(0,1),IZT, & NZT(1),ZB(0,1),IZB,NZB(1),FL(0,1,1),IFL,FR(0,1,1), & IFR,FT(0,1,1),IFT,FB(0,1,1),IFB,ANSL(1),ANSR(1), & ANST(1),ANSB(1),ZLC2,IZLC2,NZLC2,ZRC2,IZRC2,NZRC2, & ZTC2,IZTC2,NZTC2,ZBC2,IZBC2,NZBC2,FLC2,IFLC2,FRC2, & IFRC2,FTC2,IFTC2,FBC2,IFBC2,ANSLC2,ANSRC2,ANSTC2, & ANSBC2,XMESH,NMESH,IMATCH, & SL4COF,SL4BCS,IFLAG1) C WRITE (9,FMT=*) 'IMATCH,IEIGS(1) after BCHILD:',IMATCH,IEIGS(1) IF (IFLAG1.NE.0) THEN IFAIL = MAX(IFAIL,IFLAG1) C WRITE (9,FMT=*) 'Return from ALLEVS' RETURN END IF C WRITE (9,FMT=*) 'IMATCH after BCHILD:',IMATCH IEIGS2 = NINT(AIMAG(ANSLC2+ANSRC2+ANSTC2+ANSBC2)/ & (8.D0*ATAN(1.D0))) C WRITE (9,FMT=*) 'BCHILD FOUND',IEIGS2,' EIGENVALUES' IF (IEIGS1+IEIGS2.NE.IEIGS(1)) THEN C The eigenvalue counts do not match up. Return error flag. C Eventually we shall have to work out what to do in this C case. C WRITE (9,FMT=*) 'IEIGS1,IEIGS2,NEIGS:',IEIGS1,IEIGS2, C & IEIGS(1) IFAIL = MAX(23,IFAIL) C WRITE (9,FMT=*) 'Return from ALLEVS' RETURN END IF END IF C Have now formed the child boxes. C Copy the child arrays into appropriate parts of the parent C arrays. IF (NZLC1.GT.IZL.OR.NZRC1.GT.IZR.OR.NZTC1.GT.IZT. & OR.NZBC1.GT.IZB) THEN IFAIL = MAX(25,IFAIL) C WRITE (9,FMT=*) 'Return from ALLEVS' RETURN END IF IF (NZLC2.GT.IZL.OR.NZRC2.GT.IZR.OR.NZTC2.GT.IZT. & OR.NZBC2.GT.IZB) THEN IFAIL = MAX(25,IFAIL) C WRITE (9,FMT=*) 'Return from ALLEVS' RETURN END IF IF (IEIGS1.NE.0) THEN DO 21 I = 0,NZLC1 ZL(I,1) = ZLC1(I) FL(I,1,1) = FLC1(I,1) FL(I,2,1) = FLC1(I,2) 21 CONTINUE ANSL(1) = ANSLC1 NZL(1) = NZLC1 DO 31 I = 0,NZRC1 ZR(I,1) = ZRC1(I) FR(I,1,1) = FRC1(I,1) FR(I,2,1) = FRC1(I,2) 31 CONTINUE ANSR(1) = ANSRC1 NZR(1) = NZRC1 DO 41 I = 0,NZTC1 ZT(I,1) = ZTC1(I) FT(I,1,1) = FTC1(I,1) FT(I,2,1) = FTC1(I,2) 41 CONTINUE ANST(1) = ANSTC1 NZT(1) = NZTC1 DO 51 I = 0,NZBC1 ZB(I,1) = ZBC1(I) FB(I,1,1) = FBC1(I,1) FB(I,2,1) = FBC1(I,2) 51 CONTINUE ANSB(1) = ANSBC1 NZB(1) = NZBC1 END IF INDEX = IEIGS1 + 1 DO 22 I = 0,NZLC2 ZL(I,INDEX) = ZLC2(I) FL(I,1,INDEX) = FLC2(I,1) FL(I,2,INDEX) = FLC2(I,2) 22 CONTINUE ANSL(INDEX) = ANSLC2 NZL(INDEX) = NZLC2 DO 32 I = 0,NZRC2 ZR(I,INDEX) = ZRC2(I) FR(I,1,INDEX) = FRC2(I,1) FR(I,2,INDEX) = FRC2(I,2) 32 CONTINUE ANSR(INDEX) = ANSRC2 NZR(INDEX) = NZRC2 DO 42 I = 0,NZTC2 ZT(I,INDEX) = ZTC2(I) FT(I,1,INDEX) = FTC2(I,1) FT(I,2,INDEX) = FTC2(I,2) 42 CONTINUE ANST(INDEX) = ANSTC2 NZT(INDEX) = NZTC2 DO 52 I = 0,NZBC2 ZB(I,INDEX) = ZBC2(I) FB(I,1,INDEX) = FBC2(I,1) FB(I,2,INDEX) = FBC2(I,2) 52 CONTINUE ANSB(INDEX) = ANSBC2 NZB(INDEX) = NZBC2 C Very important that the next statement be executed before C any recursive calls to ALLEVS IEIGS(INDEX) = IEIGS2 IF (IEIGS1.NE.0) THEN IEIGS(1) = IEIGS1 C WRITE (13,FMT=*) 'INDEX before ALLEVS:',INDEX C IFLAG2 = 0 C WRITE (9,FMT=*) '*** Calling ALLEVS with IEIGS(INDEX)=', C & IEIGS(INDEX) CALL ALLEVS(ZL(0,INDEX),IZL,NZL(INDEX),ZR(0,INDEX),IZR, & NZR(INDEX),ZT(0,INDEX),IZT,NZT(INDEX),ZB(0,INDEX), & IZB,NZB(INDEX),FL(0,1,INDEX),IFL,FR(0,1,INDEX),IFR, & FT(0,1,INDEX),IFT,FB(0,1,INDEX),IFB,ANSL(INDEX), & ANSR(INDEX),ANST(INDEX),ANSB(INDEX),ZLC1,IZLC1, & ZRC1,IZRC1,ZTC1,IZTC1,ZBC1,IZBC1, & FLC1,IFLC1,FRC1,IFRC1,FTC1,IFTC1,FBC1,IFBC1, & ZLC2,IZLC2,ZRC2,IZRC2,ZTC2,IZTC2,ZBC2,IZBC2, & FLC2,IFLC2,FRC2,IFRC2,FTC2,IFTC2,FBC2,IFBC2, & XMESH,NMESH,IMATCH,SL4COF,SL4BCS,TOL, & EIGVAL(INDEX),NEIGS,IEIGS(INDEX),IFAIL) C IF (IFLAG2.NE.0) THEN C WRITE (9,FMT=*) 'IFLAG2 was ',IFLAG2 C END IF C IFLAG1 = 0 C WRITE (9,FMT=*) '** Calling ALLEVS with IEIGS(1)=',IEIGS(1) C WRITE (9,FMT=*) 'Check IMATCH before in:',IMATCH CALL ALLEVS(ZL(0,1),IZL,NZL(1),ZR(0,1),IZR,NZR(1),ZT(0,1),IZT, & NZT(1),ZB(0,1),IZB,NZB(1),FL(0,1,1),IFL,FR(0,1,1), & IFR,FT(0,1,1),IFT,FB(0,1,1),IFB,ANSL(1),ANSR(1), & ANST(1),ANSB(1),ZLC1,IZLC1,ZRC1,IZRC1,ZTC1,IZTC1, & ZBC1,IZBC1,FLC1,IFLC1,FRC1,IFRC1,FTC1,IFTC1,FBC1, & IFBC1,ZLC2,IZLC2,ZRC2,IZRC2,ZTC2,IZTC2,ZBC2,IZBC2, & FLC2,IFLC2,FRC2,IFRC2,FTC2,IFTC2,FBC2,IFBC2, & XMESH,NMESH,IMATCH,SL4COF,SL4BCS,TOL, & EIGVAL,NEIGS,IEIGS,IFAIL) C IF (IFLAG1.NE.0) THEN C WRITE (9,FMT=*) 'IFLAG1 was ',IFLAG1 C END IF C IFAIL = MAX(IFLAG1,IFLAG2) C WRITE (9,FMT=*) 'Return from ALLEVS' RETURN ELSE C WRITE (13,FMT=*) 'INDEX before ALLEVS:',INDEX C IFLAG2 = 0 C WRITE (9,FMT=*) '*** Calling ALLEVS with IEIGS(INDEX)=', C & IEIGS(INDEX) CALL ALLEVS(ZL(0,INDEX),IZL,NZL(INDEX),ZR(0,INDEX),IZR, & NZR(INDEX),ZT(0,INDEX),IZT,NZT(INDEX),ZB(0,INDEX), & IZB,NZB(INDEX),FL(0,1,INDEX),IFL,FR(0,1,INDEX),IFR, & FT(0,1,INDEX),IFT,FB(0,1,INDEX),IFB,ANSL(INDEX), & ANSR(INDEX),ANST(INDEX),ANSB(INDEX),ZLC1,IZLC1, & ZRC1,IZRC1,ZTC1,IZTC1,ZBC1,IZBC1, & FLC1,IFLC1,FRC1,IFRC1,FTC1,IFTC1,FBC1,IFBC1, & ZLC2,IZLC2,ZRC2,IZRC2,ZTC2,IZTC2,ZBC2,IZBC2, & FLC2,IFLC2,FRC2,IFRC2,FTC2,IFTC2,FBC2,IFBC2, & XMESH,NMESH,IMATCH,SL4COF,SL4BCS,TOL, & EIGVAL(INDEX),NEIGS,IEIGS(INDEX),IFAIL) C IF (IFLAG2.NE.0) THEN C WRITE (9,FMT=*) 'IFLAG2 was ',IFLAG2 C END IF C IFAIL = MAX(IFAIL,IFLAG2) RETURN END IF END IF END C ----------------------------------------------------------------------- SUBROUTINE STRIP(XMESH,NMESH,IMATCH,SL4COF,SL4BCS,SLSMPC,SLSMBC, & RELAM,EPS,K,TOL,RLOW,RUP,DL,DU,IFAIL) C Subroutine to get a strip containing the `K-th'eigenvalue of a C non-selfadjoint SL problem C C XMESH(0:NMESH): input. This array specifies a mesh covering C the interval on which the problem is defined. C C NMESH: the number of intervals in the mesh. C C IMATCH: the matchpoint index, between 0 and NMESH. C C SL4COF: external SUBROUTINE defining the coefficients in C the differential equation. C C SL4BCS: external SUBROUTINE defining the boundary conditions C for the problem. C C SLSMPC: external SUBROUTINE defining the CONSTANT coefficients in C a simple differential equation used for asymptotics / comparison C purposes. C C SLSMBC: external SUBROUTINE defining the boundary conditions C for the simple comparison problem. C C RELAM: initial estimate of the real part of the eigenvalue(s) C sought. Need not be at all accurate. C C EPS: estimate of error in RELAM. Thus, the routine will look C initially for eigenvalues in the strip C C RELAM - EPS < Re(Lambda) < RELAM + EPS C C but this strip will be expanded or shifted if necessary C to locate eigenvalues. The final strip in which the C eigenvalue(s) is(are) located is C C RLOW < Re(Lambda) < RUP C C (see below). C C K: the routine will attempt to ensure that there are precisely C K eigenvalues with Re(Lambda) < RLOW and precisely K+1 C with Re(Lambda) < RUP. However for some values of K this C may not be possible because of the presence of eigenvalues C of multiplicity greater than 1. In such cases, the routine C will try to narrow the strip until RUP-RLOW < TOL (see below C for description of TOL). C C TOL: In case the eigenvalue sought turns out to have real part C very close to the real part of another eigenvalue, TOL C specifies the width of the strip required around the C eigenvalue. Must be set on entry; unchanged on exit. C C RLOW, RUP: If IFAIL = 0 or IFAIL > 9 on exit, then RLOW and C RUP specify a strip RLOW < Re(Lambda) < RUP containing C at least one eigenvalue of the problem. C C DL, DU: IF IFAIL = 0 or IFAIL > 9 on exit, then C C DL = -0.5 + no. of eigenvalues with real part < RLOW C DU = -0.5 + no. of eigenvalues with real part < RUP C C and hence DU - DL is the number of eigenvalues in the C strip. C C IFAIL: the error flag. Should be set to 0 on entry. On exit: C C IFAIL = 0: successful exit. RLOW and RUP specify a strip containing C an eigenvalue and DU - DL specifies the number of eigenvalues in C this strip. C IFAIL between 1 and 7: failure in the routine NZERO. See NZERO C for the meaning of this value if IFAIL. C IFAIL = 8: the routine was able to determine RLOW but not RUP. C IFAIL = 9: the routine was able to determine RUP but not RLOW. C IFAIL = 10: the strip RLOW < Re(Lambda) < RUP contains more C than one eigenvalue but the code has not managed to narrow C its width to be less than TOL, due to an internal restriction C on the number of calls to NZERO. C IFAIL = 13: the routine located RLOW and RUP with RUP-RLOW < TOL C but the strip appears to contain several eigenvalues. C C .. Parameters .. INTEGER N,LDFJAC,LR PARAMETER (N=2,LDFJAC=N,LR=N* (N+1)/2+1) C .. C .. Scalar Arguments .. DOUBLE PRECISION DL,DU,EPS,RELAM,RLOW,RUP,TOL INTEGER IFAIL,IMATCH,K,NMESH C .. C .. Local Scalars .. DOUBLE COMPLEX FMIN,FUN,ZMIN DOUBLE PRECISION CONST,EPSFCN,FACTOR,RELOC,W1,W2,XTOL INTEGER IFLAG,ILOW,IREVCM,ISECT,IUP,ML,MODE,MU,NZ,IFEVAL LOGICAL DONELO,DONEUP,NONVAL C .. C .. Local Arrays .. DOUBLE COMPLEX ZL(2),ZMAX(2),ZU(2) DOUBLE PRECISION ARGQH(2),ARGQL(2),ARGQU(2),DIAG(N), & FJAC(LDFJAC,N),FVEC(N),QTF(N),R(LR),W(N,4),X(N) C .. C .. External Functions .. DOUBLE PRECISION X02AJF EXTERNAL X02AJF C .. C .. External Subroutines .. EXTERNAL C05NDF,NZERO,SLMISS C .. C .. Intrinsic Functions .. INTRINSIC ABS,CMPLX,DBLE,AIMAG,MAX,REAL,SQRT C .. C .. Array Arguments .. DOUBLE PRECISION XMESH(0:NMESH) C .. C .. Subroutine Arguments .. EXTERNAL SL4BCS,SL4COF,SLSMBC,SLSMPC C .. C First use bisection to find a strip a < Re(lambda) < b containing C just one eigenvalue, if this is possible. DONEUP = .FALSE. DONELO = .FALSE. IUP = 0 ILOW = 0 RELOC = RELAM IFLAG = 0 CALL NZERO(XMESH,NMESH,IMATCH,RELOC,SL4COF,SL4BCS,SLSMPC, & SLSMBC,NZ,ZMIN,FMIN,ZMAX,ARGQH,IFLAG) IF (IFLAG.NE.0) THEN IFAIL = IFLAG RETURN END IF IF (DBLE(NZ-K)-0.5D0.LT.0.D0) THEN C RELOC is less than Re(lambda) RLOW = RELOC DL = DBLE(NZ-K) - 0.5D0 ARGQL(1) = ARGQH(1) ARGQL(2) = ARGQH(2) ZL(1) = ZMAX(1) ZL(2) = ZMAX(2) DONELO = .TRUE. ELSE RUP = RELOC DU = DBLE(NZ-K) - 0.5D0 ARGQU(1) = ARGQH(1) ARGQU(2) = ARGQH(2) ZU(1) = ZMAX(1) ZU(2) = ZMAX(2) DONEUP = .TRUE. END IF IF (.NOT.DONEUP) THEN C Look for upper bound on eigenvalue RELOC = RELAM + EPS 10 IFLAG = 0 CALL NZERO(XMESH,NMESH,IMATCH,RELOC,SL4COF,SL4BCS,SLSMPC, & SLSMBC,NZ,ZMIN,FMIN,ZMAX,ARGQH,IFLAG) IF (IFLAG.NE.0) THEN IFAIL = IFLAG RETURN END IF IF (DBLE(NZ-K)-0.5D0.LT.0.D0) THEN C Upper bound not found; tighten lower bound RLOW = RELOC DL = DBLE(NZ-K) - 0.5D0 ARGQL(1) = ARGQH(1) ARGQL(2) = ARGQH(2) ZL(1) = ZMAX(1) ZL(2) = ZMAX(2) EPS = 2.D0*EPS RELOC = RELOC + EPS IUP = IUP + 1 IF (IUP.GT.40) THEN IFAIL = 8 RETURN END IF GO TO 10 END IF C If we are here then we have found an upper bound DU = DBLE(NZ-K) - 0.5D0 RUP = RELOC ARGQU(1) = ARGQH(1) ARGQU(2) = ARGQH(2) ZU(1) = ZMAX(1) ZU(2) = ZMAX(2) DONEUP = .TRUE. ELSE C IF (.NOT.DONELO) THEN C Look for lower bound on eigenvalue RELOC = RELAM - EPS 20 IFLAG = 0 CALL NZERO(XMESH,NMESH,IMATCH,RELOC,SL4COF,SL4BCS,SLSMPC, & SLSMBC,NZ,ZMIN,FMIN,ZMAX,ARGQH,IFLAG) IF (IFLAG.NE.0) THEN IFAIL = IFLAG RETURN END IF IF (DBLE(NZ-K)-0.5D0.GT.0.D0) THEN C Lower bound not found; tighten upper bound RUP = RELOC DU = DBLE(NZ-K) - 0.5D0 ARGQU(1) = ARGQH(1) ARGQU(2) = ARGQH(2) ZU(1) = ZMAX(1) ZU(2) = ZMAX(2) EPS = 2.D0*EPS RELOC = RELOC - EPS ILOW = ILOW + 1 IF (ILOW.GT.40) THEN IFAIL = 9 RETURN END IF GO TO 20 END IF C If we are here then we have found a lower bound DL = DBLE(NZ-K) - 0.5D0 RLOW = RELOC ARGQL(1) = ARGQH(1) ARGQL(2) = ARGQH(2) ZL(1) = ZMAX(1) ZL(2) = ZMAX(2) DONELO = .TRUE. END IF C We now have both upper and lower bounds on the eigenvalue. ISECT = 0 C NONVAL = ARGQL(1) .EQ. 0.D0 .OR. ARGQU(2) .EQ. 0.D0 C 30 IF ((NONVAL.OR.ABS(DU-DL).GT.1.5D0) .AND. C & ABS(RUP-RLOW).GT.TOL) THEN 30 IF (ABS(DU-DL).GT.1.5D0 .AND. & ABS(RUP-RLOW).GT.TOL) THEN C We want to narrow down the strip containing the eigenvalue RELOC = (ABS(DL)*RUP+ABS(DU)*RLOW)/ (ABS(DL)+ABS(DU)) IFLAG = 0 CALL NZERO(XMESH,NMESH,IMATCH,RELOC,SL4COF,SL4BCS,SLSMPC, & SLSMBC,NZ,ZMIN,FMIN,ZMAX,ARGQH,IFLAG) IF (IFLAG.NE.0) THEN IFAIL = IFLAG RETURN END IF IF (DBLE(NZ-K)-0.5D0.GT.0.D0) THEN RUP = RELOC DU = DBLE(NZ-K) - 0.5D0 ARGQU(1) = ARGQH(1) ARGQU(2) = ARGQH(2) ZU(1) = ZMAX(1) ZU(2) = ZMAX(2) ELSE RLOW = RELOC DL = DBLE(NZ-K) - 0.5D0 ARGQL(1) = ARGQH(1) ARGQL(2) = ARGQH(2) ZL(1) = ZMAX(1) ZL(2) = ZMAX(2) END IF C NONVAL = ARGQL(1) .EQ. 0.D0 .OR. ARGQU(2) .EQ. 0.D0 ISECT = ISECT + 1 IF (ISECT.GT.40) THEN IFAIL = 10 RETURN END IF GO TO 30 END IF C If we are here then we have either found a strip of width C less than TOL containing the eigenvalue, plus maybe some other C eigenvalues as well, or else we have found a strip containing C just one eigenvalue. IF (ABS(DU-DL).GT.1.5D0) THEN IFAIL = 13 ELSE IFAIL = 0 END IF RETURN END C ----------------------------------------------------------------------- SUBROUTINE FNDBOX(XL,XR,ICOUNT,IDIRCT,EBR,ETR,ETL, & EBL,ZL,IZL,NZL,ZR,IZR,NZR,ZT,IZT,NZT,ZB, & IZB,NZB,FL,IFL,FR,IFR,FT,IFT,FB,IFB,ANSL, & ANSR,ANST,ANSB,ANS,XMESH,NMESH,IMATCH, & SL4COF,SL4BCS,IFAIL) INTEGER ICOUNT,IDIRCT,IFAIL,IMATCH,NMESH,IFB,IFL,IFR,IFT, & IZB,IZL,IZR,IZT,NZB,NZL,NZR,NZT DOUBLE PRECISION XL,XR,YB,YBOLD,YT DOUBLE PRECISION XMESH(0:NMESH) DOUBLE COMPLEX ANSB,ANSL,ANSR,ANST,ANS,EBL,EBR,ETL,ETR DOUBLE COMPLEX FB(0:IFB,1:2),FL(0:IFL,1:2),FR(0:IFR,1:2), & FT(0:IFT,1:2),ZB(0:IZB),ZL(0:IZL),ZR(0:IZR), & ZT(0:IZT) EXTERNAL MKBOX,SL4BCS,SL4COF,X01AAF DOUBLE PRECISION X01AAF INTRINSIC CMPLX,AIMAG,NINT YB = 0.D0 YT = SQRT(X01AAF(0.D0)/2.D0) IDIRCT = 1 ICOUNT = 0 ANS = CMPLX(0.D0,0.D0) 30 EBR = CMPLX(XR,YB) ETR = CMPLX(XR,YT) ETL = CMPLX(XL,YT) EBL = CMPLX(XL,YB) ICOUNT = ICOUNT + 1 CALL MKBOX(EBR,ETR,ETL,EBL,ZL,IZL,NZL,ZR,IZR,NZR,ZT,IZT,NZT, & ZB,IZB,NZB,FL,IFL,FR,IFR,FT,IFT,FB,IFB,ANSL, & ANSR,ANST,ANSB,ANS,XMESH,NMESH,IMATCH,SL4COF,SL4BCS, & IFAIL) WRITE(6,FMT=*) 'ICOUNT =' ,ICOUNT WRITE(6,FMT=*) 'ANS =' ,ANS WRITE(6,FMT=*) 'EBR =' ,EBR WRITE(6,FMT=*) 'ETR =' ,ETR WRITE(6,FMT=*) 'ETL =' ,ETL WRITE(6,FMT=*) 'EBL =' ,EBL IF (NINT(AIMAG(ANS)).EQ.0) THEN IF (IDIRCT.EQ.1) THEN YBOLD = YB YB = -YT YT = -YBOLD IDIRCT = 2 GOTO 30 ELSE YBOLD = YB YB = -YB YT = YB + 2.D0*(YT - YBOLD) IDIRCT = 1 GOTO 30 ENDIF ELSE WRITE(6,FMT=*) 'The current box has corners:' WRITE(6,FMT=*) 'EBR = ' ,EBR WRITE(6,FMT=*) 'ETR = ' ,ETR WRITE(6,FMT=*) 'ETL = ' ,ETL WRITE(6,FMT=*) 'EBL = ' ,EBL WRITE(6,FMT=*) 'The number of eigenvalues in the box is' & ,NINT(AIMAG(ANS)) ENDIF RETURN END C --------------------------------------------------------------------- SUBROUTINE RFNBOX(ZL,IZL,NZL,ZR,IZR,NZR,ZT,IZT,NZT,ZB,IZB,NZB,FL, & IFL,FR,IFR,FT,IFT,FB,IFB,ANSL,ANSR,ANST,ANSB, & ZLC,IZLC,NZLC,ZRC,IZRC,NZRC,ZTC,IZTC,NZTC,ZBC, & IZBC,NZBC,FLC,IFLC,FRC,IFRC,FTC,IFTC,FBC,IFBC, & ANSLC,ANSRC,ANSTC,ANSBC,XMESH,NMESH,IMATCH, & SL4COF,SL4BCS,TOL,MODE,EIGVAL,IFAIL) C C This routine refines a box until one of the following is C true: C C IF (MODE.EQ.'B') THEN the code bisects the box until both sides C are of length less than TOL. C C IF (MODE.EQ.'C') THEN the code bisects the box until it appears C to be safe to use a linear approximation to the miss-distance C function to find the eigenvalue. The eigenvalue is located C to a tolerance TOL. C C C Get `length' of box C .. Scalar Arguments .. DOUBLE COMPLEX ANSB,ANSBC,ANSL,ANSLC,ANSR,ANSRC,ANST,ANSTC,EIGVAL DOUBLE PRECISION TOL INTEGER IFAIL,IFB,IFBC,IFL,IFLC,IFR,IFRC,IFT,IFTC,IMATCH,IZB,IZBC, & IZL,IZLC,IZR,IZRC,IZT,IZTC,NMESH,NZB,NZBC,NZL,NZLC,NZR, & NZRC,NZT,NZTC CHARACTER MODE C .. C .. Array Arguments .. DOUBLE COMPLEX FB(0:IFB,1:2),FBC(0:IFBC,1:2),FL(0:IFL,1:2), & FLC(0:IFLC,1:2),FR(0:IFR,1:2),FRC(0:IFRC,1:2), & FT(0:IFT,1:2),FTC(0:IFTC,1:2),ZB(0:IZB), & ZBC(0:IZBC),ZL(0:IZL),ZLC(0:IZLC),ZR(0:IZR), & ZRC(0:IZRC),ZT(0:IZT),ZTC(0:IZTC) DOUBLE PRECISION XMESH(0:NMESH) C .. C .. Subroutine Arguments .. EXTERNAL SL4BCS,SL4COF C .. C .. Local Scalars .. DOUBLE COMPLEX ELAM,F,FD1,FD2,ZNEW,ZNEW1,ZNEW2 DOUBLE PRECISION CMAX,CONST,ERROR,LENGTH INTEGER I,IEVAL LOGICAL BISECT,NONLIN,OUT C .. C .. External Subroutines .. EXTERNAL NXTBOX,SLMISS C .. C .. Intrinsic Functions .. INTRINSIC ABS,DBLE,AIMAG,EXP,MAX,MIN C .. IEVAL = 0 10 LENGTH=MAX( & ABS(ZL(0)-ZB(0))/ & MAX(1.D0,MIN(ABS(AIMAG(ZL(0))),ABS(AIMAG(ZB(0))))), & ABS(ZB(0)-ZR(0))/ & MAX(1.D0,MIN(ABS(DBLE(ZR(0))),ABS(DBLE(ZB(0)))))) C C WRITE (9,FMT=*) 'Bottom right corner:',ZR(0) C WRITE (9,FMT=*) 'Top right corner:',ZT(0) C WRITE (9,FMT=*) 'Top left corner:',ZL(0) C WRITE (9,FMT=*) 'Bottom left corner:',ZB(0) C Decide whether or not miss-distance is well-approximated C by a linear function. CALL SLMISS(XMESH,NMESH,IMATCH,ZL(0),SL4COF,SL4BCS,F,CONST,IFAIL) C WRITE (7,FMT=*) 'F,FL(0,1),CONST,FL(0,2):',F,FL(0,1),CONST,FL(0,2) CALL SLMISS(XMESH,NMESH,IMATCH,ZR(0),SL4COF,SL4BCS,F,CONST,IFAIL) C WRITE (7,FMT=*) 'F,FR(0,1),CONST,FR(0,2):',F,FR(0,1),CONST,FR(0,2) CALL SLMISS(XMESH,NMESH,IMATCH,ZT(0),SL4COF,SL4BCS,F,CONST,IFAIL) C WRITE (7,FMT=*) 'F,FT(0,1),CONST,FT(0,2):',F,FT(0,1),CONST,FT(0,2) CALL SLMISS(XMESH,NMESH,IMATCH,ZB(0),SL4COF,SL4BCS,F,CONST,IFAIL) C WRITE (7,FMT=*) 'F,FB(0,1),CONST,FB(0,2):',F,FB(0,1),CONST,FB(0,2) CMAX = MAX(DBLE(FL(0,2)),DBLE(FR(0,2))) CMAX = MAX(CMAX,DBLE(FT(0,2))) CMAX = MAX(CMAX,DBLE(FB(0,2))) FD1 = (FL(0,1)*EXP(DBLE(FL(0,2))-CMAX)- & FR(0,1)*EXP(DBLE(FR(0,2))-CMAX))/ (ZL(0)-ZR(0)) FD2 = (FT(0,1)*EXP(DBLE(FT(0,2))-CMAX)- & FB(0,1)*EXP(DBLE(FB(0,2))-CMAX))/ (ZT(0)-ZB(0)) ZNEW1 = ZL(0) - FL(0,1)*EXP(DBLE(FL(0,2))-CMAX)/FD1 ZNEW2 = ZT(0) - FT(0,1)*EXP(DBLE(FT(0,2))-CMAX)/FD2 NONLIN = ABS(FD1-FD2) .GT. 0.1D0*MIN(ABS(FD1),ABS(FD2)) C WRITE (9,FMT=*) 'Initial NONLIN:',NONLIN OUT = MAX(DBLE(ZNEW1),DBLE(ZNEW2)) .GT. DBLE(ZR(0)) + TOL OUT = OUT .OR. MIN(DBLE(ZNEW1),DBLE(ZNEW2)) .LT. DBLE(ZL(0)) - TOL OUT = OUT .OR. MAX(AIMAG(ZNEW1),AIMAG(ZNEW2)) .GT. & AIMAG(ZT(0)) + TOL OUT = OUT .OR. MIN(AIMAG(ZNEW1),AIMAG(ZNEW2)) .LT. & AIMAG(ZB(0)) - TOL C WRITE (9,FMT=*) 'OUT:',OUT C WRITE (7,FMT=*) 'ZR(0),ZT(0),ZL(0),ZB(0):',ZR(0),ZT(0),ZL(0),ZB(0) C WRITE (7,FMT=*) 'ZR(NZR),ZT(NZT),ZL(NZL),ZB(NZB):',ZR(NZR), C & ZT(NZT),ZL(NZL),ZB(NZB) NONLIN = NONLIN .OR. OUT C WRITE (9,FMT=*) 'FD1,FD2:',FD1,FD2 C WRITE (9,FMT=*) 'ABS(FD1-FD2),MIN(ABS(FD1),ABS(FD2)):', C & ABS(FD1-FD2),MIN(ABS(FD1),ABS(FD2)) C WRITE (9,FMT=*) 'ZNEW1,ZNEW2:',ZNEW1,ZNEW2 C NONLIN = .TRUE. BISECT = (LENGTH.GT.TOL .AND. MODE.EQ.'B') .OR. & (NONLIN .AND. MODE.EQ.'C' .AND. LENGTH.GT.TOL) C WRITE (9,FMT=*) 'MODE,NONLIN,BISECT:',MODE,NONLIN,BISECT IF (.NOT.BISECT) THEN GO TO 70 END IF C We do the bisection of the box. 20 CALL NXTBOX(ZL,IZL,NZL,ZR,IZR,NZR,ZT,IZT,NZT,ZB,IZB,NZB,FL,IFL,FR, & IFR,FT,IFT,FB,IFB,ANSL,ANSR,ANST,ANSB,ZLC,IZLC,NZLC, & ZRC,IZRC,NZRC,ZTC,IZTC,NZTC,ZBC,IZBC,NZBC,FLC,IFLC, & FRC,IFRC,FTC,IFTC,FBC,IFBC,ANSLC,ANSRC,ANSTC,ANSBC, & XMESH,NMESH,IMATCH,SL4COF,SL4BCS,IFAIL) IF (IFAIL.NE.0) RETURN C Copy new box onto old box NZL = NZLC NZR = NZRC NZT = NZTC NZB = NZBC ANSL = ANSLC ANSR = ANSRC ANST = ANSTC ANSB = ANSBC DO 30 I = 0,NZL ZL(I) = ZLC(I) FL(I,1) = FLC(I,1) FL(I,2) = FLC(I,2) 30 CONTINUE DO 40 I = 0,NZR ZR(I) = ZRC(I) FR(I,1) = FRC(I,1) FR(I,2) = FRC(I,2) 40 CONTINUE DO 50 I = 0,NZT ZT(I) = ZTC(I) FT(I,1) = FTC(I,1) FT(I,2) = FTC(I,2) 50 CONTINUE DO 60 I = 0,NZB ZB(I) = ZBC(I) FB(I,1) = FBC(I,1) FB(I,2) = FBC(I,2) 60 CONTINUE C This trap deals with the case where bisection has C been re-entered after a failed attempt to use Newton's C method. IF (IEVAL.GT.5) THEN BISECT = .TRUE. IEVAL = IEVAL - 1 GO TO 20 END IF IEVAL = 0 GO TO 10 70 CONTINUE C If we are here then we have finished bisection. IF (MODE.EQ.'B' .OR. LENGTH.LT.TOL) THEN EIGVAL = (ZL(0)+ZR(0)+ZT(0)+ZB(0))/4.D0 IFAIL = 0 IF (MODE.EQ.'C' .AND. LENGTH.LT.TOL) THEN C WRITE (9,FMT=*) 'Bisected to convergence' END IF RETURN END IF C If we are here then we should find the eigenvalue by Newton's C method. IEVAL = 0 ELAM = 0.5D0* (ZNEW1+ZNEW2) 80 CALL SLMISS(XMESH,NMESH,IMATCH,ELAM,SL4COF,SL4BCS,F,CONST,IFAIL) IF (IFAIL.NE.0) RETURN ZNEW = ELAM - 2.D0*F*EXP(CONST-CMAX)/ (FD1+FD2) C WRITE (9,FMT=*) 'Best approx so far:',ZNEW C WRITE (9,FMT=*) 'IEVAL=',IEVAL IF (ABS(ZNEW-ELAM).GT.TOL .AND. IEVAL.LT.20) THEN IEVAL = IEVAL + 1 ELAM = ZNEW IF (DBLE(ELAM).GT.DBLE(ZR(0))+TOL .OR. & DBLE(ELAM).LT.DBLE(ZL(0))-TOL .OR. & AIMAG(ELAM).GT.AIMAG(ZT(0))+TOL .OR. & AIMAG(ELAM).LT.AIMAG(ZB(0))-TOL) THEN BISECT = .TRUE. GO TO 20 END IF GO TO 80 END IF ERROR = ABS(DBLE(ZNEW-ELAM))/ & MAX(1.D0,MIN(ABS(DBLE(ELAM)),ABS(DBLE(ZNEW)))) & + ABS(AIMAG(ZNEW-ELAM))/ & MAX(1.D0,MIN(ABS(AIMAG(ELAM)),ABS(AIMAG(ZNEW)))) IF (ERROR.GT.TOL .AND. IEVAL.GE.5) THEN C EIGVAL = ZNEW C IFAIL = 22 BISECT = .TRUE. IEVAL = IEVAL - 1 C WRITE (9,FMT=*) 'Return to bisection' C WRITE (8,FMT=*) 'Return to bisection' C WRITE (8,FMT=*) 'NZL,NZR,NZT,NZB:',NZL,NZR,NZT,NZB C WRITE (8,FMT=*) 'ZL(NZL),ZB(0):',ZL(NZL),ZB(0) C WRITE (8,FMT=*) 'ZB(NZB),ZR(0):',ZB(NZB),ZR(0) C WRITE (8,FMT=*) 'ZR(NZR),ZT(0):',ZR(NZR),ZT(0) C WRITE (8,FMT=*) 'ZT(NZT),ZL(0):',ZT(NZT),ZL(0) GO TO 20 C RETURN END IF C If we are here then we have found an eigenvalue EIGVAL = ZNEW IFAIL = 0 RETURN END C --------------------------------------------------------------------- SUBROUTINE NXTBOX(ZL,IZL,NZL,ZR,IZR,NZR,ZT,IZT,NZT,ZB,IZB,NZB,FL, & IFL,FR,IFR,FT,IFT,FB,IFB,ANSL,ANSR,ANST,ANSB, & ZLC,IZLC,NZLC,ZRC,IZRC,NZRC,ZTC,IZTC,NZTC,ZBC, & IZBC,NZBC,FLC,IFLC,FRC,IFRC,FTC,IFTC,FBC,IFBC, & ANSLC,ANSRC,ANSTC,ANSBC,XMESH,NMESH,IMATCH, & SL4COF,SL4BCS,IFAIL) C Given a box containing one eigenvalue, this routine bisects the box C in one of two directions, and forms a new box containing one C eigenvalue. C C First decide which way to bisect C .. Scalar Arguments .. DOUBLE COMPLEX ANSB,ANSBC,ANSL,ANSLC,ANSR,ANSRC,ANST,ANSTC INTEGER IFAIL,IFB,IFBC,IFL,IFLC,IFR,IFRC,IFT,IFTC,IMATCH,IZB,IZBC, & IZL,IZLC,IZR,IZRC,IZT,IZTC,NMESH,NZB,NZBC,NZL,NZLC,NZR, & NZRC,NZT,NZTC C .. C .. Array Arguments .. DOUBLE COMPLEX FB(0:IFB,1:2),FBC(0:IFBC,1:2),FL(0:IFL,1:2), & FLC(0:IFLC,1:2),FR(0:IFR,1:2),FRC(0:IFRC,1:2), & FT(0:IFT,1:2),FTC(0:IFTC,1:2),ZB(0:IZB), & ZBC(0:IZBC),ZL(0:IZL),ZLC(0:IZLC),ZR(0:IZR), & ZRC(0:IZRC),ZT(0:IZT),ZTC(0:IZTC) DOUBLE PRECISION XMESH(0:NMESH) C .. C .. Subroutine Arguments .. EXTERNAL SL4BCS,SL4COF C .. C .. Local Scalars .. INTEGER IEIGS DOUBLE PRECISION RLENT1,RLENT2 C .. C .. External Subroutines .. EXTERNAL BCHILD,LCHILD,LRCHLD,RCHILD,TCHILD,TBCHLD C .. C .. Intrinsic Functions .. INTRINSIC ABS,ATAN,AIMAG,NINT C .. RLENT1 = MIN(ABS(DBLE(ZL(NZL))),ABS(DBLE(ZR(0)))) RLENT1 = MAX(1.D0,SQRT(RLENT1)) RLENT2 = MIN(ABS(AIMAG(ZT(0))),ABS(AIMAG(ZB(NZB)))) RLENT2 = MAX(1.D0,SQRT(RLENT2)) IF (ABS(ZL(NZL)-ZR(0))/RLENT1.GT. & ABS(ZT(0)-ZB(NZB))/RLENT2) THEN C We have a wide, short box. Form the left and right children. IFAIL = 0 CALL LCHILD(ZL,IZL,NZL,ZR,IZR,NZR,ZT,IZT,NZT,ZB,IZB,NZB,FL, & IFL,FR,IFR,FT,IFT,FB,IFB,ANSL,ANSR,ANST,ANSB,ZLC, & IZLC,NZLC,ZRC,IZRC,NZRC,ZTC,IZTC,NZTC,ZBC,IZBC, & NZBC,FLC,IFLC,FRC,IFRC,FTC,IFTC,FBC,IFBC,ANSLC, & ANSRC,ANSTC,ANSBC,XMESH,NMESH,IMATCH,SL4COF, & SL4BCS,IFAIL) IF (IFAIL.NE.0) RETURN IEIGS = NINT(AIMAG(ANSLC+ANSRC+ANSTC+ANSBC)/ & (8.D0*ATAN(1.D0))) C WRITE (9,FMT=*) 'ANSLC,ANSRC,ANSTC,ANSBC:', C & ANSLC,ANSRC,ANSTC,ANSBC IF (IEIGS.GT.0) THEN C The new box is the left child. GO TO 10 ELSE C The new box is the right child IFAIL = 0 CALL RCHILD(ZL,IZL,NZL,ZR,IZR,NZR,ZT,IZT,NZT,ZB,IZB,NZB, & FL,IFL,FR,IFR,FT,IFT,FB,IFB,ANSL,ANSR,ANST, & ANSB,ZLC,IZLC,NZLC,ZRC,IZRC,NZRC,ZTC,IZTC, & NZTC,ZBC,IZBC,NZBC,FLC,IFLC,FRC,IFRC,FTC,IFTC, & FBC,IFBC,ANSLC,ANSRC,ANSTC,ANSBC,XMESH,NMESH, & IMATCH,SL4COF,SL4BCS,IFAIL) IF (IFAIL.NE.0) RETURN IEIGS = NINT(AIMAG(ANSLC+ANSRC+ANSTC+ANSBC)/ & (8.D0*ATAN(1.D0))) C WRITE (9,FMT=*) 'ANSLC,ANSRC,ANSTC,ANSBC:', C & ANSLC,ANSRC,ANSTC,ANSBC IF (IEIGS.LE.0) THEN C Try doing an LRCHLD: IFAIL = 0 CALL LRCHLD(ZL,IZL,NZL,ZR,IZR,NZR,ZT,IZT,NZT,ZB,IZB,NZB, & FL,IFL,FR,IFR,FT,IFT,FB,IFB,ANSL,ANSR,ANST, & ANSB,ZLC,IZLC,NZLC,ZRC,IZRC,NZRC,ZTC,IZTC, & NZTC,ZBC,IZBC,NZBC,FLC,IFLC,FRC,IFRC,FTC,IFTC, & FBC,IFBC,ANSLC,ANSRC,ANSTC,ANSBC,XMESH,NMESH, & IMATCH,SL4COF,SL4BCS,IFAIL) IF (IFAIL.NE.0) RETURN IEIGS = NINT(AIMAG(ANSLC+ANSRC+ANSTC+ANSBC)/ & (8.D0*ATAN(1.D0))) C WRITE (9,FMT=*) 'ANSLC,ANSRC,ANSTC,ANSBC:', C & ANSLC,ANSRC,ANSTC,ANSBC IF (IEIGS.LE.0) THEN IFAIL = 20 RETURN END IF END IF END IF ELSE C We have a long, thin box. Form the top and bottom children. IFAIL = 0 CALL TCHILD(ZL,IZL,NZL,ZR,IZR,NZR,ZT,IZT,NZT,ZB,IZB,NZB,FL, & IFL,FR,IFR,FT,IFT,FB,IFB,ANSL,ANSR,ANST,ANSB,ZLC, & IZLC,NZLC,ZRC,IZRC,NZRC,ZTC,IZTC,NZTC,ZBC,IZBC, & NZBC,FLC,IFLC,FRC,IFRC,FTC,IFTC,FBC,IFBC,ANSLC, & ANSRC,ANSTC,ANSBC,XMESH,NMESH,IMATCH,SL4COF, & SL4BCS,IFAIL) IF (IFAIL.NE.0) RETURN IEIGS = NINT(AIMAG(ANSLC+ANSRC+ANSTC+ANSBC)/ & (8.D0*ATAN(1.D0))) C WRITE (9,FMT=*) 'ANSLC,ANSRC,ANSTC,ANSBC:', C & ANSLC,ANSRC,ANSTC,ANSBC IF (IEIGS.GT.0) THEN C The new box is the left child. GO TO 10 ELSE C The new box is the right child IFAIL = 0 CALL BCHILD(ZL,IZL,NZL,ZR,IZR,NZR,ZT,IZT,NZT,ZB,IZB,NZB, & FL,IFL,FR,IFR,FT,IFT,FB,IFB,ANSL,ANSR,ANST, & ANSB,ZLC,IZLC,NZLC,ZRC,IZRC,NZRC,ZTC,IZTC, & NZTC,ZBC,IZBC,NZBC,FLC,IFLC,FRC,IFRC,FTC,IFTC, & FBC,IFBC,ANSLC,ANSRC,ANSTC,ANSBC,XMESH,NMESH, & IMATCH,SL4COF,SL4BCS,IFAIL) IF (IFAIL.NE.0) RETURN IEIGS = NINT(AIMAG(ANSLC+ANSRC+ANSTC+ANSBC)/ & (8.D0*ATAN(1.D0))) C WRITE (9,FMT=*) 'ANSLC,ANSRC,ANSTC,ANSBC:', C & ANSLC,ANSRC,ANSTC,ANSBC IF (IEIGS.LE.0) THEN C Try doing a TBCHLD: IFAIL = 0 CALL TBCHLD(ZL,IZL,NZL,ZR,IZR,NZR,ZT,IZT,NZT,ZB,IZB,NZB, & FL,IFL,FR,IFR,FT,IFT,FB,IFB,ANSL,ANSR,ANST, & ANSB,ZLC,IZLC,NZLC,ZRC,IZRC,NZRC,ZTC,IZTC, & NZTC,ZBC,IZBC,NZBC,FLC,IFLC,FRC,IFRC,FTC,IFTC, & FBC,IFBC,ANSLC,ANSRC,ANSTC,ANSBC,XMESH,NMESH, & IMATCH,SL4COF,SL4BCS,IFAIL) IF (IFAIL.NE.0) RETURN IEIGS = NINT(AIMAG(ANSLC+ANSRC+ANSTC+ANSBC)/ & (8.D0*ATAN(1.D0))) C WRITE (9,FMT=*) 'ANSLC,ANSRC,ANSTC,ANSBC:', C & ANSLC,ANSRC,ANSTC,ANSBC IF (IEIGS.LE.0) THEN IFAIL = 21 RETURN END IF END IF END IF END IF 10 IFAIL = 0 RETURN END C --------------------------------------------------------------------- SUBROUTINE LCHILD(ZL,IZL,NZL,ZR,IZR,NZR,ZT,IZT,NZT,ZB,IZB,NZB,FL, & IFL,FR,IFR,FT,IFT,FB,IFB,ANSL,ANSR,ANST,ANSB, & ZLC,IZLC,NZLC,ZRC,IZRC,NZRC,ZTC,IZTC,NZTC,ZBC, & IZBC,NZBC,FLC,IFLC,FRC,IFRC,FTC,IFTC,FBC,IFBC, & ANSLC,ANSRC,ANSTC,ANSBC,XMESH,NMESH,IMATCH, & SL4COF,SL4BCS,IFAIL) C Given a box specified by ZL,ZR,ZT,ZB,FL,FR,FT,FB,ANSL,ANSR,ANST, C ANSB, this routine computes the arrays for the child-box which C is the left half of the input box. These arrays are denoted C by ZLC,ZRC,ZTC,ZBC,FLC,FRC,FTC,FBC, and the corresponding integrals C are denoted ANSLC,ANSRC,ANSTC,ANSBC. C C As this routine generates the left-half child, the left-half C arrays are unchanged. C .. Scalar Arguments .. DOUBLE COMPLEX ANSB,ANSBC,ANSL,ANSLC,ANSR,ANSRC,ANST,ANSTC INTEGER IFAIL,IFB,IFBC,IFL,IFLC,IFR,IFRC,IFT,IFTC,IMATCH,IZB,IZBC, & IZL,IZLC,IZR,IZRC,IZT,IZTC,NMESH,NZB,NZBC,NZL,NZLC,NZR, & NZRC,NZT,NZTC C .. C .. Array Arguments .. DOUBLE COMPLEX FB(0:IFB,1:2),FBC(0:IFBC,1:2),FL(0:IFL,1:2), & FLC(0:IFLC,1:2),FR(0:IFR,1:2),FRC(0:IFRC,1:2), & FT(0:IFT,1:2),FTC(0:IFTC,1:2),ZB(0:IZB), & ZBC(0:IZBC),ZL(0:IZL),ZLC(0:IZLC),ZR(0:IZR), & ZRC(0:IZRC),ZT(0:IZT),ZTC(0:IZTC) DOUBLE PRECISION XMESH(0:NMESH) C .. C .. Subroutine Arguments .. EXTERNAL SL4BCS,SL4COF C .. C .. Local Scalars .. DOUBLE COMPLEX EBR,ELAM,ETR,F DOUBLE PRECISION CONST INTEGER I,IFLAG,IND,IREVCM,J C .. C .. External Subroutines .. EXTERNAL BXSIDE,SLMISS C .. C .. Intrinsic Functions .. INTRINSIC CMPLX,DBLE,LOG,MIN C .. IFAIL = 0 IF (MIN(IZLC,IFLC).LT.NZL) THEN IFAIL = 18 RETURN END IF NZLC = NZL DO 10 I = 0,NZL ZLC(I) = ZL(I) FLC(I,1) = FL(I,1) FLC(I,2) = FL(I,2) ANSLC = ANSL 10 CONTINUE C That does the left. Now do the right, which has to be generated C from scratch. IREVCM = 1 IFLAG = 0 C Get top and bottom endpoints of the right side of the child box: C these are, respectively, the midpoints of the top and bottom C sides of the parent box. ETR = 0.5D0* (ZT(0)+ZT(NZT)) EBR = 0.5D0* (ZB(0)+ZB(NZB)) C Do the integral from EBR to ETR: 20 CALL BXSIDE(EBR,ETR,ELAM,F,CONST,IREVCM,ZRC,IZRC,NZRC,FRC,IFRC, & ANSRC,IFLAG) IF (IREVCM.NE.0) THEN C Function value required: SLMISS will put it in the variable F, C with logarithmic scaling in CONST. CALL SLMISS(XMESH,NMESH,IMATCH,ELAM,SL4COF,SL4BCS,F,CONST, & IFLAG) IF (IFLAG.NE.0) THEN IFAIL = IFLAG RETURN END IF GO TO 20 END IF IF (IFLAG.NE.0) THEN IFAIL = IFLAG RETURN END IF C We have now done the sides. We must do the top and bottom. C Do bottom. I = 0 30 IF (DBLE(ZB(I)).LT.DBLE(EBR)) THEN ZBC(I) = ZB(I) FBC(I,1) = FB(I,1) FBC(I,2) = FB(I,2) I = I + 1 NZBC = I GO TO 30 END IF ZBC(NZBC) = ZRC(0) FBC(NZBC,1) = FRC(0,1) FBC(NZBC,2) = FRC(0,2) ANSBC = CMPLX(0.D0,0.D0) DO 40 I = 1,NZBC ANSBC = ANSBC + LOG(FBC(I,1)/FBC(I-1,1)) 40 CONTINUE C Do top. I = 0 50 IF (DBLE(ZT(I)).GE.DBLE(ETR)) THEN I = I + 1 GO TO 50 END IF ZTC(0) = ZRC(NZRC) FTC(0,1) = FRC(NZRC,1) FTC(0,2) = FRC(NZRC,2) NZTC = NZT - I + 1 DO 60 J = 1,NZTC IND = J + I - 1 ZTC(J) = ZT(IND) FTC(J,1) = FT(IND,1) FTC(J,2) = FT(IND,2) 60 CONTINUE ANSTC = CMPLX(0.D0,0.D0) DO 70 I = 1,NZTC ANSTC = ANSTC + LOG(FTC(I,1)/FTC(I-1,1)) 70 CONTINUE C All sides are now done. RETURN END C --------------------------------------------------------------------- SUBROUTINE RCHILD(ZL,IZL,NZL,ZR,IZR,NZR,ZT,IZT,NZT,ZB,IZB,NZB,FL, & IFL,FR,IFR,FT,IFT,FB,IFB,ANSL,ANSR,ANST,ANSB, & ZLC,IZLC,NZLC,ZRC,IZRC,NZRC,ZTC,IZTC,NZTC,ZBC, & IZBC,NZBC,FLC,IFLC,FRC,IFRC,FTC,IFTC,FBC,IFBC, & ANSLC,ANSRC,ANSTC,ANSBC,XMESH,NMESH,IMATCH, & SL4COF,SL4BCS,IFAIL) C Given a box specified by ZL,ZR,ZT,ZB,FL,FR,FT,FB,ANSL,ANSR,ANST, C ANSB, this routine computes the arrays for the child-box which C is the right half of the input box. These arrays are denoted C by ZLC,ZRC,ZTC,ZBC,FLC,FRC,FTC,FBC, and the corresponding integrals C are denoted ANSLC,ANSRC,ANSTC,ANSBC. C C As this routine generates the right-half child, the right-half C arrays are unchanged. C .. Scalar Arguments .. DOUBLE COMPLEX ANSB,ANSBC,ANSL,ANSLC,ANSR,ANSRC,ANST,ANSTC INTEGER IFAIL,IFB,IFBC,IFL,IFLC,IFR,IFRC,IFT,IFTC,IMATCH,IZB,IZBC, & IZL,IZLC,IZR,IZRC,IZT,IZTC,NMESH,NZB,NZBC,NZL,NZLC,NZR, & NZRC,NZT,NZTC C .. C .. Array Arguments .. DOUBLE COMPLEX FB(0:IFB,1:2),FBC(0:IFBC,1:2),FL(0:IFL,1:2), & FLC(0:IFLC,1:2),FR(0:IFR,1:2),FRC(0:IFRC,1:2), & FT(0:IFT,1:2),FTC(0:IFTC,1:2),ZB(0:IZB), & ZBC(0:IZBC),ZL(0:IZL),ZLC(0:IZLC),ZR(0:IZR), & ZRC(0:IZRC),ZT(0:IZT),ZTC(0:IZTC) DOUBLE PRECISION XMESH(0:NMESH) C .. C .. Subroutine Arguments .. EXTERNAL SL4BCS,SL4COF C .. C .. Local Scalars .. DOUBLE COMPLEX EBL,ELAM,ETL,F DOUBLE PRECISION CONST INTEGER I,IFLAG,IND,IREVCM,J C .. C .. External Subroutines .. EXTERNAL BXSIDE,SLMISS C .. C .. Intrinsic Functions .. INTRINSIC CMPLX,DBLE,LOG,MIN C .. IFAIL = 0 IF (MIN(IZRC,IFRC).LT.NZR) THEN IFAIL = 18 RETURN END IF NZRC = NZR DO 10 I = 0,NZR ZRC(I) = ZR(I) FRC(I,1) = FR(I,1) FRC(I,2) = FR(I,2) ANSRC = ANSR 10 CONTINUE C That does the right. Now do the left, which has to be generated C from scratch. IREVCM = 1 IFLAG = 0 C Get top and bottom endpoints of the left side of the child box: C these are, respectively, the midpoints of the top and bottom C sides of the parent box. ETL = 0.5D0* (ZT(0)+ZT(NZT)) EBL = 0.5D0* (ZB(0)+ZB(NZB)) C Do the integral from ETL to EBL: 20 CALL BXSIDE(ETL,EBL,ELAM,F,CONST,IREVCM,ZLC,IZLC,NZLC,FLC,IFLC, & ANSLC,IFLAG) IF (IREVCM.NE.0) THEN C Function value required: SLMISS will put it in the variable F, C with logarithmic scaling in CONST. CALL SLMISS(XMESH,NMESH,IMATCH,ELAM,SL4COF,SL4BCS,F,CONST, & IFLAG) IF (IFLAG.NE.0) THEN IFAIL = IFLAG RETURN END IF GO TO 20 END IF IF (IFLAG.NE.0) THEN IFAIL = IFLAG RETURN END IF C We have now done the sides. We must do the top and bottom. C Do top. I = 0 30 IF (DBLE(ZT(I)).GT.DBLE(ETL)) THEN ZTC(I) = ZT(I) FTC(I,1) = FT(I,1) FTC(I,2) = FT(I,2) I = I + 1 NZTC = I GO TO 30 END IF ZTC(NZTC) = ZLC(0) FTC(NZTC,1) = FLC(0,1) FTC(NZTC,2) = FLC(0,2) ANSTC = CMPLX(0.D0,0.D0) DO 40 I = 1,NZTC ANSTC = ANSTC + LOG(FTC(I,1)/FTC(I-1,1)) 40 CONTINUE C Do bottom. I = 0 50 IF (DBLE(ZB(I)).LE.DBLE(EBL)) THEN I = I + 1 GO TO 50 END IF ZBC(0) = ZLC(NZLC) FBC(0,1) = FLC(NZLC,1) FBC(0,2) = FLC(NZLC,2) NZBC = NZB - I + 1 DO 60 J = 1,NZBC IND = J + I - 1 ZBC(J) = ZB(IND) FBC(J,1) = FB(IND,1) FBC(J,2) = FB(IND,2) 60 CONTINUE ANSBC = CMPLX(0.D0,0.D0) DO 70 I = 1,NZBC ANSBC = ANSBC + LOG(FBC(I,1)/FBC(I-1,1)) 70 CONTINUE C All sides are now done. RETURN END C --------------------------------------------------------------------- SUBROUTINE LRCHLD(ZL,IZL,NZL,ZR,IZR,NZR,ZT,IZT,NZT,ZB,IZB,NZB,FL, & IFL,FR,IFR,FT,IFT,FB,IFB,ANSL,ANSR,ANST,ANSB, & ZLC,IZLC,NZLC,ZRC,IZRC,NZRC,ZTC,IZTC,NZTC,ZBC, & IZBC,NZBC,FLC,IFLC,FRC,IFRC,FTC,IFTC,FBC,IFBC, & ANSLC,ANSRC,ANSTC,ANSBC,XMESH,NMESH,IMATCH, & SL4COF,SL4BCS,IFAIL) C C Given a box specified by ZL,ZR,ZT,ZB,FL,FR,FT,FB,ANSL,ANSR,ANST, C ANSB, this routine computes the arrays for the child-box which C is the middle vertical half of the input box. These arrays are denoted C by ZLC,ZRC,ZTC,ZBC,FLC,FRC,FTC,FBC, and the corresponding integrals C are denoted ANSLC,ANSRC,ANSTC,ANSBC. C C .. Scalar Arguments .. DOUBLE COMPLEX ANSB,ANSBC,ANSL,ANSLC,ANSR,ANSRC,ANST,ANSTC INTEGER IFAIL,IFB,IFBC,IFL,IFLC,IFR,IFRC,IFT,IFTC,IMATCH,IZB,IZBC, & IZL,IZLC,IZR,IZRC,IZT,IZTC,NMESH,NZB,NZBC,NZL,NZLC,NZR, & NZRC,NZT,NZTC C .. C .. Array Arguments .. DOUBLE COMPLEX FB(0:IFB,1:2),FBC(0:IFBC,1:2),FL(0:IFL,1:2), & FLC(0:IFLC,1:2),FR(0:IFR,1:2),FRC(0:IFRC,1:2), & FT(0:IFT,1:2),FTC(0:IFTC,1:2),ZB(0:IZB), & ZBC(0:IZBC),ZL(0:IZL),ZLC(0:IZLC),ZR(0:IZR), & ZRC(0:IZRC),ZT(0:IZT),ZTC(0:IZTC) DOUBLE PRECISION XMESH(0:NMESH) C .. C .. Subroutine Arguments .. EXTERNAL SL4BCS,SL4COF C .. C .. Local Scalars .. DOUBLE COMPLEX EBR,ELAM,ETR,F,EBL,ETL DOUBLE PRECISION CONST INTEGER I,IFLAG,IND,IREVCM,J C .. C .. External Subroutines .. EXTERNAL BXSIDE,SLMISS C .. C .. Intrinsic Functions .. INTRINSIC CMPLX,DBLE,LOG,MIN C .. IFAIL = 0 IREVCM = 1 IFLAG = 0 C C Get top and bottom endpoints of the right side of the child box: C ETR = 0.75D0*ZT(0)+0.25D0*ZT(NZT) EBR = 0.25D0*ZB(0)+0.75D0*ZB(NZB) C Do the integral from EBR to ETR: 20 CALL BXSIDE(EBR,ETR,ELAM,F,CONST,IREVCM,ZRC,IZRC,NZRC,FRC,IFRC, & ANSRC,IFLAG) IF (IREVCM.NE.0) THEN C Function value required: SLMISS will put it in the variable F, C with logarithmic scaling in CONST. CALL SLMISS(XMESH,NMESH,IMATCH,ELAM,SL4COF,SL4BCS,F,CONST, & IFLAG) IF (IFLAG.NE.0) THEN IFAIL = IFLAG RETURN END IF GO TO 20 END IF IF (IFLAG.NE.0) THEN IFAIL = IFLAG RETURN END IF C Done the right side. Now do the left side. ETL = 0.25D0*ZT(0)+0.75D0*ZT(NZT) EBL = 0.75D0*ZB(0)+0.25D0*ZB(NZB) C Do the integral from ETL to EBL: 21 CALL BXSIDE(ETL,EBL,ELAM,F,CONST,IREVCM,ZLC,IZLC,NZLC,FLC,IFLC, & ANSLC,IFLAG) IF (IREVCM.NE.0) THEN C Function value required: SLMISS will put it in the variable F, C with logarithmic scaling in CONST. CALL SLMISS(XMESH,NMESH,IMATCH,ELAM,SL4COF,SL4BCS,F,CONST, & IFLAG) IF (IFLAG.NE.0) THEN IFAIL = IFLAG RETURN END IF GO TO 21 END IF IF (IFLAG.NE.0) THEN IFAIL = IFLAG RETURN END IF C We have now done the sides. We must do the top and bottom. C Do bottom. I = 0 ZBC(I) = ZLC(NZLC) FBC(I,1) = FLC(NZLC,1) FBC(I,2) = FLC(NZLC,2) I = 1 30 IF (DBLE(ZB(I)).LE.DBLE(ZBC(0))) THEN I = I + 1 GOTO 30 END IF J = 1 31 IF (DBLE(ZB(I)).LT.DBLE(ZRC(0))) THEN ZBC(J) = ZB(I) FBC(J,1) = FB(I,1) FBC(J,2) = FB(I,2) I = I + 1 J = J + 1 NZBC = J GO TO 31 END IF ZBC(NZBC) = ZRC(0) FBC(NZBC,1) = FRC(0,1) FBC(NZBC,2) = FRC(0,2) ANSBC = CMPLX(0.D0,0.D0) DO 40 I = 1,NZBC ANSBC = ANSBC + LOG(FBC(I,1)/FBC(I-1,1)) 40 CONTINUE C Do top. I = 0 ZTC(I) = ZRC(NZRC) FTC(I,1) = FRC(NZRC,1) FTC(I,2) = FRC(NZRC,2) I = 1 50 IF (DBLE(ZT(I)).GE.DBLE(ZRC(NZRC))) THEN I = I + 1 GO TO 50 END IF J = 1 51 IF (DBLE(ZT(I)).GT.DBLE(ZLC(0))) THEN ZTC(J) = ZT(I) FTC(J,1) = FT(I,1) FTC(J,2) = FT(I,2) I = I + 1 J = J + 1 NZTC = J GOTO 51 END IF ZTC(NZTC) = ZLC(0) FTC(NZTC,1) = FLC(0,1) FTC(NZTC,2) = FLC(0,2) ANSTC = CMPLX(0.D0,0.D0) DO 70 I = 1,NZTC ANSTC = ANSTC + LOG(FTC(I,1)/FTC(I-1,1)) 70 CONTINUE C All sides are now done. RETURN END C --------------------------------------------------------------------- SUBROUTINE BCHILD(ZL,IZL,NZL,ZR,IZR,NZR,ZT,IZT,NZT,ZB,IZB,NZB,FL, & IFL,FR,IFR,FT,IFT,FB,IFB,ANSL,ANSR,ANST,ANSB, & ZLC,IZLC,NZLC,ZRC,IZRC,NZRC,ZTC,IZTC,NZTC,ZBC, & IZBC,NZBC,FLC,IFLC,FRC,IFRC,FTC,IFTC,FBC,IFBC, & ANSLC,ANSRC,ANSTC,ANSBC,XMESH,NMESH,IMATCH, & SL4COF,SL4BCS,IFAIL) C Given a box specified by ZL,ZR,ZT,ZB,FL,FR,FT,FB,ANSL,ANSR,ANST, C ANSB, this routine computes the arrays for the child-box which C is the bottom half of the input box. These arrays are denoted C by ZLC,ZRC,ZTC,ZBC,FLC,FRC,FTC,FBC, and the corresponding integrals C are denoted ANSLC,ANSRC,ANSTC,ANSBC. C C As this routine generates the bottom-half child, the bottom-half C arrays are unchanged. C .. Scalar Arguments .. DOUBLE COMPLEX ANSB,ANSBC,ANSL,ANSLC,ANSR,ANSRC,ANST,ANSTC INTEGER IFAIL,IFB,IFBC,IFL,IFLC,IFR,IFRC,IFT,IFTC,IMATCH,IZB,IZBC, & IZL,IZLC,IZR,IZRC,IZT,IZTC,NMESH,NZB,NZBC,NZL,NZLC,NZR, & NZRC,NZT,NZTC C .. C .. Array Arguments .. DOUBLE COMPLEX FB(0:IFB,1:2),FBC(0:IFBC,1:2),FL(0:IFL,1:2), & FLC(0:IFLC,1:2),FR(0:IFR,1:2),FRC(0:IFRC,1:2), & FT(0:IFT,1:2),FTC(0:IFTC,1:2),ZB(0:IZB), & ZBC(0:IZBC),ZL(0:IZL),ZLC(0:IZLC),ZR(0:IZR), & ZRC(0:IZRC),ZT(0:IZT),ZTC(0:IZTC) DOUBLE PRECISION XMESH(0:NMESH) C .. C .. Subroutine Arguments .. EXTERNAL SL4BCS,SL4COF C .. C .. Local Scalars .. DOUBLE COMPLEX ELAM,ETL,ETR,F DOUBLE PRECISION CONST INTEGER I,IFLAG,IND,IREVCM,J C .. C .. External Subroutines .. EXTERNAL BXSIDE,SLMISS C .. C .. Intrinsic Functions .. INTRINSIC CMPLX,AIMAG,LOG,MIN C .. IFAIL = 0 IF (MIN(IZBC,IFBC).LT.NZB) THEN IFAIL = 18 RETURN END IF NZBC = NZB DO 10 I = 0,NZB ZBC(I) = ZB(I) FBC(I,1) = FB(I,1) FBC(I,2) = FB(I,2) ANSBC = ANSB 10 CONTINUE C That does the bottom. Now do the top, which has to be generated C from scratch. IREVCM = 1 IFLAG = 0 C Get left and right endpoints of the top of the child box: C these are, respectively, the midpoints of the left and right C sides of the parent box. ETL = 0.5D0* (ZL(0)+ZL(NZL)) ETR = 0.5D0* (ZR(0)+ZR(NZR)) C Do the integral from ETR to ETL: 20 CALL BXSIDE(ETR,ETL,ELAM,F,CONST,IREVCM,ZTC,IZTC,NZTC,FTC,IFTC, & ANSTC,IFLAG) IF (IREVCM.NE.0) THEN C Function value required: SLMISS will put it in the variable F, C with logarithmic scaling in CONST. CALL SLMISS(XMESH,NMESH,IMATCH,ELAM,SL4COF,SL4BCS,F,CONST, & IFLAG) IF (IFLAG.NE.0) THEN IFAIL = IFLAG RETURN END IF GO TO 20 END IF IF (IFLAG.NE.0) THEN IFAIL = IFLAG RETURN END IF C We have now done top and bottom. We must do the sides. C Do right side. I = 0 30 IF (AIMAG(ZR(I)).LT.AIMAG(ETR)) THEN ZRC(I) = ZR(I) FRC(I,1) = FR(I,1) FRC(I,2) = FR(I,2) I = I + 1 NZRC = I GO TO 30 END IF ZRC(NZRC) = ZTC(0) FRC(NZRC,1) = FTC(0,1) FRC(NZRC,2) = FTC(0,2) ANSRC = CMPLX(0.D0,0.D0) DO 40 I = 1,NZRC ANSRC = ANSRC + LOG(FRC(I,1)/FRC(I-1,1)) 40 CONTINUE C Do left side. I = 0 50 IF (AIMAG(ZL(I)).GE.AIMAG(ETL)) THEN I = I + 1 GO TO 50 END IF ZLC(0) = ZTC(NZTC) FLC(0,1) = FTC(NZTC,1) FLC(0,2) = FTC(NZTC,2) NZLC = NZL - I + 1 DO 60 J = 1,NZLC IND = J + I - 1 ZLC(J) = ZL(IND) FLC(J,1) = FL(IND,1) FLC(J,2) = FL(IND,2) 60 CONTINUE ANSLC = CMPLX(0.D0,0.D0) DO 70 I = 1,NZLC ANSLC = ANSLC + LOG(FLC(I,1)/FLC(I-1,1)) 70 CONTINUE C All sides are now done. RETURN END C --------------------------------------------------------------------- SUBROUTINE TCHILD(ZL,IZL,NZL,ZR,IZR,NZR,ZT,IZT,NZT,ZB,IZB,NZB,FL, & IFL,FR,IFR,FT,IFT,FB,IFB,ANSL,ANSR,ANST,ANSB, & ZLC,IZLC,NZLC,ZRC,IZRC,NZRC,ZTC,IZTC,NZTC,ZBC, & IZBC,NZBC,FLC,IFLC,FRC,IFRC,FTC,IFTC,FBC,IFBC, & ANSLC,ANSRC,ANSTC,ANSBC,XMESH,NMESH,IMATCH, & SL4COF,SL4BCS,IFAIL) C Given a box specified by ZL,ZR,ZT,ZB,FL,FR,FT,FB,ANSL,ANSR,ANST, C ANSB, this routine computes the arrays for the child-box which C is the top half of the input box. These arrays are denoted C by ZLC,ZRC,ZTC,ZBC,FLC,FRC,FTC,FBC, and the corresponding integrals C are denoted ANSLC,ANSRC,ANSTC,ANSBC. C C As this routine generates the top-half child, the top-half C arrays are unchanged. C .. Scalar Arguments .. DOUBLE COMPLEX ANSB,ANSBC,ANSL,ANSLC,ANSR,ANSRC,ANST,ANSTC INTEGER IFAIL,IFB,IFBC,IFL,IFLC,IFR,IFRC,IFT,IFTC,IMATCH,IZB,IZBC, & IZL,IZLC,IZR,IZRC,IZT,IZTC,NMESH,NZB,NZBC,NZL,NZLC,NZR, & NZRC,NZT,NZTC C .. C .. Array Arguments .. DOUBLE COMPLEX FB(0:IFB,1:2),FBC(0:IFBC,1:2),FL(0:IFL,1:2), & FLC(0:IFLC,1:2),FR(0:IFR,1:2),FRC(0:IFRC,1:2), & FT(0:IFT,1:2),FTC(0:IFTC,1:2),ZB(0:IZB), & ZBC(0:IZBC),ZL(0:IZL),ZLC(0:IZLC),ZR(0:IZR), & ZRC(0:IZRC),ZT(0:IZT),ZTC(0:IZTC) DOUBLE PRECISION XMESH(0:NMESH) C .. C .. Subroutine Arguments .. EXTERNAL SL4BCS,SL4COF C .. C .. Local Scalars .. DOUBLE COMPLEX EBL,EBR,ELAM,F DOUBLE PRECISION CONST INTEGER I,IFLAG,IND,IREVCM,J C .. C .. External Subroutines .. EXTERNAL BXSIDE,SLMISS C .. C .. Intrinsic Functions .. INTRINSIC CMPLX,AIMAG,LOG,MIN C .. IFAIL = 0 IF (MIN(IZTC,IFTC).LT.NZT) THEN IFAIL = 18 RETURN END IF NZTC = NZT DO 10 I = 0,NZT ZTC(I) = ZT(I) FTC(I,1) = FT(I,1) FTC(I,2) = FT(I,2) ANSTC = ANST 10 CONTINUE C That does the top. Now do the bottom, which has to be generated C from scratch. IREVCM = 1 IFLAG = 0 C Get left and right endpoints of the bottom of the child box: C these are, respectively, the midpoints of the left and right C sides of the parent box. EBL = 0.5D0* (ZL(0)+ZL(NZL)) EBR = 0.5D0* (ZR(0)+ZR(NZR)) C Do the integral from EBL to EBR: 20 CALL BXSIDE(EBL,EBR,ELAM,F,CONST,IREVCM,ZBC,IZBC,NZBC,FBC,IFBC, & ANSBC,IFLAG) IF (IREVCM.NE.0) THEN C Function value required: SLMISS will put it in the variable F, C with logarithmic scaling in CONST. CALL SLMISS(XMESH,NMESH,IMATCH,ELAM,SL4COF,SL4BCS,F,CONST, & IFLAG) IF (IFLAG.NE.0) THEN IFAIL = IFLAG RETURN END IF GO TO 20 END IF IF (IFLAG.NE.0) THEN IFAIL = IFLAG RETURN END IF C We have now done top and bottom. We must do the sides. C Do left side. I = 0 30 IF (AIMAG(ZL(I)).GT.AIMAG(EBL)) THEN ZLC(I) = ZL(I) FLC(I,1) = FL(I,1) FLC(I,2) = FL(I,2) I = I + 1 NZLC = I GO TO 30 END IF ZLC(NZLC) = ZBC(0) FLC(NZLC,1) = FBC(0,1) FLC(NZLC,2) = FBC(0,2) ANSLC = CMPLX(0.D0,0.D0) DO 40 I = 1,NZLC ANSLC = ANSLC + LOG(FLC(I,1)/FLC(I-1,1)) 40 CONTINUE C Do right side. I = 0 50 IF (AIMAG(ZR(I)).LE.AIMAG(EBR)) THEN I = I + 1 GO TO 50 END IF ZRC(0) = ZBC(NZBC) FRC(0,1) = FBC(NZBC,1) FRC(0,2) = FBC(NZBC,2) NZRC = NZR - I + 1 DO 60 J = 1,NZRC IND = J + I - 1 ZRC(J) = ZR(IND) FRC(J,1) = FR(IND,1) FRC(J,2) = FR(IND,2) 60 CONTINUE ANSRC = CMPLX(0.D0,0.D0) DO 70 I = 1,NZRC ANSRC = ANSRC + LOG(FRC(I,1)/FRC(I-1,1)) 70 CONTINUE C All sides are now done. RETURN END C ----------------------------------------------------------- SUBROUTINE TBCHLD(ZL,IZL,NZL,ZR,IZR,NZR,ZT,IZT,NZT,ZB,IZB,NZB,FL, & IFL,FR,IFR,FT,IFT,FB,IFB,ANSL,ANSR,ANST,ANSB, & ZLC,IZLC,NZLC,ZRC,IZRC,NZRC,ZTC,IZTC,NZTC,ZBC, & IZBC,NZBC,FLC,IFLC,FRC,IFRC,FTC,IFTC,FBC,IFBC, & ANSLC,ANSRC,ANSTC,ANSBC,XMESH,NMESH,IMATCH, & SL4COF,SL4BCS,IFAIL) C Given a box specified by ZL,ZR,ZT,ZB,FL,FR,FT,FB,ANSL,ANSR,ANST, C ANSB, this routine computes the arrays for the child-box which C is the middle horizontal half of the input box. These arrays are denoted C by ZLC,ZRC,ZTC,ZBC,FLC,FRC,FTC,FBC, and the corresponding integrals C are denoted ANSLC,ANSRC,ANSTC,ANSBC. C C .. Scalar Arguments .. DOUBLE COMPLEX ANSB,ANSBC,ANSL,ANSLC,ANSR,ANSRC,ANST,ANSTC INTEGER IFAIL,IFB,IFBC,IFL,IFLC,IFR,IFRC,IFT,IFTC,IMATCH,IZB,IZBC, & IZL,IZLC,IZR,IZRC,IZT,IZTC,NMESH,NZB,NZBC,NZL,NZLC,NZR, & NZRC,NZT,NZTC C .. C .. Array Arguments .. DOUBLE COMPLEX FB(0:IFB,1:2),FBC(0:IFBC,1:2),FL(0:IFL,1:2), & FLC(0:IFLC,1:2),FR(0:IFR,1:2),FRC(0:IFRC,1:2), & FT(0:IFT,1:2),FTC(0:IFTC,1:2),ZB(0:IZB), & ZBC(0:IZBC),ZL(0:IZL),ZLC(0:IZLC),ZR(0:IZR), & ZRC(0:IZRC),ZT(0:IZT),ZTC(0:IZTC) DOUBLE PRECISION XMESH(0:NMESH) C .. C .. Subroutine Arguments .. EXTERNAL SL4BCS,SL4COF C .. C .. Local Scalars .. DOUBLE COMPLEX EBL,EBR,ELAM,F,ETL,ETR DOUBLE PRECISION CONST INTEGER I,IFLAG,IND,IREVCM,J C .. C .. External Subroutines .. EXTERNAL BXSIDE,SLMISS C .. C .. Intrinsic Functions .. INTRINSIC CMPLX,AIMAG,LOG,MIN C .. IFAIL = 0 C C Do the bottom: C IREVCM = 1 IFLAG = 0 C Get left and right endpoints of the bottom of the child box: EBL = 0.25D0*ZL(0)+0.75D0*ZL(NZL) EBR = 0.75D0*ZR(0)+0.25D0*ZR(NZR) C Do the integral from EBL to EBR: 20 CALL BXSIDE(EBL,EBR,ELAM,F,CONST,IREVCM,ZBC,IZBC,NZBC,FBC,IFBC, & ANSBC,IFLAG) IF (IREVCM.NE.0) THEN C Function value required: SLMISS will put it in the variable F, C with logarithmic scaling in CONST. CALL SLMISS(XMESH,NMESH,IMATCH,ELAM,SL4COF,SL4BCS,F,CONST, & IFLAG) IF (IFLAG.NE.0) THEN IFAIL = IFLAG RETURN END IF GO TO 20 END IF IF (IFLAG.NE.0) THEN IFAIL = IFLAG RETURN END IF C C Do the top: C IREVCM = 1 IFLAG = 0 C Get left and right endpoints of the top of the child box: ETL = 0.75D0*ZL(0)+0.25D0*ZL(NZL) ETR = 0.25D0*ZR(0)+0.75D0*ZR(NZR) C Do the integral from ETR to ETL: 21 CALL BXSIDE(ETR,ETL,ELAM,F,CONST,IREVCM,ZTC,IZTC,NZTC,FTC,IFTC, & ANSTC,IFLAG) IF (IREVCM.NE.0) THEN C Function value required: SLMISS will put it in the variable F, C with logarithmic scaling in CONST. CALL SLMISS(XMESH,NMESH,IMATCH,ELAM,SL4COF,SL4BCS,F,CONST, & IFLAG) IF (IFLAG.NE.0) THEN IFAIL = IFLAG RETURN END IF GO TO 21 END IF IF (IFLAG.NE.0) THEN IFAIL = IFLAG RETURN END IF C We have now done top and bottom. We must do the sides. C Do left side. ZLC(0) = ZTC(NZTC) FLC(0,1) = FTC(NZTC,1) FLC(0,2) = FTC(NZTC,2) I = 1 30 IF (AIMAG(ZL(I)).GE.AIMAG(ZTC(NZTC))) THEN I = I + 1 GOTO 30 END IF J = 1 31 IF (AIMAG(ZL(I)).GT.AIMAG(ZBC(0))) THEN ZLC(J) = ZL(I) FLC(J,1) = FL(I,1) FLC(J,2) = FL(I,2) I = I + 1 J = J + 1 GO TO 31 END IF NZLC = J ZLC(NZLC) = ZBC(0) FLC(NZLC,1) = FBC(0,1) FLC(NZLC,2) = FBC(0,2) ANSLC = CMPLX(0.D0,0.D0) DO 40 I = 1,NZLC ANSLC = ANSLC + LOG(FLC(I,1)/FLC(I-1,1)) 40 CONTINUE C Do right side. ZRC(0) = ZBC(NZBC) FRC(0,1) = FBC(NZBC,1) FRC(0,2) = FBC(NZBC,2) I = 1 50 IF (AIMAG(ZR(I)).LE.AIMAG(ZBC(NZBC))) THEN I = I + 1 GO TO 50 END IF J = 1 51 IF (AIMAG(ZR(I)).LT.AIMAG(ZTC(0))) THEN ZRC(J) = ZR(I) FRC(J,1) = FR(I,1) FRC(J,2) = FR(I,2) I = I + 1 J = J + 1 GO TO 51 END IF NZRC = J ZRC(NZRC) = ZTC(0) FRC(NZRC,1) = FTC(0,1) FRC(NZRC,2) = FTC(0,2) ANSRC = CMPLX(0.D0,0.D0) DO 70 I = 1,NZRC ANSRC = ANSRC + LOG(FRC(I,1)/FRC(I-1,1)) 70 CONTINUE C All sides are now done. RETURN END C ----------------------------------------------------------- SUBROUTINE MKBOX(EBR,ETR,ETL,EBL,ZL,IZL,NZL,ZR,IZR,NZR,ZT,IZT,NZT, & ZB,IZB,NZB,FL,IFL,FR,IFR,FT,IFT,FB,IFB,ANSL,ANSR, & ANST,ANSB,ANSWER,XMESH,NMESH,IMATCH,SL4COF, & SL4BCS,IFAIL) C This routine creates a box in the complex plane with corners at C EBR (bottom right), ETR (top right), ETL (top left) and EBL (bottom C left). It also computes the integrals of f'/f round the sides of the C box (right, top, left and bottom sides) taken in anti-clockwise sense C round the boundary. The integrals are in ANSR (right side), ANST (top), C ANSL (left side) and ANSB (bottom). The meshpoints used for these C integrals are (respectively) in C C ZR(I), I = 0,..,NZR; ZT(I), I = 0,..,NZT; C ZL(I), I = 0,..,NZL; ZB(I), I = 0,..,NZB. C C Note that ZR(0) = EBR, ZR(NZR) = ETR = ZT(0), C ZT(NZT) = ETL = ZL(0), ZL(NZL) = EBL = ZB(0), C ZB(NZB) = EBR = ZR(0). C C The function values are stored in C C FR(I,1), I = 0,..,NZR; FT(I,1), I = 0,..,NZT; C FL(I,1), I = 0,..,NZL; FB(I,1), I = 0,..,NZB. C C and the scaling factors are stored in C C FR(I,2), I = 0,..,NZR; FT(I,2), I = 0,..,NZT; C FL(I,2), I = 0,..,NZL; FB(I,2), I = 0,..,NZB. C C (the value of f(z(i)) is F(i,1)*exp(F(i,2))). C C C .. Scalar Arguments .. DOUBLE COMPLEX ANSB,ANSL,ANSR,ANST,ANSWER,EBL,EBR,ETL,ETR INTEGER IFAIL,IFB,IFL,IFR,IFT,IMATCH,IZB,IZL,IZR,IZT,NMESH,NZB, & NZL,NZR,NZT C .. C .. Array Arguments .. DOUBLE COMPLEX FB(0:IFB,1:2),FL(0:IFL,1:2),FR(0:IFR,1:2), & FT(0:IFT,1:2),ZB(0:IZB),ZL(0:IZL),ZR(0:IZR), & ZT(0:IZT) DOUBLE PRECISION XMESH(0:NMESH) C .. C .. Subroutine Arguments .. EXTERNAL SL4BCS,SL4COF C .. C .. Local Scalars .. DOUBLE COMPLEX ELAM,F DOUBLE PRECISION CONST INTEGER IFLAG,IREVCM C .. C .. External Subroutines .. EXTERNAL BXSIDE,SLMISS C .. C .. Intrinsic Functions .. INTRINSIC ATAN C .. C WRITE (6,FMT=*) 'Successfully called MKBOX:' C Do right side: IREVCM = 1 IFLAG = 0 C WRITE (9,FMT=*) 'Doing right side:' 10 CALL BXSIDE(EBR,ETR,ELAM,F,CONST,IREVCM,ZR,IZR,NZR,FR,IFR,ANSR, & IFLAG) IF (IREVCM.NE.0) THEN C Function value required C WRITE (9,FMT=*) 'EBR,ETR,ELAM:',EBR,ETR,ELAM CALL SLMISS(XMESH,NMESH,IMATCH,ELAM,SL4COF,SL4BCS,F,CONST, & IFLAG) C WRITE (6,FMT=*) 'Returned from call to SLMISS' IF (IFLAG.NE.0) THEN IFAIL = IFLAG RETURN END IF C WRITE (9,FMT=*) 'EBR,ETR,ELAM,F:',EBR,ETR,ELAM,F GO TO 10 END IF IF (IFLAG.NE.0) THEN IFAIL = IFLAG RETURN END IF C Done right side. C Do top side: IREVCM = 1 IFLAG = 0 C WRITE (9,FMT=*) 'Doing top:' 20 CALL BXSIDE(ETR,ETL,ELAM,F,CONST,IREVCM,ZT,IZT,NZT,FT,IFT,ANST, & IFLAG) IF (IREVCM.NE.0) THEN C Function value required CALL SLMISS(XMESH,NMESH,IMATCH,ELAM,SL4COF,SL4BCS,F,CONST, & IFLAG) IF (IFLAG.NE.0) THEN IFAIL = IFLAG RETURN END IF C WRITE (9,FMT=*) ELAM,F GO TO 20 END IF IF (IFLAG.NE.0) THEN IFAIL = IFLAG RETURN END IF C Done top side. C Do left side: IREVCM = 1 IFLAG = 0 C WRITE (9,FMT=*) 'Doing left side:' 30 CALL BXSIDE(ETL,EBL,ELAM,F,CONST,IREVCM,ZL,IZL,NZL,FL,IFL,ANSL, & IFLAG) IF (IREVCM.NE.0) THEN C Function value required CALL SLMISS(XMESH,NMESH,IMATCH,ELAM,SL4COF,SL4BCS,F,CONST, & IFLAG) IF (IFLAG.NE.0) THEN IFAIL = IFLAG RETURN END IF C WRITE (9,FMT=*) ELAM,F GO TO 30 END IF IF (IFLAG.NE.0) THEN IFAIL = IFLAG RETURN END IF C Done left side C Do bottom side: IREVCM = 1 IFLAG = 0 C WRITE (9,FMT=*) 'Doing bottom:' 40 CALL BXSIDE(EBL,EBR,ELAM,F,CONST,IREVCM,ZB,IZB,NZB,FB,IFB,ANSB, & IFLAG) IF (IREVCM.NE.0) THEN C Function value required CALL SLMISS(XMESH,NMESH,IMATCH,ELAM,SL4COF,SL4BCS,F,CONST, & IFLAG) IF (IFLAG.NE.0) THEN IFAIL = IFLAG RETURN END IF C WRITE (9,FMT=*) ELAM,F GO TO 40 END IF IF (IFLAG.NE.0) THEN IFAIL = IFLAG RETURN END IF C Have now done all sides ANSWER = (ANSL+ANSR+ANST+ANSB)/ (8.D0*ATAN(1.D0)) IFAIL = 0 RETURN END C ----------------------------------------------------------- SUBROUTINE BXSIDE(ELAM1,ELAM2,ELAM,F,CONST,IREVCM,ZARRAY,NZMAX, & NZACT,FARRAY,NFMAX,ANSWER,IFAIL) C C This routine performs the integral of f'(lambda)/f(lambda) C over the line segment from ELAM1 to ELAM2. Generally this C line segment will be parallel to either the X-axis or the C Y-axis in the complex lambda-plane. Function values are C obtained by reverse communication, for flexibility. C C Documentation: C C ELAM1, ELAM2: The endpoints of the line segment. Must be C set on input. Unchanged on output. C C ELAM: output. On intermediate return with IREVCM.NE.0, ELAM C specifies the point on the line segment at which a value C of f is required. This value of f must be assigned to the C variable F before calling BXSIDE again with all other C parameters UNCHANGED. C C F: input. On first call with IREVCM.EQ.1 F need not be assigned. C On subsequent calls with IREVCM greater than 1, F must be C assigned the value f(ELAM)*exp(-CONST). The value of CONST C (see below) should be chosen by the user to avoid overflow C in the evaluation of F. C C CONST: input. On first call with IREVCM.EQ.1 CONST need not be C assigned. On subsequent calls with IREVCM greater than 1, C CONST should be such that f(ELAM) = F*exp(CONST). C C IREVCM: input/output. On first call, IREVCM must be assigned the C value 1. Thereafter, the user must not change the value of C IREVCM. On an intermediate exit with IREVCM > 1 the user C must assign F and CONST as described above and call BXSIDE C again. On exit with IREVCM = 0, the routine has finished the C integration. C C ZARRAY(0:NZMAX): output. On final exit with IREVCM.EQ.0 this array C stores all the points at which function evaluations were made, C with ZARRAY(0) = ELAM1 and ZARRAY(NZACT) = ELAM2. C C NMAX: input. The dimension of ZARRAY as declared in the calling C (sub)program. C C NZACT: output. On final exit, NZACT specifies the number of mesh C intervals used to perform the integration. C C FARRAY(0:NFMAX,1:2): output. On final exit with IREVCM.EQ.0 the C elements of FARRAY are assigned as follows: C C FARRAY(I,1) = f(ZARRAY(I))*exp(-CONST(ZARRAY(I))) C FARRAY(I,2) = CONST(ZARRAY(I)) C C for I = 0,1,..,NZACT. C C NFMAX: input. The first dimension of FARRAY as declared in the C calling (sub)program. C C ANSWER: output. On final exit with IFAIL.EQ.0, ANSWER specifies C the value of the integral. C C IFAIL: input/output. On entry IFAIL should be assigned the value C 0. On output IFAIL.EQ.0 unless a failure has occurred. The C three nonzero values of IFAIL which may be returned are: C C IFAIL = 15: cannot manage to reduce the stepsize to a small C enough size in a reasonable number of attempts. May C indicate that f is very close to zero, or that there C is an error in some of the user-supplied values of f. C C IFAIL = 16: the stepsize has become so small that reasonable C progress cannot be made. Possible causes: same as for C IFAIL = 15. C C IFAIL = 17: one of NZMAX, NFMAX is too small to contain all C the points and function values required to do the C integral. C C .. Scalar Arguments .. DOUBLE COMPLEX ANSWER,ELAM,ELAM1,ELAM2,F DOUBLE PRECISION CONST INTEGER IFAIL,IREVCM,NFMAX,NZACT,NZMAX C .. C .. Array Arguments .. DOUBLE COMPLEX FARRAY(0:NFMAX,1:2),ZARRAY(0:NZMAX) C .. C .. Local Scalars .. DOUBLE COMPLEX FOLD,HMAX,HSTEP,LGQUOT,QUOT,Z,ZOLD DOUBLE PRECISION ARGQ,ARGQH,COLD,EPS INTEGER I,ISTEP,NQUAD,NMIN LOGICAL GTFAR1,GTFAR2,GTFAR3,GTFAR4 C .. C .. External Functions .. DOUBLE PRECISION X02AJF EXTERNAL X02AJF C .. C .. Intrinsic Functions .. INTRINSIC ABS,ATAN2,CMPLX,DBLE,AIMAG,LOG,MAX,MIN C .. C .. Save statement .. SAVE I,NQUAD,HMAX,HSTEP,ZOLD,Z,FOLD,COLD,QUOT,LGQUOT,ISTEP,EPS C .. C WRITE (6,FMT=*) 'Successfully called BXSIDE' IF (IREVCM.EQ.0) RETURN GO TO (10,20,50) IREVCM 10 ZOLD = ELAM1 ZARRAY(0) = ELAM1 I = 0 ISTEP = 0 NQUAD = 100 NMIN = 10 EPS = 1000.D0*X02AJF(0.D0) HMAX = (ELAM2-ELAM1)/NMIN C WRITE (9,FMT=*) 'ELAM2,ELAM1,HSTEP:',ELAM2,ELAM1,HSTEP ELAM = ZOLD IREVCM = 2 C WRITE (9,FMT=*) 'ZOLD before return',ZOLD RETURN 20 FARRAY(0,1) = F FARRAY(0,2) = CONST FOLD = F COLD = CONST ANSWER = CMPLX(0.D0,0.D0) C Get initial stepsize HSTEP = (ELAM2-ELAM1)/NQUAD 30 IF (ABS(HSTEP).GT.1.D-2) THEN HSTEP = HSTEP/10.D0 C WRITE (9,FMT=*) 'HSTEP (provisional):',HSTEP GO TO 30 END IF 40 Z = ZOLD + HSTEP C WRITE (9,FMT=*) 'ZOLD after return',ZOLD C WRITE (9,FMT=*) 'HSTEP:',HSTEP C WRITE (9,FMT=*) 'ZOLD + HSTEP:',ZOLD + HSTEP C Trap possible overshooting of the interval of integration GTFAR1 = DBLE(ELAM1) .LT. DBLE(ELAM2) .AND. & DBLE(Z) .GT. (DBLE(ELAM2)+EPS) GTFAR2 = DBLE(ELAM1) .GT. DBLE(ELAM2) .AND. & DBLE(Z) .LT. (DBLE(ELAM2)-EPS) GTFAR3 = AIMAG(ELAM1) .LT. AIMAG(ELAM2) .AND. & AIMAG(Z) .GT. (AIMAG(ELAM2)+EPS) GTFAR4 = AIMAG(ELAM1) .GT. AIMAG(ELAM2) .AND. & AIMAG(Z) .LT. (AIMAG(ELAM2)-EPS) C WRITE (7,FMT=*) 'GTFAR1,GTFAR2,GTFAR3,GTFAR4:',GTFAR1,GTFAR2, C & GTFAR3,GTFAR4 IF (GTFAR1 .OR. GTFAR2 .OR. GTFAR3 .OR. GTFAR4) THEN C We have stepped too far: reset Z to the end of the integration C range. C WRITE (9,FMT=*) 'We have stepped too far' C WRITE (9,FMT=*) 'Z is now',ELAM2 Z = ELAM2 END IF ELAM = Z IREVCM = 3 RETURN 50 QUOT = F/FOLD LGQUOT = LOG(QUOT) ARGQ = ATAN2(AIMAG(QUOT),DBLE(QUOT)) ARGQ = ABS(ARGQ) WRITE (9,FMT=*) 'ZOLD:',ZOLD WRITE (9,FMT=*) 'FOLD,COLD:',FOLD,COLD WRITE (9,FMT=*) 'Z:',Z WRITE (9,FMT=*) 'F,CONST:',F,CONST WRITE (9,FMT=*) 'QUOT,ARGQ:',QUOT,ARGQ WRITE (9,FMT=*) ARGQH = ARGQ/ABS(HSTEP) IF (ARGQ.GT.1.0D0) THEN C Step was too large C WRITE (9,FMT=*) 'Step too large' HSTEP = HSTEP/ (1.4D0*ARGQ) ISTEP = ISTEP + 1 IF (ISTEP.LT.20) GO TO 40 IFAIL = 15 IREVCM = 0 RETURN END IF IF (ABS(HSTEP).LT.1.D-14) THEN C The stepsize is too small to make decent progress C IFAIL = 16 IREVCM = 0 RETURN END IF C If we are here then the step was small enough C Can we increase the stepsize? IF (ISTEP.EQ.0 .AND. ARGQ.LT.0.6D0) THEN C We may consider increasing the stepsize for the next step C WRITE (6,FMT=*) 'HSTEP,HMAX:',HSTEP,HMAX HSTEP = HSTEP/MAX(0.2D0,1.2D0*ABS(ARGQ)) IF (ABS(HSTEP).GT.ABS(HMAX)) & HSTEP = HSTEP*ABS(HMAX)/ABS(HSTEP) END IF ANSWER = ANSWER + LGQUOT I = I + 1 ZARRAY(I) = Z FARRAY(I,1) = F FARRAY(I,2) = CMPLX(COLD,0.D0) C WRITE (9,FMT=*) 'ELAM,F,CONST:',Z,F,REAL(COLD) NZACT = I ZOLD = Z FOLD = F COLD = CONST IF (ISTEP.GT.0) ISTEP = ISTEP - 1 IF (ABS(Z-ELAM2).GT.EPS) THEN C We have not yet reached the end of the range of integration IF (I+1.LT.MIN(NZMAX,NFMAX)) GO TO 40 C Out of storage! IFAIL = 17 IREVCM = 0 RETURN END IF C If we are here then we have reached the end of the range C of integration. IREVCM = 0 IFAIL = 0 RETURN END C ------------------------------------------------------------------ SUBROUTINE NZERO(XMESH,NMESH,IMATCH,RELAM,SL4COF,SL4BCS, & SLSMPC,SLSMBC,NZ,ZMIN,FMIN,ZMAX,ARGQHM, & IFAIL) C C This routine computes the number of eigenvalues of our C SLP whose real part is less than RELAM C C .. Parameters .. INTEGER ILOC PARAMETER (ILOC=10) C .. C .. Scalar Arguments .. DOUBLE PRECISION ARGQHM(2),RELAM INTEGER IFAIL,IMATCH,NMESH,NZ DOUBLE COMPLEX ZMIN,FMIN,ZMAX(2) C .. C .. Array Arguments .. DOUBLE PRECISION XMESH(0:NMESH) C .. C .. Subroutine Arguments .. EXTERNAL SL4BCS,SL4COF,SLSMPC,SLSMBC C .. C .. Local Scalars .. DOUBLE COMPLEX CSUM,CSUM2,ELAM,FUN,FUNSMP,ZEVAL, & ZMAX1(2),ZMAX2(2),FUN1 DOUBLE PRECISION ARGQH1(2),ARGQH2(2),CNSMP,CONST, & RHIGH,RLOW,AFMIN INTEGER I,IMLOC,IREVCM,IZC,NEVAL2,NEVALS C .. C .. Local Arrays .. DOUBLE PRECISION XLOC(0:10) C .. C .. External Subroutines .. EXTERNAL SLMISS,ZCOUNT C .. C .. External Functions .. DOUBLE PRECISION X01AAF EXTERNAL X01AAF C .. C .. Intrinsic Functions .. INTRINSIC DBLE,AIMAG,EXP,NINT C .. C Remark: On exit, ARGQHM(1) = most negative value of ARG(Q)/H, C ARGQHM(2) = most positive value of ARG(Q)/H, and ZMAX(1) = value C of lambda where ARGQHM(1) occurs, ZMAX(2) = value of lambda C where ARGQHM(2) occurs. C ZMIN = CMPLX(RELAM,0.D0) AFMIN = 1.D40 FMIN = CMPLX(AFMIN,AFMIN) XLOC(0) = XMESH(0) XLOC(ILOC) = XMESH(NMESH) DO 10 I = 1,ILOC XLOC(I) = (DBLE(ILOC-I)*XLOC(0)+DBLE(I)*XLOC(ILOC))/DBLE(ILOC) 10 CONTINUE IMLOC = NINT(DBLE(ILOC)/2.D0) RLOW = 0.D0 C Set least acceptable approximation to +infinity: RHIGH = 1.D8 IREVCM = 1 IFAIL = 0 20 CALL ZCOUNT(RLOW,RHIGH,RELAM,ZEVAL,FUN,IREVCM,CSUM,IZC,NEVALS, & ZMAX1,ARGQH1,IFAIL) IF (IFAIL.NE.0) THEN RETURN END IF IF (IREVCM.NE.0) THEN ELAM = ZEVAL CALL SLMISS(XMESH,NMESH,IMATCH,ELAM,SL4COF,SL4BCS,FUN,CONST, & IFAIL) IF (IFAIL.NE.0) RETURN C The reference problem has constant coefficients so C we can get away with using a really coarse mesh to C find its miss-distance function. CALL SLMISS(XLOC,ILOC,IMLOC,ELAM,SLSMPC,SLSMBC,FUNSMP,CNSMP, & IFAIL) IF (IFAIL.NE.0) RETURN FUN = (FUN/FUNSMP) C WRITE (7,FMT=*) ABS(FMIN),ABS(FUN),ELAM,ZMIN IF (ABS(FUN).LT.AFMIN) THEN ZMIN = ELAM FMIN = FUN AFMIN = ABS(FMIN) END IF C WRITE (9,FMT=*) 'Z,F(Z):',REAL(ELAM),REAL(AIMAG(ELAM)), C & REAL(FUN),REAL(AIMAG(FUN)) C WRITE (7,FMT=*) GO TO 20 END IF RLOW = 0.D0 C Set greatest acceptable approximation to -infinity: RHIGH = -1.D8 IREVCM = 1 IFAIL = 0 30 CALL ZCOUNT(RLOW,RHIGH,RELAM,ZEVAL,FUN1,IREVCM,CSUM2,IZC,NEVAL2, & ZMAX2,ARGQH2,IFAIL) IF (IFAIL.NE.0) THEN RETURN END IF IF (IREVCM.NE.0) THEN ELAM = ZEVAL C WRITE (9,FMT=*) 'Calling SLMISS (1)' CALL SLMISS(XMESH,NMESH,IMATCH,ELAM,SL4COF,SL4BCS,FUN1,CONST, & IFAIL) IF (IFAIL.NE.0) RETURN C The reference problem has constant coefficients so C we can get away with using a really coarse mesh to C find its miss-distance function. C WRITE (9,FMT=*) 'Call SLMISS (2)' CALL SLMISS(XLOC,ILOC,IMLOC,ELAM,SLSMPC,SLSMBC,FUNSMP,CNSMP, & IFAIL) IF (IFAIL.NE.0) RETURN FUN1 = FUN1/FUNSMP C WRITE (9,FMT=*) 'Z,F(Z):',REAL(ELAM),REAL(AIMAG(ELAM)), C & REAL(FUN1),REAL(AIMAG(FUN1)) C WRITE (7,FMT=*) ABS(FMIN),ABS(FUN1),ELAM,ZMIN IF (ABS(FUN1).LT.AFMIN) THEN C WRITE (7,FMT=*) 'Updating ZMIN' ZMIN = ELAM FMIN = FUN1 AFMIN = ABS(FMIN) END IF GO TO 30 END IF C NZ = NINT(AIMAG(CSUM-CSUM2)) NZ = NINT(AIMAG(CSUM-CSUM2+0.5D0*LOG(FUN1/FUN)/X01AAF(0.D0))) C WRITE (9,FMT=*) 'CSUM-CSUM2 + LOG(FUN1/FUN)/2Pi:', C & CSUM-CSUM2 + 0.5D0*LOG(FUN1/FUN)/X01AAF(0.D0) C WRITE (9,FMT=*) 'No of fn evals:',NEVALS + NEVAL2 C WRITE (9,FMT=*) '0.5D0*LOG(FUN1/FUN)/X01AAF(0.D0):', C & 0.5D0*LOG(FUN1/FUN)/X01AAF(0.D0) IF (ARGQH1(1).LT.ARGQH2(1)) THEN ZMAX(1) = ZMAX1(1) ARGQHM(1) = ARGQH1(1) ELSE ZMAX(1) = ZMAX2(1) ARGQHM(1) = ARGQH2(1) END IF IF (ARGQH1(2).GT.ARGQH2(2)) THEN ZMAX(2) = ZMAX1(2) ARGQHM(2) = ARGQH1(2) ELSE ZMAX(2) = ZMAX2(2) ARGQHM(2) = ARGQH2(2) END IF C WRITE (9,FMT=*) C WRITE (9,FMT=*) 'Re(Lambda):',RELAM C WRITE (9,FMT=*) 'NZ:',NZ IFAIL = 0 RETURN END C ------------------------------------------------------------------ SUBROUTINE ZCOUNT(RLOW,RHIGH,RELAM,ZEVAL,FUN,IREVCM,CSUM,IZC, & NEVALS,ZMAX,ARGQHM,IFAIL) C C Program to do zero count based on integration along C the imaginary axis C C .. Scalar Arguments .. DOUBLE COMPLEX CSUM,FUN,ZEVAL DOUBLE PRECISION RELAM,RHIGH,RLOW INTEGER IFAIL,IREVCM,IZC,NEVALS C .. C .. Local Scalars .. DOUBLE COMPLEX CHIGH,CLOW,F,FOLD,LGQUOT,QUOT,Z,ZOLD DOUBLE PRECISION ARGQ,ARGQH,HSTEP,IMBIT,IMBITN,REBIT INTEGER ISTEP,NQUAD LOGICAL GFENUF C .. C .. Intrinsic Functions .. INTRINSIC ABS,ATAN,ATAN2,CMPLX,DBLE,AIMAG,LOG,MAX,NINT,REAL C .. C .. Save statement .. SAVE NQUAD,REBIT,TRNDIF,CLOW,CHIGH,CTERM,FOLD,F,ZOLD,Z,LGQUOT,ZTR, & ZTRN,ZN,QUOT,IMBIT,IMBITN,HSTEP,ARGQH,ISTEP C .. C .. Array Arguments .. DOUBLE COMPLEX ZMAX(2) DOUBLE PRECISION ARGQHM(2) C .. IF (IREVCM.EQ.0) RETURN C WRITE (7,FMT=*) 'Z,ZOLD',Z,ZOLD C WRITE (7,FMT=*) ' *** REBIT = ***',REBIT GO TO (10,20,40) IREVCM 10 REBIT = RELAM CLOW = CMPLX(REBIT,RLOW) CHIGH = CMPLX(REBIT,RHIGH) ARGQHM(1) = 0.D0 ARGQHM(2) = 0.D0 ZMAX(1) = CLOW ZMAX(2) = CLOW ZOLD = CLOW IMBIT = RLOW C IMBIT = RLOW/ (1.D0+ABS(RLOW)) C IMBIT = 1.D2*RLOW/ (1.D0+1.D2*ABS(RLOW)) IREVCM = 2 ZEVAL = ZOLD RETURN 20 FOLD = FUN C FOLD = FUN(ZOLD) CSUM = CMPLX(0.D0,0.D0) NQUAD = 100 NEVALS = 0 HSTEP = 1.D0/NQUAD IF (RLOW.GT.RHIGH) HSTEP = -HSTEP ISTEP = 0 30 IMBITN = IMBIT + HSTEP C* C Avoid overshooting ends: C* IF ((HSTEP.GT.0.D0.AND.IMBITN.GT.RHIGH) .OR. C* & (HSTEP.LT.0.D0.AND.IMBITN.LT.RHIGH)) THEN C* Z = CHIGH C* ELSE Z = CMPLX(REBIT,IMBITN) C* C Z = CMPLX(REBIT,IMBITN/ (1.D0-ABS(IMBITN))) C* C Z = CMPLX(REBIT,1.D-2*IMBITN/ (1.D0-ABS(IMBITN))) C* C* END IF C* IF ((AIMAG(Z).GT.RHIGH.AND.RLOW.LT.RHIGH) .OR. C* & (AIMAG(Z).LT.RHIGH.AND.RLOW.GT.RHIGH)) THEN C* Z = CMPLX(REBIT,RHIGH) C* END IF ZEVAL = Z IREVCM = 3 RETURN 40 F = FUN C F = FUN(Z) NEVALS = NEVALS + 1 QUOT = F/FOLD LGQUOT = LOG(QUOT) C WRITE (7,FMT=*) 'QUOT=',QUOT C WRITE (7,FMT=*) 'REAL(QUOT):',REAL(QUOT) C WRITE (7,FMT=*) 'IM(QUOT):',AIMAG(QUOT) C WRITE (7,FMT=*) 'ATAN2(0.00001,1.0000):', C & ATAN2(1.D-5,1.D0) ARGQ = ATAN2(AIMAG(QUOT),DBLE(REAL(QUOT))) ARGQH = ARGQ/HSTEP IF (ARGQH.LT.ARGQHM(1)) THEN ZMAX(1) = CMPLX(REBIT,AIMAG(Z+ZOLD)*0.5D0) ARGQHM(1) = ARGQH END IF IF (ARGQH.GT.ARGQHM(2)) THEN ZMAX(2) = CMPLX(REBIT,AIMAG(Z+ZOLD)*0.5D0) ARGQHM(2) = ARGQH END IF C WRITE (7,FMT=*) 'ZMAX(1),ZMAX(2):',ZMAX(1),ZMAX(2) C WRITE (10,FMT=*) 'Z,ARG(Q):',Z,ARGQ C WRITE (7,FMT=*) 'Z,ARG(Q):',Z, C & ARGQ ARGQ = ABS(ARGQ) C WRITE (7,FMT=*) REAL(AIMAG(Z)),REAL(F),REAL(AIMAG(F)) C WRITE (7,FMT=*) 'Z,ZOLD,F,FOLD:',Z,ZOLD,F,FOLD C WRITE (7,FMT=*) 'HSTEP=',HSTEP IF (ARGQ.GT.1.0D0) THEN C WRITE (7,FMT=*) 'Stepsize was too large' C WRITE (7,FMT=*) 'New stepsize is',HSTEP/(1.4D0*ARGQ) C Step was too large HSTEP = HSTEP/ (1.4D0*ARGQ) ISTEP = ISTEP + 1 IF (ISTEP.LT.10) GO TO 30 IFAIL = 5 IREVCM = 0 RETURN END IF IF (ABS(HSTEP).LT.1.D-14) THEN C The stepsize is too small to make decent progress C IFAIL = 6 IREVCM = 0 RETURN END IF C If we are here then the step was small enough C Can we increase the stepsize? C WRITE (7,FMT=*) 'Stepsize is',HSTEP C WRITE (7,FMT=*) 'ISTEP=',ISTEP IF (ISTEP.EQ.0 .AND. ABS(ARGQ).LT.0.6D0) THEN C We may consider increasing the stepsize for the next step C WRITE (7,FMT=*) 'Increasing stepsize' HSTEP = HSTEP/MAX(0.2D0,1.2D0*ABS(ARGQ)) C WRITE (7,FMT=*) 'New stepsize is',HSTEP C IF (ABS(HSTEP).GT.0.1D0) HSTEP = SIGN(1.D0,HSTEP)*0.1D0 END IF CSUM = CSUM + LGQUOT ZOLD = Z FOLD = F IMBIT = IMBITN IF (ISTEP.GT.0) ISTEP = ISTEP - 1 C Check to see if we have gone far enough: GFENUF = (RLOW.LT.RHIGH.AND.AIMAG(Z).GT.RHIGH).OR. & (RLOW.GT.RHIGH.AND.AIMAG(Z).LT.RHIGH) GFENUF = GFENUF.AND.(ABS(ARGQ).LT.0.1D0) IF (.NOT.GFENUF) GO TO 30 C** IF (ABS(RHIGH-AIMAG(Z)).GT.1.D-6) GO TO 30 C We have now finished the quadrature C Divide by 2*Pi CSUM = CSUM/ (8.D0*ATAN(1.D0)) IREVCM = 0 IZC = NINT(AIMAG(CSUM)) RETURN END