c------------------------------------------------------------------- C SUBROUTINE SLEIGN C ********** C SEPTEMBER 1, 2005; P.B. BAILEY, W.N. EVERITT AND A. ZETTL C VERSION 1.3 C ********** SUBROUTINE SLEIGN(A,B,INTAB,P0ATA,QFATA,P0ATB,QFATB,A1,A2,B1,B2, + NUMEIG,EIG,TOL,IFLAG,ISLFUN,SLFUN,NCA,NCB) C ********** C C This subroutine is designed for the calculation of a specified C eigenvalue, EIG, of a Sturm-Liouville problem for the equation C C -(p(x)*y'(x))' + q(x)*y(x) = eig*w(x)*y(x) on (a,b) C C with user-supplied coefficient functions p, q, and w, C and with separated boundary conditions. (For coupled C boundary conditions, see the companion subroutine SLCOUP.) C The problem may be either nonsingular or singular. In C the nonsingular case, boundary conditions are of the form C C A1*y(a) + A2*p(a)*y'(a) = 0 C B1*y(b) + B2*p(b)*y'(b) = 0, C C but are of the form C C A1*[y,u](a) + A2*[y,v](a) = 0 C B1*[y,U](b) + B2*[y,V](a) = 0, C C when the endpoints are singular, of type Limit Circle. C In either case the boundary conditions are prescribed by C specifying the numbers A1, A2, B1, B2. In the singular C case the user must also supply the "boundary condition C functions" u, v near a and/or U, V near b, whichever is C singular, of type Limit Circle. C The index of the desired eigenvalue is specified in NUMEIG C and its requested accuracy in TOL. Initial data for the C associated eigenfunction are also computed along with values C at selected points, if desired, in array SLFUN. C C In addition to the coefficient functions p, q, and w, the user C must supply subroutine UV to describe the boundary condition C when the problem is limit circle. UV can be a dummy subroutine C if the problem is not limit circle. C C The SUBROUTINE statement is C C SUBROUTINE SLEIGN(A,B,INTAB,P0ATA,QFATA,P0ATB,QFATB,A1,A2,B1,B2, C 1 NUMEIG,EIG,TOL,IFLAG,ISLFUN,SLFUN,NCA,NCB) C C where C C A and B are input variables defining the interval. If the C interval is finite, A must be less than B. (See INTAB below.) C C INTAB is an integer input variable specifying the nature of the C interval. It can have four values. C C INTAB = 1 - A and B are finite. C INTAB = 2 - A is finite and B is infinite (+). C INTAB = 3 - A is infinite (-) and B is finite. C INTAB = 4 - A is infinite (-) and B is infinite (+). C C If either A or B is infinite, it is classified singular and C its value is ignored. C C P0ATA, QFATA, P0ATB, and QFATB are input variables set to C 1.0 or -1.0 as the following properties of p, q, and w at C the interval endpoints are true or false, respectively. C C P0ATA - p(a) is zero. (If true, A is singular.) C QFATA - q(a) and w(a) are finite. (If false, A is singular.) C P0ATB - p(b) is zero. (If true, B is singular.) C QFATB - q(b) and w(b) are finite. (If false, B is singular.) C C A1 and A2 are input variables set to prescribe the boundary C condition at A. C C B1 and B2 are input variables set to prescribe the boundary C condition at B. C C NUMEIG is an integer variable. On input, it should be set to C the index of the desired eigenvalue (increasing sequence where C index 0 corresponds to the lowest eigenvalue -- if the C eigenvalues are bounded below -- or to the smallest nonegative C eigenvalue otherwise). On output, it is unchanged unless the C problem (apparently) lacks eigenvalue NUMEIG, in which case it C is reset to the index of the largest eigenvalue that seems to C exist. C C EIG is a variable set on input to 0.0 or to an initial guess of C the eigenvalue. If EIG is set to 0.0, SLEIGN2 will generate C the initial guess. On output, EIG holds the calculated C eigenvalue if IFLAG (see below) signals success. C C TOL is a variable set on input to the desired accuracy of the C eigenvalue. On output, TOL is reset to the accuracy estimated C to have been achieved if IFLAG (see below) signals success. C This accuracy estimate is absolute if EIG is less than one C in magnitude, and relative otherwise. In addition, prefixing C TOL with a negative sign, removed after interrogation, serves C as a flag to request trace output from the calculation. C C IFLAG is an integer output variable set as follows: C C IFLAG = 0 - improper input parameters. C IFLAG = 1 - successful problem solution, within tolerance. C IFLAG = 2 - best problem result, not within tolerance. C IFLAG = 3 - NUMEIG exceeds actual highest eigenvalue index. C IFLAG = 4 - RAY and EIG fail to agree after 5 tries. C IFLAG = 6 - in SECANT-METHOD, ABS(DE) .LT. EPSMIN . C IFLAG = 7 - iterations are stuck in a loop. C IFLAG = 8 - number of iterations has reached the set limit. C IFLAG = 9 - residual truncation error dominates. C IFLAG = 10 - integrator tolerance cannot be reduced. C IFLAG = 11 - no more improvement. C IFLAG = 13 - AA cannot be moved in any further. C IFLAG = 14 - BB cannot be moved in any further. c iflag = 15 - Bad behavior of coefficients at an endpoint. C IFLAG = 16 - Could not get started. C IFLAG = 17 - Failed to get a bracket. C IFLAG = 18 - Estimator failed. C IFLAG = 51 - integration failure after 1st call to INTEG. C IFLAG = 52 - integration failure after 2nd call to INTEG. C IFLAG = 53 - integration failure after 3rd call to INTEG. C IFLAG = 54 - integration failure after 4th call to INTEG. C C ISLFUN is an integer input variable set to the number of C selected eigenfunction values desired. If no values are C desired, set ISLFUN to zero. c (If ISLFUN is set to -1, the result will be that SLEIGN2 c will return to the calling program directly after sampling c the coefficients p,q,w . This device is used only once, c in SUBROUTINE PERIO.) C C SLFUN is an array of length at least 9. On output, the first 9 C locations contain the integration interval and initial data C that completely determine the eigenfunction. C C SLFUN(1) - point where two pieces of eigenfunction Y match. C SLFUN(2) - left endpoint XAA of the (truncated) interval. C SLFUN(3) - value of THETA at XAA. (Y = RHO*sin(THETA)) C SLFUN(4) - value of F at XAA. (RHO = exp(F)) C SLFUN(5) - right endpoint XBB of the (truncated) interval. C SLFUN(6) - value of THETA at XBB. C SLFUN(7) - value of F at XBB. C SLFUN(8) - final value of integration accuracy parameter EPS. C SLFUN(9) - the constant Z in the polar form transformation. C C F(XAA) and F(XBB) are chosen so that the eigenfunction is C continuous in the interval (XAA,XBB) and has weighted (by W) C L2-norm of 1.0 on the interval. If ISLFUN is positive, then C on input the further ISLFUN locations of SLFUN specify the C points, in ascending order, where the eigenfunction values C are desired and on output contain the values themselves. C c nca & ncb are integers which indicate the nature of the c endpoints a & b, respectively. Namely: c nca = 1 : Endpoint a is REGULAR . C 2 : " WEAKLY REGULAR . C 3 : " LIMIT CIRCLE, NON-OSC . C 4 : " LIMIT CIRCLE, OSC . C 5 : " LIMIT POINT, REGULAR, AT FINITE PT. C 6 : " LIMIT POINT, DEFAULT . C 7 : " LIMIT POINT, AT INFINITY OR IRREG. C 8 : " LIMIT POINT, BAD BEHAVIOR AT ENDPT. C ncb = 1 : Endpoint b is REGULAR . C etc. as for nca . C C ********** C INPUT QUANTITIES: A,B,INTAB,P0ATA,QFATA,P0ATB,QFATB,A1,A2,B1,B2, C NUMEIG,EIG,TOL,IFLAG,ISLFUN,SLFUN,NCA,NCB,PR C OUTPUT QUANTITIES: NUMEIG,EIG,TOL,IFLAG,SLFUN, C MFS,MLS,PI,TWOPI,HPI,EPSMIN,Z,JAY,ZEE, C AA,TMID,BB,DTHDAA,DTHDBB,MDTHZ,ADDD C C .. Scalar Arguments .. REAL A,A1,A2,B,B1,B2,EIG,P0ATA,P0ATB,QFATA,QFATB,TOL INTEGER IFLAG,INTAB,ISLFUN,NCA,NCB,NUMEIG C .. C .. Array Arguments .. REAL SLFUN(9) C .. C .. Scalars in Common .. REAL A1S,A2S,AA,ASAV,B1S,B2S,BB,BSAV,DTHDAA,DTHDBB, + EIGSAV,EPSMIN,HPI,P0ATAS,P0ATBS,PI, + QFATAS,QFATBS,TMID,TSAVEL,TSAVER,TWOPI,Z, + LPWA,LPQA,FA,GWA,GQA,LPWB,LPQB,FB,GWB,GQB INTEGER IND,INTSAV,ISAVE,MDTHZ,MFS,MLS,MMWD,T21 LOGICAL ADDD,PR C .. C .. Arrays in Common .. REAL TEE(100),TT(7,2),YS(200),YY(7,3,2),ZEE(100) INTEGER JAY(100),MMW(100),NT(2) C .. C .. Local Scalars .. REAL AA1,AAA,AAF,AAL,AAS,ALFA,ATHETA,BALLPK,BB1,BBB, + BBF,BBL,BBS,BESTAA,BESTBB,BETA,BSTEIG,BSTEPS, + BSTEST,BSTMID,CHLIM,CHNG,CONVC,DE,DEDW,DIST, + DTHDA,DTHDB,DTHDE,DTHDEA,DTHDEB,DTHETA,DTHOLD, + DTHOLY,EEE,EIGLO,EIGLT,EIGPI,EIGRT,EIGUP,EL, + ELIMUP,EMAX,EMIN,EOLD,EPS,EPSL,EPSM,ER1,ER1M, + ESTERR,FLO,FLOUP,FMAX,FUP,GUESS,OLDEST,OLDRAY, + OLRAYS,ONE,PIN,PT2,PT3,RAY,RLX,SAVAA,SAVBB, + SAVERR,SUM,T1,T2,T3,TAU,TAU0,TAUM,TMID1,TMP,U,V, + WL,ZZ INTEGER I,IMAX,IMID,IMIN,J,JFLAG,JJL,JJR,K,LOOP2,LOOP3,MF,ML,NEIG, + NEIGST,NITER,NRAY,NTMP LOGICAL AOK,BOK,BRACKT,BRS,BRSS,CHEPS,CONVRG,ENDA,ENDB,EXIT, + FIRSTT,GESS0,IOSC,LCOA,LCOB,LIMUP,LOGIC,NEWTF,NEWTON, + OLNEWT,ONEDIG,THEGT0,THELT0,TRUNKA,TRUNKB C .. C .. Local Arrays .. REAL DELT(100),DS(100),EIGEST(200,4),PS(100),PSS(100), + QS(100),WS(100),XS(100) C .. C .. External Functions .. REAL EPSLON EXTERNAL EPSLON C .. C .. External Subroutines .. EXTERNAL AABB,EFDATA,EIGFCN,ESTIM,OBTAIN,PRELIM,RESET,SAMPLE,SAVE C .. C .. Intrinsic Functions .. INTRINSIC ABS,ATAN,INT,MAX,MIN,SIGN C .. C .. Common blocks .. COMMON /ALBE/LPWA,LPQA,FA,GWA,GQA,LPWB,LPQB,FB,GWB,GQB COMMON /BCDATA/A1S,A2S,P0ATAS,QFATAS,B1S,B2S,P0ATBS,QFATBS COMMON /DATADT/ASAV,BSAV,INTSAV COMMON /DATAF/EIGSAV,IND COMMON /LP/MFS,MLS COMMON /PASS/YS,MMW,MMWD COMMON /PIE/PI,TWOPI,HPI COMMON /PRIN/PR,T21 COMMON /RNDOFF/EPSMIN COMMON /TDATA/AA,TMID,BB,DTHDAA,DTHDBB,MDTHZ,ADDD COMMON /TEEZ/TEE COMMON /TEMP/TT,YY,NT COMMON /TSAVE/TSAVEL,TSAVER,ISAVE COMMON /Z1/Z COMMON /ZEEZ/JAY,ZEE C .. C To produce printout, set PR = .TRUE. PR = .true. C IFLAG = 1 ONE = 1.0D0 EPSMIN = EPSLON(ONE) if(epsmin .le. 1.0d-12) epsmin = 1.0d-12 PI = 4.0D0*ATAN(ONE) TWOPI = 2.0D0*PI HPI = 0.5D0*PI Z = 1.0D0 NEIG = NUMEIG - 1 C LOGIC = 1 .LE. INTAB .AND. INTAB .LE. 4 .AND. + P0ATA*QFATA*P0ATB*QFATB .NE. 0.0D0 IF (INTAB.EQ.1) LOGIC = LOGIC .AND. A .LT. B IF (.NOT.LOGIC) THEN IFLAG = 0 GO TO 150 END IF C C (JFLAG = 0 INDICATES FAILURE OF THE ESTIMATOR) JFLAG = 1 C CALL SAVE(A,B,INTAB,P0ATA,QFATA,P0ATB,QFATB,A1,A2,B1,B2) C DO 5 I = 1,100 JAY(I) = 0 ZEE(I) = 1.0D0 5 CONTINUE C CALL SAMPLE(NCA,NCB,MF,ML,XS,PS,QS,WS,DS,DELT,EMIN,EMAX,IMIN,IMAX, + AAA,BBB) C (MAKE MF,ML AVALABLE THROUGH COMMON/LP/ ) MFS = MF MLS = ML C EIGPI = NUMEIG*PI PIN = EIGPI + PI TAU = ABS(TOL) TAU0 = 0.001D0 TAUM = MAX(TAU,EPSMIN) LIMUP = .FALSE. ELIMUP = EMAX GUESS = EIG GESS0 = ABS(EIG) .LE. (1D-7) LCOA = NCA .EQ. 4 LCOB = NCB .EQ. 4 IOSC = LCOA .OR. LCOB IF (GESS0) THEN CALL ESTIM(IOSC,PIN,MF,ML,PS,QS,WS,PSS,DS,DELT,TAU0,JJL,JJR, + SUM,LIMUP,ELIMUP,EMAX,IMAX,IMIN,EL,WL,DEDW,BALLPK, + EEE,JFLAG) END IF IF (JFLAG.EQ.0) THEN IF (PR) WRITE (T21,FMT=*) ' ESTIMATOR FAILED ' IFLAG = 18 RETURN END IF CALL PRELIM(MF,ML,EEE,BALLPK,XS,QS,WS,DS,DELT,PS,PSS,TAU0,JJL,JJR, + LIMUP,ELIMUP,EL,WL,DEDW,GUESS,EIG,AAA,AA,BBB,BB,NCA, + NCB,EIGPI,ALFA,BETA,DTHDEA,DTHDEB,IMID,TMID) IF (JAY(3).EQ.0 .AND. ZEE(2).NE.1.0D0) THEN IF (PR) WRITE (T21,FMT=*) ' COULD NOT GET STARTED. ' IFLAG = 16 RETURN END IF C IF (NCA.EQ.8 .OR. NCB.EQ.8) THEN C (WE CAN'T HANDLE THIS KIND OF LIMIT POINT) IFLAG = 15 RETURN END IF C IF (ISLFUN.EQ.-1) THEN SLFUN(1) = TMID SLFUN(2) = AA SLFUN(3) = ALFA SLFUN(5) = BB SLFUN(6) = BETA + EIGPI SLFUN(9) = Z ADDD = .FALSE. MDTHZ = 0 IFLAG = 1 RETURN END IF C C SET LOGICAL VARIABLES: AOK = INTAB .LT. 3.D0 .AND. P0ATA .LT. 0.0D0 .AND. + QFATA .GT. 0.0D0 BOK = (INTAB.EQ.1 .OR. INTAB.EQ.3) .AND. P0ATB .LT. 0.0D0 .AND. + QFATB .GT. 0.0D0 LCOA = NCA .EQ. 4 LCOB = NCB .EQ. 4 TRUNKA = LCOA .OR. (NCA.EQ.6) TRUNKB = LCOB .OR. (NCB.EQ.6) C C END PRELIMINARY WORK, BEGIN MAIN TASK OF COMPUTING EIG. C C LOGICAL VARIABLES HAVE THE FOLLOWING MEANINGS IF TRUE. C AOK - ENDPOINT A IS NOT SINGULAR. C BOK - ENDPOINT B IS NOT SINGULAR. C BRACKT - EIG HAS BEEN BRACKETED. C CONVRG - CONVERGENCE TEST FOR EIG HAS BEEN SUCCESSFULLY PASSED. C NEWTON - NEWTON ITERATION MAY BE EMPLOYED. C THELT0 - LOWER BOUND FOR EIG HAS BEEN FOUND. C THEGT0 - UPPER BOUND FOR EIG HAS BEEN FOUND. C LIMIT - UPPER BOUND EXISTS WITH BOUNDARY CONDITIONS SATISFIED. C ONEDIG - MOST SIGNIFICANT DIGIT CAN BE EXPECTED TO BE CORRECT. C C INITIALIZE SOME OF THE CONTROL VARIABLES EXIT = .FALSE. FIRSTT = .TRUE. LOOP2 = 0 LOOP3 = 0 PT2 = 0.0D0 PT3 = 0.0D0 ENDA = .FALSE. ENDB = .FALSE. EPSM = EPSMIN CHEPS = .FALSE. NEWTF = .FALSE. BRS = .FALSE. BRSS = .FALSE. SAVERR = 1.0D+9 ATHETA = 1.0D+9 BSTEST = 1.0D+9 OLDEST = 1.0D+9 OLDRAY = 1.0D+9 NRAY = 1 AAL = AA BBL = BB AAS = AA BBS = BB AAF = AAA BBF = BBB NEIGST = 0 EIG = EEE BSTEIG = EIG EPS = 0.0001D0 EPSL = EPS BESTAA = AA BESTBB = BB BSTMID = TMID BSTEPS = EPS C 110 CONTINUE C (INITIAL-IZE) BRACKT = .FALSE. CONVRG = .FALSE. THELT0 = .FALSE. THEGT0 = .FALSE. NEWTON = .FALSE. EIGLO = EMIN - 1.0D0 FLO = -5.0D0 FUP = 5.0D0 EIGLT = 0.0D0 EIGRT = 0.0D0 EIGUP = EMAX + 1.0D0 IF (LIMUP) EIGUP = MIN(EMAX,ELIMUP) DTHOLD = 1.0D0 C IF (PR) WRITE (T21,FMT=*) IF (PR) WRITE (T21,FMT=*) + '---------------------------------------------' IF (PR) WRITE (T21,FMT=*) ' INITIAL GUESS FOR EIG = ',EIG IF (PR) WRITE (T21,FMT=*) ' AA,BB = ',AA,BB C (UNTIL(CONVRG .OR. EXIT) DO 120 NITER = 1,40 IF (PR) WRITE (*,FMT=*) IF (PR) WRITE (*,FMT=*) ' ******************** ' IF (PR) WRITE (*,FMT=*) ' EIGENVALUE ',NUMEIG C CALL RESET(AA,BB,ALFA,BETA,EIG,GUESS,MF,ML,QS,WS,IMID,TMID, + LIMUP,ELIMUP,DTHDEA,DTHDEB,THELT0,EIGLO,EIGUP, + BRACKT,NCA,NCB,IFLAG) IF (IFLAG.EQ.11) THEN c EXIT = .TRUE. c GO TO 130 GO TO 333 END IF C AA1 = AA BB1 = BB TMID1 = TMID IFLAG = 1 CALL OBTAIN(NCA,ALFA,DTHDEA,NCB,BETA,DTHDEB,AA1,BB1,EIG,EIGPI, + TMID1,EPS,DTHETA,DTHDE,ONEDIG,DTHDA,DTHDB,ER1, + ER1M,IFLAG) IF (IFLAG.EQ.2) THEN AAS = AA BBS = BB END IF ATHETA = ABS(DTHETA) IF (51.LE.IFLAG .AND. IFLAG.LE.54) THEN EXIT = .TRUE. GO TO 130 ELSE FIRSTT = .FALSE. END IF C----------------------------------------------------------- CHEPS = .FALSE. CONVRG = .FALSE. OLNEWT = NEWTON NEWTON = ABS(DTHETA) .LT. 0.06D0 .AND. BRACKT IF (NEWTON) ONEDIG = ONEDIG .OR. + ABS(DTHETA+ER1) .LT. 0.5D0*DTHOLD IF (.NOT.ONEDIG .AND. BRS) THEN EXIT = .TRUE. IF (PR) WRITE (T21,FMT=*) ' NOT ONEDIG ' GO TO 130 END IF C IF (PR) WRITE (*,FMT=*) ' SET BRACKET ' IF (DTHETA.GT.0.0D0) THEN IF (.NOT.THEGT0 .OR. EIG.LE.EIGUP) THEN THEGT0 = .TRUE. EIGUP = EIG FUP = DTHETA EIGRT = EIG - DTHETA/DTHDE END IF ELSE IF (.NOT.THELT0 .OR. EIG.GE.EIGLO) THEN THELT0 = .TRUE. EIGLO = EIG FLO = DTHETA EIGLT = EIG - DTHETA/DTHDE END IF END IF C C EIG IS BRACKETED WHEN BOTH THEGT0=.TRUE. AND THELT0=.TRUE. C BRACKT = THELT0 .AND. THEGT0 IF (BRACKT) LOOP2 = 0 C (TEST-FOR-CONVERGENCE) C C MEASURE CONVERGENCE AFTER ADDING SEPARATE CONTRIBUTIONS TO ERROR. C FLOUP = MIN(ABS(FLO),ABS(FUP)) T1 = (ABS(DTHETA)+ER1M)/ABS(DTHDE) IF (TRUNKA) THEN T2 = (1.0D0+AA)*ABS(DTHDA)/ABS(DTHDE) PT2 = (AAF-AA)*DTHDA/DTHDE ELSE T2 = 0.0D0 END IF IF (TRUNKB) THEN T3 = (1.0D0-BB)*ABS(DTHDB)/ABS(DTHDE) PT3 = (BBF-BB)*DTHDB/DTHDE ELSE T3 = 0.0D0 END IF IF (PR) WRITE (T21,FMT=*) ' FLO,FUP,FLOUP = ',FLO,FUP,FLOUP IF (PR) WRITE (T21,FMT=*) ' DTHDE,DTHDA,DTHDB = ',DTHDE,DTHDA, + DTHDB ESTERR = T1 + T2 + T3 CONVC = T1 + T2 + T3 IF (BRACKT) THEN TMP = EIGUP - EIGLO c ESTERR = MIN(ESTERR,TMP) CONVC = MIN(ESTERR,TMP) END IF 333 CONTINUE ESTERR = ESTERR/MAX(ONE,ABS(EIG)) CONVC = CONVC/MAX(ONE,ABS(EIG)) NEIGST = NEIGST + 1 EIGEST(NEIGST,1) = EIG EIGEST(NEIGST,2) = DTHETA EIGEST(NEIGST,3) = ESTERR EIGEST(NEIGST,4) = 0.0D0 c CONVRG = ESTERR .LE. TAUM .AND. NEWTON CONVRG = CONVC .LE. TAUM .AND. NEWTON IF (PR) WRITE (T21,FMT=*) ' T1,T2,T3 = ',T1,T2,T3 IF (PR) WRITE (T21,FMT=*) ' PT2,PT3 = ',PT2,PT3 IF (PR) WRITE (T21,FMT=*) ' TMID,EPS = ',TMID,EPS IF (PR) WRITE (T21,FMT=*) ' ONEDIG,BRACKT,NEWTON,CONVRG = ', + ONEDIG,BRACKT,NEWTON,CONVRG IF (PR) WRITE (T21,FMT=*) ' EIG,DTHETA,ESTERR = ',EIG,DTHETA, + ESTERR IF (BRACKT .AND. (ESTERR.LT.BSTEST.OR..NOT.BRS) .AND. + T1.LT.0.1D0) THEN BESTAA = AA BESTBB = BB BSTMID = TMID BSTEPS = EPS BSTEIG = EIG BSTEST = ESTERR BRS = BRACKT IF (BRS) BRSS = BRS IF (PR) WRITE (T21,FMT=*) ' BSTEIG,BSTEST = ',BSTEIG, + BSTEST IF (PR) WRITE (T21,FMT=*) ' BRS = ',BRS END IF IF (THEGT0 .AND. PR) WRITE (T21,FMT=*) ' EIGUP = ', + EIGUP IF (THELT0 .AND. PR) WRITE (T21,FMT=*) ' EIGLO = ', + EIGLO IF (PR) WRITE (T21,FMT=*) + '-------------------------------------------' IF (CONVRG) THEN IF (PR) WRITE (*,FMT=*) ' NUMBER OF ITERATIONS WAS ',NITER IF (PR) WRITE (T21,FMT=*) ' NUMBER OF ITERATIONS WAS ', + NITER IF (PR) WRITE (*,FMT=*) + '-----------------------------------------' GO TO 130 ELSE IF (NEWTON) THEN IF (OLNEWT .AND. ATHETA.GT.0.8D0*ABS(DTHOLD)) THEN IF (PR) WRITE (T21,FMT=*) ' ATHETA,DTHOLD = ', + ATHETA,DTHOLD IF (PR) WRITE (T21,FMT=* + ) ' NEWTON DID NOT IMPROVE EIG ' NEWTF = .TRUE. LOOP3 = LOOP3 + 1 ELSE IF (TRUNKA .OR. TRUNKB) THEN ENDA = ATHETA .LT. 1.0D0 .AND. TRUNKA .AND. + ABS(PT2) .GT. MAX(TAUM,T1) ENDB = ATHETA .LT. 1.0D0 .AND. TRUNKB .AND. + ABS(PT3) .GT. MAX(TAUM,T1) IF (ENDA .OR. ENDB) THEN NEWTON = .FALSE. ELSE IF (((T2+T3).GT.T1) .AND. + (ATHETA.LT.1.0D0) .AND. + (AA.LE.AAF.AND.BB.GE.BBF)) THEN IF (PR) WRITE (*,FMT=* + ) ' RESIDUAL TRUNCATION ERROR DOMINATES ' EXIT = .TRUE. IFLAG = 9 IF (PR) WRITE (T21,FMT=*) ' IFLAG = 9 ' GO TO 130 END IF END IF IF (NEWTF .OR. ENDA .OR. ENDB) THEN IF (PR) WRITE (T21,FMT=*) ' NEWTF,ENDA,ENDB = ', + NEWTF,ENDA,ENDB EXIT = .TRUE. GO TO 130 END IF C (NEWTON'S-METHOD) IF (PR) WRITE (*,FMT=*) ' NEWTON''S METHOD ' RLX = 1.2D0 IF (BRACKT) RLX = 1.0D0 EIG = EIG - RLX*DTHETA/DTHDE IF (EIG.LE.EIGLO .OR. EIG.GE.EIGUP) EIG = 0.5D0* + (EIGLO+EIGUP) IF (PR) WRITE (T21,FMT=*) ' NEWTON: EIG = ',EIG ELSE IF (BRACKT) THEN IF (PR) WRITE (*,FMT=*) ' BRACKET ' C (SECANT-METHOD) IF (PR) WRITE (*,FMT=*) ' DO SECANT METHOD ' FMAX = MAX(-FLO,FUP) EOLD = EIG EIG = 0.5D0* (EIGLO+EIGUP) IF (FMAX.LE.1.5D0) THEN U = -FLO/ (FUP-FLO) DIST = EIGUP - EIGLO EIG = EIGLO + U*DIST V = MIN(EIGLT,EIGRT) IF (EIG.LE.V) EIG = 0.5D0* (EIG+V) V = MAX(EIGLT,EIGRT) IF (EIG.GE.V) EIG = 0.5D0* (EIG+V) DE = EIG - EOLD IF (ABS(DE).LT.EPSMIN) THEN TOL = ABS(DE)/MAX(ONE,ABS(EIG)) IFLAG = 6 EXIT = .TRUE. GO TO 130 END IF END IF IF (PR) WRITE (T21,FMT=*) ' SECANT: EIG = ',EIG ELSE C (TRY-FOR-BRACKET) LOOP2 = LOOP2 + 1 IF (LOOP2.GT.9 .AND. .NOT.LIMUP) THEN IFLAG = 12 IF (PR) WRITE (T21,FMT=*) ' IFLAG = 12 ' EXIT = .TRUE. GO TO 130 END IF IF (EIG.EQ.EEE) THEN IF (GUESS.NE.0.0D0) DEDW = 1.0D0/DTHDE CHNG = -0.6D0* (DEDW+1.0D0/DTHDE)*DTHETA IF (EIG.NE.0.0D0 .AND. ABS(CHNG).GT. + 0.1D0*ABS(EIG)) CHNG = -0.1D0*SIGN(EIG,DTHETA) ELSE CHNG = -1.2D0*DTHETA/DTHDE IF (CHNG.EQ.0.D0) CHNG = 0.1D0*MAX(ONE,ABS(EIG)) IF (PR) WRITE (T21,FMT=* + ) ' IN BRACKET, 1,CHNG = ',CHNG C C LIMIT CHANGE IN EIG TO A FACTOR OF 2. C IF (ABS(CHNG).GT. (1.0D0+2.0D0*ABS(EIG))) THEN TMP = 1.0D0+2.0D0*ABS(EIG) CHNG = SIGN(TMP,CHNG) IF (PR) WRITE (T21,FMT=* + ) ' IN BRACKET, 2,CHNG = ',CHNG ELSE IF (ABS(EIG).GE.1.0D0 .AND. + ABS(CHNG).LT.0.1D0*ABS(EIG)) THEN CHNG = 0.1D0*SIGN(EIG,CHNG) IF (PR) WRITE (T21,FMT=* + ) ' IN BRACKET, 3,CHNG = ',CHNG END IF IF (DTHOLD.LT.0.0D0 .AND. LIMUP .AND. + CHNG.GT. (ELIMUP-EIG)) THEN CHNG = 0.95D0* (ELIMUP-EIG) IF (PR) WRITE (T21,FMT=* + ) ' ELIMUP,EIG,CHNG = ',ELIMUP,EIG,CHNG IF (CHNG.LT.EPSMIN) THEN IF (PR) WRITE (*,FMT=*) ' ELIMUP,EIG = ', + ELIMUP,EIG IF (PR) WRITE (T21,FMT=* + ) ' IN BRACKET, CHNG.LT.EPSMIN ' ENDA = TRUNKA .AND. AA .GT. AAF .AND. + (0.5D0*DTHETA+ (AAF-AA)*DTHDA) .GT. + 0.0D0 ENDB = TRUNKB .AND. BB .LT. BBF .AND. + (0.5D0*DTHETA- (BBF-BB)*DTHDB) .GT. + 0.0D0 IF (.NOT. (ENDA.OR.ENDB)) THEN NUMEIG = NEIG - INT(-DTHETA/PI) IF (PR) WRITE (*,FMT=* + ) ' NEW NUMEIG = ',NUMEIG IF (PR) WRITE (T21, + FMT=*) ' NEW NUMEIG = ',NUMEIG END IF IFLAG = 3 GO TO 150 END IF END IF END IF EOLD = EIG CHLIM = 2.0D0*ESTERR*MAX(ONE,ABS(EIG)) IF (ATHETA.LT.0.06D0 .AND. ABS(CHNG).GT.CHLIM .AND. + CHLIM.NE.0.0D0) CHNG = SIGN(CHLIM,CHNG) IF ((THELT0.AND.CHNG.LT.0.0D0) .OR. + (THEGT0.AND.CHNG.GT.0.0D0)) CHNG = -CHNG EIG = EIG + CHNG IF (PR) WRITE (T21,FMT=*) ' BRACKET: EIG = ',EIG END IF END IF IF (IFLAG.EQ.3) GO TO 130 IF (NITER.GE.3 .AND. DTHOLY.EQ.DTHETA) THEN IFLAG = 7 IF (PR) WRITE (T21,FMT=*) ' IFLAG = 7 ' BSTEIG = EIG BSTEST = ESTERR EXIT = .TRUE. GO TO 130 END IF DTHOLY = DTHOLD DTHOLD = DTHETA IF (PR) WRITE (*,FMT=*) ' NUMBER OF ITERATIONS WAS ',NITER IF (PR) WRITE (*,FMT=*) + '-----------------------------------------------' 120 CONTINUE IFLAG = 8 IF (PR) WRITE (T21,FMT=*) ' IFLAG = 8 ' EXIT = .TRUE. 130 CONTINUE IF (AA.EQ.AAL .AND. BB.EQ.BBL .AND. EPS.LT.EPSL .AND. + ESTERR.GE.0.5D0*SAVERR .AND. .NOT.EXIT) GO TO 140 EPSL = EPS AAL = AA BBL = BB TOL = BSTEST EIG = BSTEIG IF (EXIT) THEN IF (PR) WRITE (T21,FMT=*) ' EXIT ' IF (FIRSTT) THEN IF (IFLAG.EQ.51 .OR. IFLAG.EQ.53) THEN IF (AA.LT.-0.71D0) THEN IF (PR) WRITE (T21,FMT=* + ) ' FIRST COMPLETE INTEGRATION FAILED. ' IF (AA.EQ.-1.0D0) GO TO 150 AAF = AA CALL AABB(AA,-ONE) IF (PR) WRITE (T21,FMT=*) ' AA MOVED FROM ',AAF, + ' IN TO ',AA EXIT = .FALSE. GO TO 110 ELSE IF (PR) WRITE (T21,FMT=*) ' AA.GE.-0.71 ' IFLAG = 13 GO TO 150 END IF ELSE IF (IFLAG.EQ.52 .OR. IFLAG.EQ.54) THEN IF (BB.GT.0.71D0) THEN IF (PR) WRITE (T21,FMT=* + ) ' FIRST COMPLETE INTEGRATION FAILED. ' IF (BB.EQ.1.0D0) GO TO 150 BBF = BB CALL AABB(BB,-ONE) IF (PR) WRITE (T21,FMT=*) ' BB MOVED FROM ',BBF, + ' IN TO ',BB EXIT = .FALSE. GO TO 110 ELSE IF (PR) WRITE (T21,FMT=*) ' BB.LE.0.71 ' IFLAG = 14 GO TO 150 END IF END IF ELSE IF (IFLAG.EQ.51 .OR. IFLAG.EQ.53) THEN IF (PR) WRITE (*,FMT=*) ' A COMPLETE INTEGRATION FAILED. ' IF (PR) WRITE (T21,FMT=*) + ' A COMPLETE INTEGRATION FAILED. ' IF (CHEPS) THEN EPS = 5.0D0*EPS EPSM = EPS IF (PR) WRITE (T21,FMT=*) ' EPS INCREASED TO ',EPS ELSE AAF = AA CALL AABB(AA,-ONE) IF (PR) WRITE (T21,FMT=*) ' AA MOVED FROM ',AAF, + ' IN TO ',AA END IF EXIT = .FALSE. GO TO 110 ELSE IF (IFLAG.EQ.52 .OR. IFLAG.EQ.54) THEN IF (PR) WRITE (*,FMT=*) ' A COMPLETE INTEGRATION FAILED. ' IF (PR) WRITE (T21,FMT=*) + ' A COMPLETE INTEGRATION FAILED. ' IF (CHEPS) THEN EPS = 5.0D0*EPS EPSM = EPS IF (PR) WRITE (T21,FMT=*) ' EPS INCREASED TO ',EPS ELSE BBF = BB CALL AABB(BB,-ONE) IF (PR) WRITE (T21,FMT=*) ' BB MOVED FROM ',BBF, + ' IN TO ',BB END IF EXIT = .FALSE. GO TO 110 ELSE IF (IFLAG.EQ.6) THEN IF (PR) WRITE (*,FMT=*) ' IN SECANT, CHNG.LT.EPSMIN ' IF (PR) WRITE (T21,FMT=*) ' IN SECANT, CHNG.LT.EPSMIN ' GO TO 140 ELSE IF (IFLAG.EQ.7) THEN IF (PR) WRITE (*,FMT=*) ' DTHETA IS REPEATING ' IF (PR) WRITE (T21,FMT=*) ' DTHETA IS REPEATING ' GO TO 140 ELSE IF (IFLAG.EQ.8) THEN IF (PR) WRITE (*,FMT=*) + ' NUMBER OF ITERATIONS REACHED SET LIMIT ' IF (PR) WRITE (T21,FMT=*) + ' NUMBER OF ITERATIONS REACHED SET LIMIT ' GO TO 140 ELSE IF (IFLAG.EQ.9) THEN IF (PR) WRITE (T21,FMT=*) + ' RESIDUAL TRUNCATION ERROR DOMINATES ' GO TO 140 ELSE IF (IFLAG.EQ.11) THEN IF (PR) WRITE (*,FMT=*) + ' IN TRY FOR BRACKET, CHNG.LT.EPSMIN ' IF (PR) WRITE (T21,FMT=*) + ' IN TRY FOR BRACKET, CHNG.LT.EPSMIN ' GO TO 140 ELSE IF (IFLAG.EQ.12) THEN IF (PR) WRITE (*,FMT=*) ' FAILED TO GET A BRACKET. ' IF (PR) WRITE (T21,FMT=*) ' FAILED TO GET A BRACKET. ' GO TO 140 ELSE IF (NEWTF .OR. .NOT.ONEDIG) THEN IF (LOOP3.GE.3) THEN IF (PR) WRITE (T21,FMT=*) + ' NEWTON IS NOT GETTING ANYWHERE ' NEWTF = .FALSE. GO TO 140 END IF IF (PR) WRITE (T21,FMT=*) ' BSTEST,OLDEST = ',BSTEST, + OLDEST IF (EPS.GT.EPSM .AND. BSTEST.LT.OLDEST) THEN CHEPS = .TRUE. SAVERR = ESTERR EPS = 0.2D0*EPS IF (BSTEST.LT.0.001D0) EPS = EPSM IF (PR) WRITE (T21,FMT=*) ' EPS REDUCED TO ',EPS EXIT = .FALSE. NEWTON = .FALSE. OLDEST = BSTEST GO TO 110 ELSE IF (EPS.LE.EPSM) THEN IF (PR) WRITE (T21,FMT=* + ) ' EPS CANNOT BE REDUCED FURTHER. ' IFLAG = 2 GO TO 140 ELSE IF (PR) WRITE (*,FMT=*) ' NO MORE IMPROVEMENT ' IF (PR) WRITE (T21,FMT=*) ' NO MORE IMPROVEMENT ' GO TO 140 END IF END IF ELSE IF (ENDA) THEN CALL AABB(AA,ONE) AA = MAX(AA,AAA) IF (AA.LE.AAF) THEN IF (PR) WRITE (T21,FMT=*) ' NO MORE IMPROVEMENT ' CALL AABB(AA,-ONE) GO TO 140 END IF IF (PR) WRITE (*,FMT=*) ' AA MOVED OUT TO ',AA IF (PR) WRITE (T21,FMT=*) ' AA MOVED OUT TO ',AA EXIT = .FALSE. GO TO 110 ELSE IF (ENDB) THEN CALL AABB(BB,ONE) BB = MIN(BB,BBB) IF (BB.GE.BBF) THEN IF (PR) WRITE (T21,FMT=*) ' NO MORE IMPROVEMENT ' CALL AABB(BB,-ONE) GO TO 140 END IF IF (PR) WRITE (*,FMT=*) ' BB MOVED OUT TO ',BB IF (PR) WRITE (T21,FMT=*) ' BB MOVED OUT TO ',BB EXIT = .FALSE. GO TO 110 END IF END IF 140 CONTINUE C C IF CONVRG IS FALSE, CHECK THAT ANY TRUNCATION ERROR MIGHT POSSIBLY C BE REDUCED OR THAT THE INTEGRATIONS MIGHT BE DONE MORE ACCURATELY. C IF (.NOT.CONVRG .AND. IFLAG.LT.50 .AND. IFLAG.NE.11) THEN SAVAA = AA SAVBB = BB IF (EPS.GT.EPSM .AND. ESTERR.LT.0.5D0*SAVERR) THEN IF (PR) WRITE (T21,FMT=*) ' SAVERR,ESTERR = ',SAVERR, + ESTERR SAVERR = ESTERR EPS = 0.2D0*EPS IF (ESTERR.LT.0.001D0) EPS = EPSM IF (PR) WRITE (T21,FMT=*) ' 2,EPS REDUCED TO ',EPS EXIT = .FALSE. NEWTON = .FALSE. IF (DTHOLY.EQ.DTHETA) THEN BSTEST = ESTERR BSTEIG = EIG END IF OLDEST = BSTEST GO TO 110 ELSE IF (ABS(PT2).GT.TAUM .OR. ABS(PT3).GT.TAUM) THEN IF ((AAS-AAF).GT.2.0D0*EPSMIN .AND. ABS(PT2).GT.TAUM) THEN CALL AABB(AA,ONE) AA = MAX(AA,AAA) IF (AA.GT.AAF .AND. AA.LT.SAVAA) THEN IF (PR) WRITE (*,FMT=*) ' AA MOVED OUT TO ',AA IF (PR) WRITE (T21,FMT=*) ' 3,AA MOVED OUT TO ',AA EXIT = .FALSE. END IF END IF IF ((BBF-BBS).GT.2.0D0*EPSMIN .AND. ABS(PT3).GT.TAUM) THEN CALL AABB(BB,ONE) BB = MIN(BB,BBB) IF (BB.GT.SAVBB .AND. BB.LT.BBF) THEN IF (PR) WRITE (*,FMT=*) ' BB MOVED OUT TO ',BB IF (PR) WRITE (T21,FMT=*) ' 3,BB MOVED OUT TO ',BB EXIT = .FALSE. END IF END IF IF (.NOT.EXIT .AND. (AA.NE.SAVAA.OR.BB.NE.SAVBB)) + GO TO 110 END IF END IF IF (PR) WRITE (T21,FMT=*) 'NUMEIG = ',NUMEIG,' EIG = ',EIG, + ' TOL = ',TOL IF (BRSS .AND. BSTEST.LT.0.05D0) THEN C DO (COMPUTE-EIGENFUNCTION-DATA) EPS = BSTEPS CALL EFDATA(ALFA,DTHDEA,A1,A2,BETA,DTHDEB,B1,B2,EIG,EIGPI,EPS, + AOK,NCA,BOK,NCB,RAY,SLFUN) C C IF NEXT CONDITION IS .TRUE., THEN SOMETHING IS APPARENTLY WRONG C WITH THE ACCURACY OF EIG. BISECT AND GO THROUGH THE LOOP AGAIN. C IF (PR) WRITE (T21,FMT=*) ' EIG,RAY = ',EIG,RAY IF (ABS(RAY-EIG).GT.2.0D0*TAUM*MAX(ONE,ABS(EIG))) THEN NRAY = NRAY + 1 IF (PR) WRITE (*,FMT=*) ' NRAY = ',NRAY IF (PR) WRITE (T21,FMT=*) ' NRAY,RAY,OLDRAY = ',NRAY,RAY, + OLDRAY IF (ESTERR.GE.0.5D0*SAVERR) THEN IFLAG = 2 GO TO 150 END IF TMP = EIG EIG = 0.5D0* (EIG+RAY) OLRAYS = OLDRAY OLDRAY = RAY IF (OLDRAY.NE.OLRAYS .AND. NRAY.LT.2) THEN GO TO 110 ELSE EIG = TMP END IF END IF C DO (GENERATE-EIGENFUNCTION-VALUES) CALL EIGFCN(A1,A2,B1,B2,AOK,BOK,NCA,NCB,SLFUN,ISLFUN) IFLAG = 1 END IF C C IF THE ESTIMATED ACCURACY IMPLIES THAT THE COMPUTED VALUE C OF THE EIGENVALUE IS UNCERTAIN, SIGNAL BY IFLAG = 2. C IF ((ABS(EIG).LE.ONE.AND.TOL.GE.ABS(EIG)) .OR. + (ABS(EIG).GT.ONE.AND.TOL.GE.ONE)) IFLAG = 2 150 CONTINUE IF ((IFLAG.EQ.2) .OR. (6.LE.IFLAG.AND.IFLAG.LE.11)) THEN NTMP = 1 IF (NEIGST.GE.10) NTMP = NEIGST - 9 K = NTMP TMP = EIGEST(NTMP,4) DO 160 I = NTMP,NEIGST IF (EIGEST(I,4).LT.TMP) THEN K = I TMP = EIGEST(I,4) END IF 160 CONTINUE J = K TMP = EIGEST(K,3) DO 155 I = K,NEIGST IF (EIGEST(I,3).LT.TMP) THEN J = I TMP = EIGEST(I,3) END IF 155 CONTINUE EIG = EIGEST(J,1) TOL = EIGEST(J,3) IFLAG = 2 END IF C C TO BE SAFE, RESET AA,BB,EPS:- AA = BESTAA BB = BESTBB EPS = BSTEPS C ZZ = -100000.0D0 IF (NCA.GE.5 .OR. NCB.GE.5) THEN ZZ = 1.0D+19 IF (LIMUP) ZZ = ELIMUP END IF SLFUN(9) = ZZ IF (PR) WRITE (T21,FMT=*) ' BEST AA,BB,TMID,EPS = ',BESTAA,BESTBB, + BSTMID,BSTEPS IF (PR) WRITE (T21,FMT=*) ' BRS,BSTEIG,BSTEST = ',BRS,BSTEIG, + BSTEST IF (PR) WRITE (T21,FMT=*) ' IFLAG = ',IFLAG IF (PR) WRITE (T21,FMT=*) + '********************************************' IF (PR) WRITE (T21,FMT=*) IF (PR) WRITE (T21,FMT=*) ' EIGEST = ' DO 1211 I = 1,NEIGST IF (PR) WRITE (T21,FMT=*) EIGEST(I,1),EIGEST(I,2),EIGEST(I,3) 1211 CONTINUE C----------------------------------------------------------- C RETURN END C SUBROUTINE SAVE(A,B,INTAB,P0ATA,QFATA,P0ATB,QFATB,A1,A2,B1,B2) C THIS PROGRAM SIMPLY SAVES THE CALLING ARGUMENTS IN C TWO COMMON BLOCKS FOR USE WHEREVER NEEDED. C IT IS CALLED BY SLEIGN. C C .. Scalar Arguments .. REAL A,A1,A2,B,B1,B2,P0ATA,P0ATB,QFATA,QFATB INTEGER INTAB C .. C .. Scalars in Common .. REAL A1S,A2S,ASAV,B1S,B2S,BSAV,P0ATAS,P0ATBS,QFATAS, + QFATBS INTEGER INTSAV C .. C .. Common blocks .. COMMON /BCDATA/A1S,A2S,P0ATAS,QFATAS,B1S,B2S,P0ATBS,QFATBS COMMON /DATADT/ASAV,BSAV,INTSAV C .. ASAV = A BSAV = B INTSAV = INTAB P0ATAS = P0ATA QFATAS = QFATA P0ATBS = P0ATB QFATBS = QFATB A1S = A1 A2S = A2 B1S = B1 B2S = B2 IF (A1S.LT.0.0D0) THEN A1S = -A1S A2S = -A2S END IF IF (B1S.LT.0.0D0) THEN B1S = -B1S B2S = -B2S END IF RETURN END C SUBROUTINE THUM(MF,ML,XS) C ********** C C THIS PROGRAM DETERMINES THE NUMBER OF CHANGES OF SIGN OF C THE BOUNDARY CONDITION FUNCTION U (OF THE PAIR (U,V)) C WHICH THE USER SUPPLIES. IT IS NEEDED ONLY IF ONE OF THE C ENDPOINTS OF THE INTERVAL (A,B) IS LCO. C IT IS CALLED BY SLEIGN. C C YS IS LIKE XS, BUT HAS TWICE AS MANY POINTS. C MMW(N) IS THE VALUE OF THE INDEX I OF U(I), MF .LE. I .LE. 2*ML-2, C WHERE U FOR THE NTH TIME CHANGES SIGN FROM - TO + C AND WHERE P*U' IS POSITIVE. C MMWD IS THE NUMBER OF SUCH POINTS OF U. C THE QUANTITIES YS, MMW, MMWD ARE NEEDED BY SUBROUTINE SETTHU. C C INPUT QUANTITIES: MF,ML,XS C OUTPUT QUANTITIES: YS,MMW,MMWD C ********** C .. Scalars in Common .. INTEGER MMWD C .. C .. Arrays in Common .. REAL YS(200) INTEGER MMW(100) C .. C .. Local Scalars .. C REAL PUP,PUP1,TMP1,TMP2,TMP3,TMP4,U,U1 INTEGER I,N C .. C .. Common blocks .. COMMON /PASS/YS,MMW,MMWD C .. C .. Scalar Arguments .. INTEGER MF,ML C .. C .. Array Arguments .. REAL XS(*) C .. C .. External Subroutines .. EXTERNAL UV C .. DO 10 I = 1,99 YS(2*I-1) = XS(I) YS(2*I) = 0.5D0* (XS(I)+XS(I+1)) 10 CONTINUE YS(199) = XS(100) N = 0 U1 = 0.0D0 PUP1 = 0.0D0 DO 20 I = 2*MF - 1,2*ML - 1 CALL UV(YS(I),U,PUP,TMP1,TMP2,TMP3,TMP4) IF (U1.LT.0.0D0 .AND. U.GT.0.0D0 .AND. PUP1.GT.0.0D0) THEN N = N + 1 MMW(N) = I - 1 END IF U1 = U PUP1 = PUP 20 CONTINUE MMWD = N RETURN END C SUBROUTINE AABB(TEND,OUT) C THIS PROGRAM IS FOR THE PURPOSE OF MOVING THE TRUNCATED ENDPOINTS, C AA OR BB, EITHER FURTHER OUT TOWARDS -1 OR +1, OR CLOSER IN, C DEPENDING ON THE SIGN OF OUT. C IT IS CALLED BY SLEIGN AND BY EFDATA. C INPUT QUANTITIES: TEND,OUT C OUTPUT QUANTITIES: TEND C C .. Scalar Arguments .. REAL OUT,TEND C .. C .. Local Scalars .. REAL S,SEND INTEGER I,J C .. C .. Local Arrays .. REAL U(15) C .. C .. Intrinsic Functions .. INTRINSIC ABS C .. U(1) = 0.7D0 U(2) = 0.8D0 U(3) = 0.9D0 U(4) = 0.95D0 U(5) = 0.99D0 U(6) = 0.999D0 U(7) = 0.9999D0 U(8) = 0.99999D0 U(9) = 0.999999D0 U(10) = 0.9999999D0 U(11) = 1.0D0 c S = ABS(TEND) J = 9 DO 10 I = 1,9 IF (U(I).LT.S .AND. S.LE.U(I+1)) J = I 10 CONTINUE IF (OUT.GT.0.0D0) THEN SEND = U(J+1) IF (S.EQ.SEND .AND. J.LT.10) SEND = U(J+2) ELSE SEND = U(J) END IF IF (SEND*TEND.LT.0.0D0) SEND = -SEND TEND = SEND RETURN END C SUBROUTINE OBTAIN(NCA,ALFA,DTHDEA,NCB,BETA,DTHDEB,AA,BB,EIG,EIGPI, + TMID,EPS,DTHETA,DTHDE,ONEDIG,DTHDA,DTHDB,ER1, + ER1M,JFLAG) C C THIS PROGRAM OBTAINS THE DIFFERENCE, DTHETA, BETWEEN THE VALUES C OF THETA OBTAINED BY INTEGRATING THE INITIAL VALUE PROBLEMS FOR C THETA FROM (OR NEAR) THE TWO ENDS OF THE INTERVAL (A,B) TO XMID. C THIS DIFFERENCE VANISHES WHEN THE VALUE BEING USED FOR THE C EIGENPARAMETER, EIG, IS EQUAL TO AN EIGENVALUE FOR THE PROBLEM. C IT IS CALLED BY SLEIGN. C C INPUT QUANTITIES: NCA,ALFA,DTHDEA,NCB,BETA,DTHDEB, C AA,BB,EIG,EIGPI,TMID,EPS,PI,TWOPI,HPI, C A1,A2,P0ATA,QFATA,B1,B2,P0ATB,QFATB, C A,B,INTAB,PR C OUTPUT QUANTITIES: EIGSAV,IND,ISAVE,TSAVEL,TSAVER,DTHDAA, C DTHDBB,DTHETA,DTHDE,ONEDIG,DTHZ, C MDTHZ,TSAVEL,TSAVER,ISAVE,JFLAG, C DTHDA,DTHDB,ER1,ER1M C .. Scalar Arguments .. REAL AA,ALFA,BB,BETA,DTHDA,DTHDB,DTHDE,DTHDEA,DTHDEB, + DTHETA,EIG,EIGPI,EPS,ER1,ER1M,TMID INTEGER JFLAG,NCA,NCB LOGICAL ONEDIG C .. C .. Scalars in Common .. REAL A,A1,A2,AAS,B,B1,B2,BBS,DTHDAA,DTHDBB,EIGSAV,HPI, + P0ATA,P0ATB,PI,QFATA,QFATB,TMIDS,TSAVEL,TSAVER, + TWOPI INTEGER IND,INTAB,ISAVE,MDTHZ,T21 LOGICAL ADDD,PR C .. C .. Local Scalars .. REAL ATHETA,C,DT,DTHZ,ER2,PX,QX,REMZ,THA,THB,WX,X INTEGER IFLAG LOGICAL AOK,BOK,LCA,LCB,LCOA,LCOB,SINGA,SINGB C .. C .. Local Arrays .. REAL ERL(3),ERR(3),YL(3),YR(3),YZL(3),YZR(3) C .. C .. External Functions .. REAL P,Q,W EXTERNAL P,Q,W C .. C .. External Subroutines .. EXTERNAL DXDT,INTEG C .. C .. Intrinsic Functions .. INTRINSIC ABS,COS,EXP,MAX,SIN C .. C .. Common blocks .. COMMON /BCDATA/A1,A2,P0ATA,QFATA,B1,B2,P0ATB,QFATB COMMON /DATADT/A,B,INTAB COMMON /DATAF/EIGSAV,IND COMMON /PIE/PI,TWOPI,HPI COMMON /PRIN/PR,T21 COMMON /TDATA/AAS,TMIDS,BBS,DTHDAA,DTHDBB,MDTHZ,ADDD COMMON /TSAVE/TSAVEL,TSAVER,ISAVE C .. AOK = INTAB .LT. 3.D0 .AND. P0ATA .LT. 0.0D0 .AND. + QFATA .GT. 0.0D0 BOK = (INTAB.EQ.1 .OR. INTAB.EQ.3) .AND. P0ATB .LT. 0.0D0 .AND. + QFATB .GT. 0.0D0 LCA = NCA .EQ. 3 .OR. NCA .EQ. 4 LCB = NCB .EQ. 3 .OR. NCB .EQ. 4 LCOA = NCA .EQ. 4 LCOB = NCB .EQ. 4 SINGA = NCA .GE. 3 SINGB = NCB .GE. 3 c JFLAG = 0 jflag = 0 IND = 1 C INITIALIZE TSAVEL, TSAVER, DTHDAA TSAVEL = -1.0D0 TSAVER = 1.0D0 IF (PR) WRITE (*,FMT=*) ' OBTAIN DTHETA ' THA = ALFA DTHDAA = 0.0D0 IF (SINGA .AND. .NOT.LCA) THEN CALL DXDT(AA,DT,X) PX = P(X) QX = Q(X) WX = W(X) C = EIG*WX - QX DTHDAA = - (COS(ALFA)**2/PX+C*SIN(ALFA)**2)*DT C C TWO SPECIAL CASES FOR DTHDAA . C IF (C.GE.0.0D0 .AND. P0ATA.LT.0.0D0 .AND. + QFATA.LT.0.0D0) DTHDAA = DTHDAA + ALFA*DT/ (X-A) IF (C.GE.0.0D0 .AND. P0ATA.GT.0.0D0 .AND. + QFATA.GT.0.0D0) DTHDAA = DTHDAA + + (ALFA-0.5D0*PI)*DT/ (X-A) END IF THB = BETA DTHDBB = 0.0D0 IF (SINGB .AND. .NOT.LCB) THEN CALL DXDT(BB,DT,X) PX = P(X) QX = Q(X) WX = W(X) C = EIG*WX - QX DTHDBB = - (COS(BETA)**2/PX+C*SIN(BETA)**2)*DT C C TWO SPECIAL CASES FOR DTHDBB . C IF (C.GE.0.0D0 .AND. P0ATB.LT.0.0D0 .AND. + QFATB.LT.0.0D0) DTHDBB = DTHDBB + (PI-BETA)*DT/ (B-X) IF (C.GE.0.0D0 .AND. P0ATB.GT.0.0D0 .AND. + QFATB.GT.0.0D0) DTHDBB = DTHDBB + + (0.5D0*PI-BETA)*DT/ (B-X) END IF C C PASS EIG TO SUBROUTINE F VIA THE COMMON/DATAF/ : EIGSAV = EIG C INTEGRATE FOR YL: C C YL = (THETA,D(THETA)/D(EIG),D(THETA)/DA) C YL(1) = 0.0D0 YL(2) = 0.0D0 YL(3) = 0.0D0 C ISAVE = 0 CALL INTEG(AA,THA,DTHDAA,DTHDEA,TMID,A1,A2,EPS,YL,ERL,AOK,NCA, + IFLAG) IF (IFLAG.EQ.5) THEN JFLAG = 51 IF (PR) WRITE (T21,FMT=*) ' JFLAG = 51 ' GO TO 130 END IF C SET DTHDA: DTHDA = DTHDAA*EXP(-2.0D0*YL(3)) C INTEGRATE FOR YR: C C YR = (THETA,D(THETA)/D(EIG),D(THETA)/DB) C YR(1) = 0.0D0 YR(2) = 0.0D0 YR(3) = 0.0D0 C ISAVE = 0 CALL INTEG(BB,THB,DTHDBB,DTHDEB,TMID,B1,B2,EPS,YR,ERR,BOK,NCB, + IFLAG) IF (IFLAG.EQ.5) THEN JFLAG = 52 IF (PR) WRITE (T21,FMT=*) ' JFLAG = 52 ' GO TO 130 END IF C SET DTHDB: DTHDB = DTHDBB*EXP(-2.0D0*YR(3)) C C SET ER1, ER2, ER1M: ER1 = ERL(1) - ERR(1) ER2 = ERL(2) - ERR(2) ER1M = MAX(ABS(ERL(1)),ABS(ERR(1))) C C IF EITHER ENDPOINT IS LCO, THEN INTEGRATIONS WITH EIG = 0 C MUST ALSO BE CARRIED OUT. C THE VALUE OF THETA WHEN EIG=0 IS USED TO "CALIBRATE" THE C EIGENVALUE INDEXING -- IN THE LCO CASE, ONLY. C IN THIS CASE, SET ISAVE = 1 FOR USE IN SUBROUTINE LCO, C AND INTEGRATE FOR YZL: C IF (LCOA .OR. LCOB) THEN EIGSAV = 0.0D0 ISAVE = 1 CALL INTEG(AA,THA,DTHDAA,DTHDEA,TMID,A1,A2,EPS,YZL,ERL,AOK, + NCA,IFLAG) IF (IFLAG.EQ.5) THEN JFLAG = 53 IF (PR) WRITE (T21,FMT=*) ' JFLAG = 53 ' GO TO 130 END IF ISAVE = 1 CALL INTEG(BB,THB,DTHDBB,DTHDEB,TMID,B1,B2,EPS,YZR,ERR,BOK, + NCB,IFLAG) IF (IFLAG.EQ.5) THEN JFLAG = 54 IF (PR) WRITE (T21,FMT=*) ' JFLAG = 54 ' GO TO 130 END IF EIGSAV = EIG C SET DTHZ, MDTHZ: DTHZ = YZR(1) - YZL(1) MDTHZ = DTHZ/PI REMZ = DTHZ - MDTHZ*PI IF (DTHZ.LT.0.0D0 .AND. REMZ.LT.0.0D0) THEN MDTHZ = MDTHZ - 1 REMZ = REMZ + PI END IF IF (REMZ.GT.3.14D0) MDTHZ = MDTHZ + 1 C RESET ISAVE TO 0: ISAVE = 0 END IF C C DTHETA MEASURES THETA DIFFERENCE FROM LEFT AND RIGHT INTEGRATIONS. C C SET DTHETA, DTHDE: DTHETA = YL(1) - YR(1) - EIGPI IF (LCOA .OR. LCOB) DTHETA = DTHETA + MDTHZ*PI DTHDE = YL(2) - YR(2) IF (PR) WRITE (T21,FMT=*) ' EIG = ',EIG IF (PR) WRITE (T21,FMT=*) ' YL(1),YR(1) = ',YL(1),YR(1) IF (PR) WRITE (T21,FMT=*) ' DTHETA,DTHDE = ',DTHETA,DTHDE IF (PR) WRITE (T21,FMT=*) ' MDTHZ = ',MDTHZ C ATHETA = ABS(DTHETA) ONEDIG = (ABS(ER1).LE.0.5D0*ABS(DTHETA) .AND. + ABS(ER2).LE.0.5D0*ABS(DTHDE)) .OR. + MAX(ATHETA,ABS(ER1)) .LT. 1.0D-6 C RIGHT AFTER LEAVING THIS SUBROUTINE: C WE NEED TO SET AAS = AA, AND BBS = BB , C AND SET ATHETA = ABS(DTHETA). C C WE ALSO NEED TO SET FIRSTT = .FALSE., UNLESS C JFLAG.EQ. 51, 52, 53, OR 54, WHEN WE NEED TO "GOTO 130" C AND SET EXIT = .TRUE. C THIS SUBROUTINE NEEDS TO BE CALLED WITH THE INPUT FLAG C BEING CALLED "IFLAG" -- SO THAT EXTERNALLY THE "IFLAG" C WILL RETURN AS "IFLAG = 52, ETC". 130 CONTINUE IFLAG = JFLAG RETURN END C SUBROUTINE RESET(AA,BB,ALFA,BETA,EIG,GUESS,MF,ML,QS,WS,IMID,TMID, + LIMUP,ELIMUP,DTHDEA,DTHDEB,THELT0,EIGLO,EIGUP, + BRACKT,NCA,NCB,IFLAG) C C THIS PROGRAM IS USED TO RESET THE MATCHING POINT, TMID, AND THE C BOUNDARY VALUES, ALFA & BETA, WHEN NECESSARY. THESE C QUANTITIES CAN DEPEND ON THE CURRENT VALUE BEING USED FOR EIG. C IT IS CALLED BY SLEIGN. C C INPUT QUANTITIES: AA,BB,ALFA,BETA,EIG,GUESS,MF,ML, C QS,WS,IMID,TMID,LIMUP,ELIMUP,DTHDEA,DTHDEB, C THELT0,EIGLO,EIGUP,BRACKT,NCA,NCB, C A,B,INTAB,A1,A2,P0ATA,QFATA,B1,B2,P0ATB,QFATB, C PI,TWOPI,HPI,EPSMIN,PR C OUTPUT QUANTITIES: IFLAG,ALFA,BETA,DTHDEA,DTHDEB,EIG,IMID,TMID, C EIGUP C C .. Scalar Arguments .. REAL AA,ALFA,BB,BETA,DTHDEA,DTHDEB,EIG,EIGLO,EIGUP, + ELIMUP,GUESS,TMID INTEGER IFLAG,IMID,MF,ML,NCA,NCB LOGICAL BRACKT,LIMUP,THELT0 C .. C .. Array Arguments .. REAL QS(100),WS(100) C .. C .. Scalars in Common .. REAL A,A1,A2,B,B1,B2,EPSMIN,HPI,P0ATA,P0ATB,PI,QFATA, + QFATB,TWOPI INTEGER INTAB,T21 LOGICAL PR C .. C .. Local Scalars .. REAL DERIVL,DERIVR,TMP1,TMP2,V INTEGER JFLAG,KFLAG LOGICAL LCA,LCB,SINGA,SINGB C .. C .. External Subroutines .. EXTERNAL ALFBET,SETMID C .. C .. Intrinsic Functions .. INTRINSIC MIN C .. C .. Common blocks .. COMMON /BCDATA/A1,A2,P0ATA,QFATA,B1,B2,P0ATB,QFATB COMMON /DATADT/A,B,INTAB COMMON /PIE/PI,TWOPI,HPI COMMON /PRIN/PR,T21 COMMON /RNDOFF/EPSMIN C .. LCA = NCA .EQ. 3 .OR. NCA .EQ. 4 LCB = NCB .EQ. 3 .OR. NCB .EQ. 4 SINGA = NCA .GE. 3 SINGB = NCB .GE. 3 C (SET-TMID-AND-BOUNDARY-CONDITIONS) IF (PR) WRITE (*,FMT=*) ' SET TMID AND BOUNDARY CONDITIONS ' V = EIG*WS(IMID) - QS(IMID) IF (V.LE.0.0D0) CALL SETMID(MF,ML,EIG,QS,WS,IMID,TMID) C (RESET-BOUNDARY-CONDITIONS) DERIVL = 0.0D0 IF (SINGA) CALL ALFBET(A,INTAB,AA,A1,A2,EIG,P0ATA,QFATA,.TRUE., + ALFA,KFLAG,DERIVL) DTHDEA = DERIVL DERIVR = 0.0D0 IF (SINGB) THEN CALL ALFBET(B,INTAB,BB,B1,B2,EIG,P0ATB,QFATB,.TRUE.,BETA, + JFLAG,DERIVR) BETA = PI - BETA END IF DTHDEB = -DERIVR IF (PR) WRITE (T21,FMT='(A,E22.14,A,E22.14)') ' ALFA=',ALFA, + ' BETA=',BETA C C CHECK THAT BOUNDARY CONDITIONS CAN BE SATISFIED AT SINGULAR C ENDPOINTS. IF NOT, TRY FOR SLIGHTLY ALTERED EIG CONSISTENT C WITH BOUNDARY CONDITIONS. C LIMUP = .TRUE. MEANS THAT THE EIGENVALUES ARE BOUNDED ABOVE, C BY APPROX. ELIMUP. C IF (LIMUP .AND. EIG.NE.GUESS .AND. .NOT.BRACKT) THEN KFLAG = 1 IF (SINGA .AND. .NOT.LCA) CALL ALFBET(A,INTAB,AA,A1,A2,EIG, + P0ATA,QFATA,.TRUE.,TMP1,KFLAG,TMP2) JFLAG = 1 IF (SINGB .AND. .NOT.LCB) CALL ALFBET(B,INTAB,BB,B1,B2,EIG, + P0ATB,QFATB,.TRUE.,TMP1,JFLAG,TMP2) IFLAG = 1 IF ((KFLAG.NE.1.OR.JFLAG.NE.1) .AND. + (THELT0.AND.EIGLO.LT.ELIMUP)) THEN TMP1 = ELIMUP+2.0D0*EPSMIN EIGUP = MIN(TMP1,EIGUP) c EIGUP = MIN(ELIMUP+2.0D0*EPSMIN,EIGUP) IF (EIG.NE.EIGLO .AND. EIG.NE.EIGUP) THEN EIG = 0.05D0*EIGLO + 0.95D0*EIGUP ELSE IFLAG = 11 IF (PR) WRITE (T21,FMT=*) ' IFLAG = 11 ' GO TO 130 END IF END IF END IF C (IMMEDIATELY UPON LEAVING THIS SUBROUTINE WE NEED C TO CHECK FOR IFLAG = 11. IF SO, SET EXIT = .TRUE. C AND GOTO 130.) 130 CONTINUE RETURN END C SUBROUTINE DXDT(T,DT,X) C ********** C C THIS SUBROUTINE TRANSFORMS COORDINATES FROM T ON (-1,1) TO C X ON (A,B) . C C INPUT QUANTITIES: T,A,B INTAB C OUTPUT QUANTITIES: DT,X C ********** C .. Scalars in Common .. REAL A,B INTEGER INTAB C .. C .. Intrinsic Functions .. INTRINSIC ABS C .. C .. Common blocks .. COMMON /DATADT/A,B,INTAB C .. C .. Scalar Arguments .. REAL DT,T,X C .. GO TO (10,20,30,40) INTAB 10 CONTINUE DT = 0.5D0* (B-A) X = 0.5D0* ((B+A)+ (B-A)*T) RETURN 20 CONTINUE DT = 2.0D0/ (1.0D0-T)**2 X = A + (1.0D0+T)/ (1.0D0-T) RETURN 30 CONTINUE DT = 2.0D0/ (1.0D0+T)**2 X = B - (1.0D0-T)/ (1.0D0+T) RETURN 40 CONTINUE DT = 1.0D0/ (1.0D0-ABS(T))**2 X = T/ (1.0D0-ABS(T)) RETURN END C REAL FUNCTION TFROMX(X) C THIS FUNCTION DETERMINES THE VALUE OF T IN (-1,1) WHICH C CORRESPONDS TO ANY VALUE OF X IN (A,B). C INPUT QUANTITIES: X,A,B,INTAB C OUTPUT QUANTITIES: TFROMX C .. Scalar Arguments .. REAL X C .. C .. Scalars in Common .. REAL A,B INTEGER INTAB C .. C .. Intrinsic Functions .. INTRINSIC ABS C .. C .. Common blocks .. COMMON /DATADT/A,B,INTAB C .. GO TO (10,20,30,40) INTAB 10 CONTINUE TFROMX = (2.0D0*X- (B+A))/ (B-A) RETURN 20 CONTINUE TFROMX = (X-A-1.0D0)/ (X-A+1.0D0) RETURN 30 CONTINUE TFROMX = (1.0D0+X-B)/ (1.0D0-X+B) RETURN 40 CONTINUE TFROMX = X/ (1.0D0+ABS(X)) RETURN END C REAL FUNCTION TFROMI(I) C ********** C C THIS FUNCTION DETERMINES THE VALUE OF THE VARIABLE T IN (-1,1) C WHICH CORRESPONDS TO THE SAMPLE POINT INDEX I. C C INPUT QUANTITIES: I C OUTPUT QUANTITIES: TFROMI C ********** C .. Scalar Arguments .. INTEGER I C .. IF (I.LT.8) THEN TFROMI = -1.0D0 + 0.1D0/4.0D0** (8-I) ELSE IF (I.GT.92) THEN TFROMI = 1.0D0 - 0.1D0/4.0D0** (I-92) ELSE TFROMI = 0.0227D0* (I-50) END IF RETURN END C SUBROUTINE EXTRAP(T,TT,EIG,VALUE,DERIV,IFLAG) C ********** C C THIS SUBROUTINE IS CALLED FROM ALFBET IN DETERMINING BOUNDARY C VALUES AT A SINGULAR ENDPOINT OF THE INTERVAL FOR A C STURM-LIOUVILLE PROBLEM IN THE FORM C C -(P(X)*Y'(X))' + Q(X)*Y(X) = EIG*W(X)*Y(X) ON (A,B) C C FOR USER-SUPPLIED COEFFICIENT FUNCTIONS P, Q, AND W. C C EXTRAP, WHICH IN TURN CALLS INTPOL, EXTRAPOLATES THE FUNCTION C C ARCTAN(1.0/SQRT(-P*(EIG*W-Q))) C C FROM ITS VALUES FOR T WITHIN (-1,1) TO AN ENDPOINT. C C INPUT QUANTITIES: T,TT,EIG,PR C OUTPUT QUANTITIES: VALUE,DERIV,IFLAG C SUBPROGRAMS CALLED C C USER-SUPPLIED ..... P,Q,W C C SLEIGN-SUPPLIED .. DXDT,INTPOL C C ********** C .. Local Scalars .. REAL ANS,CTN,ERROR,PROD,PX,QX,T1,TEMP,WX,X INTEGER KGOOD C .. C .. Local Arrays .. REAL FN1(5),XN(5) C .. C .. External Functions .. REAL P,Q,W EXTERNAL P,Q,W C .. C .. External Subroutines .. EXTERNAL DXDT,INTPOL C .. C .. Intrinsic Functions .. INTRINSIC ABS,ATAN,SQRT,TAN C .. C .. Scalar Arguments .. REAL DERIV,EIG,T,TT,VALUE INTEGER IFLAG C .. C .. Scalars in Common .. INTEGER T21 LOGICAL PR C .. C .. Common blocks .. COMMON /PRIN/PR,T21 C .. C (JUST TO MAKE SURE VALUE IS DEFINED, EVEN IF NOT C (NEEDED, SET IT TO 0. VALUE = 0.0D0 C IFLAG = 1 KGOOD = 0 C HERE, COMING FROM ALFBET, TT IS (PROBABLY) AA OR BB T1 = TT XN(1) = TT 10 CONTINUE CALL DXDT(T1,TEMP,X) PX = P(X) QX = Q(X) WX = W(X) PROD = -PX* (EIG*WX-QX) IF (PROD.GT. (1.0D+10)) THEN ANS = 1.D-5 GO TO 20 END IF IF (PROD.LE.0.0D0) THEN T1 = 0.5D0* (T1+T) IF ((1.0D0+ (T1-T)**2).GT.1.0D0) GO TO 10 IF (PROD.GT.-1.0D-6) VALUE = 2.0D0*ATAN(1.0D0) IFLAG = 5 RETURN ELSE KGOOD = KGOOD + 1 XN(KGOOD) = T1 FN1(KGOOD) = ATAN(1.0D0/SQRT(PROD)) T1 = 0.5D0* (T+T1) IF (KGOOD.LT.5) GO TO 10 END IF C AT THIS POINT, THE XN(I) ARE VALUES OF T BETWEEN TT C (AA OR BB) AND 1.0 OR -1.0, OBTAINED BY BISECTION, C AND THE VALUES OF FN1 ARE THE CORRESPONDING VALUES C OF ATAN(1.0/SQRT(PROD)). C IN THIS CALL TO INTPOL, T IS 1.0 OR -1.0 AND C THE T1 SERVES AS ABSERR IN INTPOL. THE RETURNED C VALUE ANS IS THE EXTRAPOLATED VALUE OF C ATAN(1.0/SQRT(PROD)) AT 1.0 OR -1.0 . IF (KGOOD.EQ.5) THEN IFLAG = 1 ELSE IF (PR) WRITE (T21,FMT=*) ' IN EXTRAP, FAILED TO ' IF (PR) WRITE (T21,FMT=*) ' GET 5 VALUES OF XN,FN ' END IF T1 = 0.00001D0 CALL INTPOL(5,XN,FN1,T,T1,3,ANS,ERROR) 20 CONTINUE VALUE = ABS(ANS) CTN = 1.0D0/TAN(VALUE) DERIV = 0.5D0*PX*WX/CTN/ (1.0D0+CTN**2) C RESTORE TT TO ITS ORIGINAL VALUE. TT = XN(1) RETURN END C SUBROUTINE INTPOL(N,XN,FN,X,ABSERR,MAXDEG,ANS,ERROR) C ********** C C THIS SUBROUTINE FORMS AN INTERPOLATING POLYNOMIAL FOR DATA PAIRS. C IT IS CALLED FROM EXTRAP AND FROM EXTR. C IT IS BEING USED TO EXTRAPOLATE THE VALUES FN AT POINTS XN C TO THE VALUE OF F AT X. C THE PARAMETER ABSERR IS THE REQUESTED ACCURACY IN ANS . C C INPUT QUANTITIES: N,XN,FN,X,ABSERR,MAXDEG C OUTPUT QUANTITIES: ANS, ERROR C ********** C .. Local Scalars .. REAL PROD INTEGER I,I1,II,IJ,IK,IKM1,J,K,L,LIMIT C .. C .. Local Arrays .. REAL V(10,10) INTEGER INDEX(10) C .. C .. Intrinsic Functions .. INTRINSIC ABS,MIN C .. C .. Scalar Arguments .. REAL ABSERR,ANS,ERROR,X INTEGER MAXDEG,N C .. C .. Array Arguments .. REAL FN(N),XN(N) C .. L = MIN(MAXDEG,N-2) + 2 LIMIT = MIN(L,N-1) DO 10 I = 1,N V(I,1) = ABS(XN(I)-X) INDEX(I) = I 10 CONTINUE DO 30 I = 1,LIMIT DO 20 J = I + 1,N II = INDEX(I) IJ = INDEX(J) IF (V(II,1).GT.V(IJ,1)) THEN INDEX(I) = IJ INDEX(J) = II END IF 20 CONTINUE 30 CONTINUE PROD = 1.0D0 I1 = INDEX(1) ANS = FN(I1) V(1,1) = FN(I1) DO 50 K = 2,L IK = INDEX(K) V(K,1) = FN(IK) DO 40 I = 1,K - 1 II = INDEX(I) V(K,I+1) = (V(I,I)-V(K,I))/ (XN(II)-XN(IK)) 40 CONTINUE IKM1 = INDEX(K-1) PROD = (X-XN(IKM1))*PROD ERROR = PROD*V(K,K) IF (ABS(ERROR).LE.ABSERR) RETURN ANS = ANS + ERROR 50 CONTINUE ANS = ANS - ERROR RETURN END C REAL FUNCTION EPSLON(X) C C ESTIMATE UNIT ROUNDOFF IN QUANTITIES OF SIZE X. C C THIS PROGRAM SHOULD FUNCTION PROPERLY ON ALL SYSTEMS C SATISFYING THE FOLLOWING TWO ASSUMPTIONS, C 1. THE BASE USED IN REPRESENTING FLOATING POINT C NUMBERS IS NOT A POWER OF THREE. C 2. THE QUANTITY A IN STATEMENT 10 IS REPRESENTED TO C THE ACCURACY USED IN FLOATING POINT VARIABLES C THAT ARE STORED IN MEMORY. C THE STATEMENT NUMBER 10 AND THE GO TO 10 ARE INTENDED TO C FORCE OPTIMIZING COMPILERS TO GENERATE CODE SATISFYING C ASSUMPTION 2. C UNDER THESE ASSUMPTIONS, IT SHOULD BE TRUE THAT, C A IS NOT EXACTLY EQUAL TO FOUR-THIRDS, C B HAS A ZERO FOR ITS LAST BIT OR DIGIT, C C IS NOT EXACTLY EQUAL TO ONE, C EPS MEASURES THE SEPARATION OF 1.0 FROM C THE NEXT LARGER FLOATING POINT NUMBER. C C .. Scalar Arguments .. REAL X C .. C .. Local Scalars .. REAL A,B,C,EPS,FOUR,THREE C .. C .. Intrinsic Functions .. INTRINSIC ABS C .. FOUR = 4.0D0 THREE = 3.0D0 A = FOUR/THREE 10 B = A - 1.0D0 C = B + B + B EPS = ABS(C-1.0D0) IF (EPS.EQ.0.0D0) GO TO 10 EPSLON = EPS*ABS(X) RETURN END C SUBROUTINE ESTPAC(IOSC,MF,ML,EEE,SUM0,QS,WS,DS,DELT,PS,PSS,TAU, + JJL,JJR,SUM,U) C ********** C C THIS SUBROUTINE ESTIMATES THE CHANGE IN 'PHASE ANGLE' IN THE C EIGENVALUE DETERMINATION OF A STURM-LIOUVILLE PROBLEM IN THE FORM C C -(P(X)*Y'(X))' + Q(X)*Y(X) = EIG*W(X)*Y(X) ON (A,B) C C FOR USER-SUPPLIED COEFFICIENT FUNCTIONS P, Q, AND W. C C THE SUBROUTINE APPROXIMATES (BY TRAPEZOIDAL RULE) THE INTEGRAL OF C C SQRT((EIG*W-Q)/P) C C WHERE THE INTEGRAL IS TAKEN OVER THOSE X IN (A,B) FOR WHICH C C (EIG*W-Q)/P .GT. 0 C C ESTPAC IS CALLED BY SUBROUTINES ESTIM AND PRELIM. C C IN THE MAIN IF..THEN..ELSE, THE FIRST PART DOES THE ESTIMATING OF C THE PHASE ANGLE CHANGE, AND THE SECOND PART DETERMINES THE ZEE'S C TO BE USED IN THE PRUFER TRANSFORMATION BY SUBROUTINE GERKZ. C C INPUT QUANTITIES: IOSC,MF,ML,EEE,SUM0,QS,WS,DS,DELT,PS,TAU C OUTPUT QUANTITIES: PSS,JJL,JJR,SUM,U,JAY,ZEE C ********** C .. Local Scalars .. REAL C1,C2,C3,C4,DPSUM,DPSUMT,ONE,PSUM,RT,RTSAV,UT,V, + WSAV,WW,ZAV,ZAVJ INTEGER J,JJ,MF1 C .. C .. Arrays in Common .. REAL ZEE(100) INTEGER JAY(100) C .. C .. Intrinsic Functions .. INTRINSIC ABS,MAX,MIN,MOD,SIGN,SQRT C .. C .. Common blocks .. COMMON /ZEEZ/JAY,ZEE C .. C .. Scalar Arguments .. REAL EEE,SUM,SUM0,TAU,U INTEGER JJL,JJR,MF,ML LOGICAL IOSC C .. C .. Array Arguments .. REAL DELT(100),DS(100),PS(100),PSS(100),QS(100), + WS(100) C .. ONE = 1.0D0 C1 = 0.1D0 C2 = 0.25D0 C3 = 0.5D0 C4 = 0.001D0 C C SUM ACCUMULATES THE INTEGRAL APPROXIMATION. U MEASURES THE TOTAL C LENGTH OF SUBINTERVALS WHERE (EIG*W-Q)/P .GT. 0.0. ZAV IS THE C AVERAGE VALUE OF SQRT((EIG*W-Q)*P) OVER THOSE SUBINTERVALS. C IF (.NOT.IOSC) THEN DO 5 J = 1,100 PSS(J) = 0.0D0 5 CONTINUE JJL = 99 JJR = 1 SUM = 0.0D0 U = 0.0D0 UT = 0.0D0 ZAV = 0.0D0 WSAV = EEE*WS(MF) - QS(MF) IF (WSAV.GT.0.0D0) THEN RTSAV = SIGN(SQRT(WSAV),PS(MF)) ELSE RTSAV = 0.0D0 END IF DO 10 J = MF + 1,ML WW = EEE*WS(J) - QS(J) IF (WW.GT.0.0D0) THEN U = U + DS(J-1) UT = UT + DELT(J-1) RT = SIGN(SQRT(WW),PS(J)) ELSE RT = 0.0D0 END IF IF (WW.EQ.0.0D0 .OR. WSAV.EQ.0.0D0 .OR. + WW.EQ.SIGN(WW,WSAV)) THEN V = RT + RTSAV ELSE V = (WW*RT+WSAV*RTSAV)/ABS(WW-WSAV) END IF WSAV = WW RTSAV = RT PSUM = DS(J-1)*V IF (EEE.EQ.0.0D0) THEN PSS(J) = PSUM ELSE DPSUM = PSUM - PSS(J) DPSUMT = DPSUM*DELT(J-1)/DS(J-1) IF (DPSUMT.GT.C4*TAU) THEN JJL = MIN(JJL,J) JJR = MAX(JJR,J) END IF END IF SUM = SUM + PSUM IF (U.GT.0.0D0) ZAV = ZAV + DELT(J-1)*V*ABS(PS(J)+PS(J-1)) 10 CONTINUE SUM = C3*SUM - SUM0 ZAV = C2*ZAV ELSE JJ = 1 JAY(1) = MF 20 CONTINUE SUM = 0.0D0 U = 0.0D0 UT = 0.0D0 ZAV = 0.0D0 ZAVJ = 0.0D0 MF1 = JAY(JJ) WSAV = EEE*WS(MF1) - QS(MF1) IF (WSAV.GT.0.0D0) THEN RTSAV = SIGN(SQRT(WSAV),PS(MF1)) ELSE RTSAV = 0.0D0 END IF DO 30 J = MF1 + 1,ML WW = EEE*WS(J) - QS(J) IF (WW.GT.0.0D0) THEN U = U + DS(J-1) UT = UT + DELT(J-1) RT = SIGN(SQRT(WW),PS(J)) ELSE RT = 0.0D0 END IF IF (WW.EQ.0.0D0 .OR. WSAV.EQ.0.0D0 .OR. + WW.EQ.SIGN(WW,WSAV)) THEN V = RT + RTSAV ELSE V = (WW*RT+WSAV*RTSAV)/ABS(WW-WSAV) END IF WSAV = WW RTSAV = RT PSUM = DS(J-1)*V SUM = SUM + PSUM IF (U.GT.0.0D0) ZAV = ZAV + DELT(J-1)*V*ABS(PS(J)+PS(J-1)) IF (U.NE.0.0D0) THEN ZAVJ = C2*ZAV/UT IF ((MOD(J-MF1,7).EQ.0) .OR. J.EQ.ML) THEN JJ = JJ + 1 JAY(JJ) = J ZEE(JJ) = ZAVJ IF (ZEE(JJ).NE.0.0D0) ZEE(JJ) = MAX(ZAVJ,C1) IF (ZEE(JJ).EQ.0.0D0) ZEE(JJ) = ONE GO TO 40 END IF END IF 30 CONTINUE 40 CONTINUE IF (J.GT.ML) THEN JJ = JJ + 1 JAY(JJ) = ML ZEE(JJ) = ZAVJ IF (ZEE(JJ).NE.0.0D0) ZEE(JJ) = MAX(ZAVJ,C1) IF (ZEE(JJ).EQ.0.D0) ZEE(JJ) = ONE END IF IF (J.LT.ML) GO TO 20 SUM = C3*SUM ZAV = C2*ZAV END IF RETURN END C SUBROUTINE ALFBET(XEND,INTAB,TT,COEF1,COEF2,EIG,P0,QF,SING,VALUE, + IFLAG,DERIV) C ********** C C THIS SUBROUTINE COMPUTES A BOUNDARY VALUE OF THE PRUEFER ANGLE, C THETA, FOR A SPECIFIED ENDPOINT OF THE INTERVAL FOR A C STURM-LIOUVILLE PROBLEM IN THE FORM C C -(P(X)*Y'(X))' + Q(X)*Y(X) = EIG*W(X)*Y(X) ON (A,B) C C FOR USER-SUPPLIED COEFFICIENT FUNCTIONS P, Q, AND W. IT IS CALLED C FROM RESET AND PRELIM. C BOTH REGULAR AND SINGULAR ENDPOINTS ARE TREATED. C C SUBPROGRAMS CALLED C C USER-SUPPLIED ..... P,Q,W C C SLEIGN-SUPPLIED .. DXDT,EXTRAP C INPUT QUANTITIES: XEND,INTAB,TT,COEF1,COEF2,EIG,P0,QF, C SING,PI,TWOPI,HPI C OUTPUT QUANTITIES: VALUE,DERIV,IFLAG C C ********** C .. Local Scalars .. REAL C,CD,D,HH,ONE,PX,QX,T,TEMP,TTS,WX,X LOGICAL LOGIC C .. C .. External Functions .. REAL P,Q,W EXTERNAL P,Q,W C .. C .. External Subroutines .. EXTERNAL DXDT,EXTRAP C .. C .. Intrinsic Functions .. INTRINSIC ABS,ATAN,SIGN,SQRT C .. C .. Common blocks .. COMMON /PIE/PI,TWOPI,HPI C .. C SET MACHINE DEPENDENT CONSTANT. C C .. Scalar Arguments .. REAL COEF1,COEF2,DERIV,EIG,P0,QF,TT,VALUE,XEND INTEGER IFLAG,INTAB LOGICAL SING C .. C .. Scalars in Common .. REAL HPI,PI,TWOPI C .. ONE = 1.0D0 IFLAG = 1 DERIV = 0.0D0 IF (.NOT.SING) THEN VALUE = 0.5D0*PI IF (COEF1.NE.0.0D0) VALUE = ATAN(-COEF2/COEF1) LOGIC = (TT.LT.0.0D0 .AND. VALUE.LT.0.0D0) .OR. + (TT.GT.0.0D0 .AND. VALUE.LE.0.0D0) IF (LOGIC) VALUE = VALUE + PI ELSE LOGIC = (INTAB.EQ.2 .AND. TT.GT.0.0D0) .OR. + (INTAB.EQ.3 .AND. TT.LT.0.0D0) .OR. INTAB .EQ. 4 .OR. + (P0.GT.0.0D0 .AND. QF.LT.0.0D0) IF (LOGIC) THEN T = SIGN(ONE,TT) TTS = TT C HERE, T IS 1.0 OR -1.0, AND C TTS IS AA OR BB (PROBABLY) C THIS CALL TO EXTRAP EXTRAPOLATES THE VALUES OF C ATAN(1.0/SQRT(-P(EIG*W - Q)) AT POINTS T BETWEEN C AA OR BB AND -1.0 OR 1.0 TO THE ENDPOINT, -1.0 OR 1.0 C IT ALSO RETURNS THE ASSOCIATED VALUE OF DERIV AT THAT C ENDPOINT. DERIV IS D(VALUE)/D(EIG) . CALL EXTRAP(T,TTS,EIG,VALUE,DERIV,IFLAG) ELSE CALL DXDT(TT,TEMP,X) PX = P(X) QX = Q(X) WX = W(X) C = 2.0D0* (EIG*WX-QX) IF (C.LT.0.0D0) THEN VALUE = 0.0D0 IF (P0.GT.0.0D0) VALUE = 0.5D0*PI ELSE HH = ABS(XEND-X) D = 2.0D0*HH/PX CD = C*D*HH IF (P0.GT.0.0D0) THEN VALUE = C*HH IF (CD.LT.1.0D0) VALUE = VALUE/ + (1.0D0+SQRT(1.0D0-CD)) VALUE = VALUE + 0.5D0*PI ELSE VALUE = D IF (CD.LT.1.0D0) VALUE = VALUE/ + (1.0D0+SQRT(1.0D0-CD)) END IF END IF END IF END IF RETURN END C SUBROUTINE SAMPLE(NCA,NCB,MF,ML,XS,PS,QS,WS,DS,DELT,EMIN,EMAX, + IMIN,IMAX,AAA,BBB) C THIS PROGRAM IS FOR THE PURPOSE OF SAMPLING THE COEFFICIENT C FUNCTIONS P, Q, W, PRIMARILY IN ORDER TO BE ABLE TO MAKE A C FIRST ESTIMATE OF THE DESIRED EIGENVALUE. C IT IS CALLED FROM SLEIGN. C THE ARRAY XS CONTAINS THE VALUES OF X IN (A,B) WHICH ARE C USED IN THE SAMPLING PROCESS. THE ARRAY PS CONTAINS THE C CORRESPONDING VALUES OF P. BUT THE ARRAYS QS AND WS CONTAIN C THE CORRESPONDING VALUES NOT OF Q AND W BUT OF C Q/P AND W/P. C ARRAYS DS AND DELT CONTAIN THE CORRESPONDING SAMPLING POINT C INTERVALS OF X AND T. C ALSO COMPUTED ARE MINIMUM AND MAXIMUM VALUES OF Q/W, CALLED C EMIN AND EMAX, WHICH OCCUR AT THE INDEX I EQUAL TO IMIN AND IMAX. C MF AND ML ARE THE FIRST AND LAST VALUES OF THE INDEX I, C 1 .LE. MF .LT. ML .LE. 99. C C INPUT QUANTITIES: NCA,NCB,EPSMIN C OUTPUT QUANTITIES: MF,ML,XS,PS,QS,WS,DS,DELT,EMIN,EMAX, C IMIN,IMAX,AAA,BBB C .. Scalar Arguments .. REAL AAA,BBB,EMAX,EMIN INTEGER IMAX,IMIN,MF,ML,NCA,NCB C .. C .. Array Arguments .. REAL DELT(100),DS(100),PS(100),QS(100),WS(100),XS(100) C .. C .. Scalars in Common .. REAL EPSMIN C .. C .. Local Scalars .. REAL ONE,PX,QX,T,T50,THRESH,TMP,TS,WX,X,X50,XSAV INTEGER I LOGICAL FYNYT,LCOA,LCOB,SINGA,SINGB C .. C .. External Functions .. REAL P,Q,TFROMI,W EXTERNAL P,Q,TFROMI,W C .. C .. External Subroutines .. EXTERNAL DXDT,THUM C .. C .. Intrinsic Functions .. INTRINSIC ABS C .. C .. Common blocks .. COMMON /RNDOFF/EPSMIN C .. SINGA = NCA .GE. 3 SINGB = NCB .GE. 3 LCOA = NCA .EQ. 4 LCOB = NCB .EQ. 4 ONE = 1.0D0 C DO (SAMPLE-COEFFICIENTS) THRESH = 1.D+35 IF (NCA.GE.5) THRESH = 1.D+17 T50 = 0.0D0 CALL DXDT(T50,TMP,X50) XS(50) = X50 TS = T50 PX = P(X50) QX = Q(X50) WX = W(X50) PS(50) = PX QS(50) = QX/PX WS(50) = WX/PX C C EMIN = MIN(Q/W), ACHIEVED AT X FOR INDEX VALUE IMIN. C EMAX = MAX(Q/W), ACHIEVED AT X FOR INDEX VALUE IMAX. C MF AND ML ARE THE LEAST AND GREATEST INDEX VALUES, RESPECTIVELY. C XSAV = X50 EMIN = 0.0D0 IF (QX.NE.0.0D0) EMIN = QX/WX EMAX = EMIN IMIN = 50 IMAX = 50 DO 20 I = 49,1,-1 T = TFROMI(I) CALL DXDT(T,TMP,X) XS(I) = X PX = P(X) QX = Q(X) WX = W(X) PS(I) = PX QS(I) = QX/PX WS(I) = WX/PX DS(I) = XSAV - X DELT(I) = 0.5D0* (TS-T) XSAV = X TS = T C C TRY TO AVOID OVERFLOW BY STOPPING WHEN FUNCTIONS ARE LARGE NEAR A C OR WHEN W IS SMALL NEAR A. C FYNYT = (ABS(WX)+ABS(QX)+1.0D0/ABS(PX)) .LE. THRESH .AND. + WX .GT. EPSMIN IF (QX.NE.0.0D0 .AND. QX/WX.LT.EMIN) THEN EMIN = QX/WX IMIN = I END IF IF (QX.NE.0.0D0 .AND. QX/WX.GT.EMAX) THEN EMAX = QX/WX IMAX = I END IF MF = I IF (.NOT.FYNYT) GO TO 30 20 CONTINUE DO 25 I = 0,-14,-1 T = TFROMI(I) IF (T.LE.-ONE+EPSMIN) GO TO 30 CALL DXDT(T,TMP,X) PX = P(X) QX = Q(X) WX = W(X) FYNYT = (ABS(WX)+ABS(QX)+1.0D0/ABS(PX)) .LE. THRESH .AND. + WX .GT. EPSMIN IF (.NOT.FYNYT) GO TO 30 25 CONTINUE 30 CONTINUE THRESH = 1.D+35 IF (NCB.GE.5) THRESH = 1.D+17 AAA = T IF (.NOT.SINGA) AAA = -ONE TS = T50 XSAV = X50 DO 40 I = 51,99 T = TFROMI(I) CALL DXDT(T,TMP,X) XS(I) = X PX = P(X) QX = Q(X) WX = W(X) PS(I) = PX QS(I) = QX/PX WS(I) = WX/PX DS(I-1) = X - XSAV DELT(I-1) = 0.5D0* (T-TS) XSAV = X TS = T C C TRY TO AVOID OVERFLOW BY STOPPING WHEN FUNCTIONS ARE LARGE NEAR B C OR WHEN W IS SMALL NEAR B. C FYNYT = (ABS(QX)+ABS(WX)+1.0D0/ABS(PX)) .LE. THRESH .AND. + WX .GT. EPSMIN IF (QX.NE.0.0D0 .AND. QX/WX.LT.EMIN) THEN EMIN = QX/WX IMIN = I END IF IF (QX.NE.0.0D0 .AND. QX/WX.GT.EMAX) THEN EMAX = QX/WX IMAX = I END IF ML = I - 1 IF (.NOT.FYNYT) GO TO 50 40 CONTINUE XS(100) = 0.0D0 DO 45 I = 100,114 T = TFROMI(I) IF (T.GE.ONE-EPSMIN) GO TO 50 CALL DXDT(T,TMP,X) PX = P(X) QX = Q(X) WX = W(X) FYNYT = (ABS(WX)+ABS(QX)+1.0D0/ABS(PX)) .LE. THRESH .AND. + WX .GT. EPSMIN IF (.NOT.FYNYT) GO TO 50 45 CONTINUE 50 CONTINUE BBB = T IF (.NOT.SINGB) BBB = ONE C IF (LCOA .OR. LCOB) CALL THUM(MF,ML,XS) C RETURN END C SUBROUTINE ESTIM(IOSC,PIN,MF,ML,PS,QS,WS,PSS,DS,DELT,TAU,JJL,JJR, + SUM,LIMUP,ELIMUP,EMAX,IMAX,IMIN,EL,WL,DEDW, + BALLPK,EEE,JFLAG) C THIS PROGRAM MAKES A FIRST ESTIMATE OF THE DESIRED C EIGENVALUE, USING THE ARRAYS PS,QS,WS OF SAMPLED VALUES OF C THE COEFFICIENT FUNCTIONS P,Q,W OBTAINED BY SUBROUTINE SAMPLE. C IT IS CALLED BY SLEIGN. C THIS ESTIMATE IS RETURNED IN THE VARIABLE EEE. C JJL AND JJR ARE THE MIN AND MAX OF THE INDEX I FOR WHICH C EIG*W - Q IS NONNEGATIVE (OF THE SAMPLED VALUES). C IMIN AND IMAX ARE THE SAMPLE INDICES WHERE THE FUNCTION QS/WS C ATTAINS ITS MINIMUM AND MAXIMUM OF EMIN AND EMAX. C WHEN THE ESTIMATING PROCESS BREAKS DOWN, A "BALLPARK" ESTIMATE, C BALLPK, IS DETERMINED. C C BASICALLY, THE ESTIMATE, EEE, IS THE SOLUTION OF THE C EQUATION C |b C INTEGRAL|SQRT((EEE*W-Q)/P)*DX = (NUMEIG+1)*PI C |a C C WHERE THE INTEGRATION IS OVER THOSE X FOR WHICH THE INTEGRAND C IS REAL. C INPUT QUANTITIES: IOSC,PIN,MF,ML,PS,QS,WS,DS,DELT,TAU, C JJL,JJR,LIMUP,ELIMUP,EMIN,EMAX,IMIN,IMAX C OUTPUT QUANTITIES: PSS,JJL,JJR,SUM,IMIN,IMAX,EL,WL,DEDW, C BALLPK,EEE,JFLAG C C .. Scalar Arguments .. REAL BALLPK,DEDW,EEE,EL,ELIMUP,EMAX,PIN,SUM,TAU,WL INTEGER IMAX,IMIN,JFLAG,JJL,JJR,MF,ML LOGICAL IOSC,LIMUP C .. C .. Array Arguments .. REAL DELT(100),DS(100),PS(100),PSS(100),QS(100), + WS(100) C .. C .. Scalars in Common .. REAL A,B INTEGER INTAB,T21 LOGICAL PR C .. C .. Local Scalars .. REAL EU,FNEW,FOLD,ONE,SUM0,U,ULO,UUP,WU INTEGER IE,JLOOP LOGICAL LOGIC C .. C .. External Subroutines .. EXTERNAL ESTPAC C .. C .. Intrinsic Functions .. INTRINSIC ABS,MAX,MIN C .. C .. Common blocks .. COMMON /DATADT/A,B,INTAB COMMON /PRIN/PR,T21 C .. JFLAG = 1 ONE = 1.0D0 SUM0 = 0.0D0 C ------------------------ C WHEN ESTPAC IS CALLED, THE RETURNED VALUE OF SUM (OR SUM0) IS C THE APPROXIMATION TO THE ABOVE INTEGRAL, C THE VALUE OF U IS THE TOTAL LENGTH OF THOSE SUBINTERVALS FOR C WHICH THE INTEGRAND IS REAL. IF (IOSC) THEN EEE = 0.0D0 CALL ESTPAC(.FALSE.,MF,ML,EEE,SUM0,QS,WS,DS,DELT,PS,PSS,TAU, + JJL,JJR,SUM,U) SUM0 = SUM END IF C ------------------------ EEE = MIN(ELIMUP,EMAX) CALL ESTPAC(.FALSE.,MF,ML,EEE,SUM0,QS,WS,DS,DELT,PS,PSS,TAU,JJL, + JJR,SUM,U) C ------------------------ EU = EEE WU = SUM UUP = U C -------------------------- C PIN IS NORMALLY SET = (NUMEIG+1)*PI, EXCEPT WHEN IOSC = .TRUE. 55 CONTINUE IF (.NOT.LIMUP .AND. ABS(SUM).GE.10.0D0*MAX(ONE,ABS(PIN))) THEN IF (SUM.GE.10.0D0*PIN) THEN IF (EEE.GE.ONE) THEN EEE = EEE/10.0D0 ELSE IF (EEE.LT.-ONE) THEN EEE = 10.0D0*EEE ELSE EEE = EEE - ONE END IF ELSE IF (EEE.LE.-ONE) THEN EEE = EEE/10.0D0 ELSE IF (EEE.GT.ONE) THEN EEE = 10.0D0*EEE ELSE EEE = EEE + ONE END IF END IF CALL ESTPAC(.FALSE.,MF,ML,EEE,SUM0,QS,WS,DS,DELT,PS,PSS,TAU, + JJL,JJR,SUM,U) GO TO 55 END IF C ------------------------------- C THE LOCAL VARIABLES EL, WL ARE USED TO INDICATE VALUE OF EIG C AND VALUE OF SUM WHEN SUM .LT. PIN; C EU, WU ARE USED TO INDICATE VALUE OF EIG AND VALUE OF SUM C WHEN SUM .GT. PIN. C ULO, UUP ARE CORRESPONDING VALUES FOR U, THE TOTAL LENGTH C OF INTERVAL FOR WHICH THE INTEGRAND IS POSITIVE. EU = EEE WU = SUM UUP = U IF (SUM.GE.PIN) THEN EL = EU WL = WU ULO = UUP C ---------------------- JLOOP = 0 60 CONTINUE IF (WL.GE.PIN) THEN C REDUCE EEE: EU = EL WU = WL UUP = ULO EEE = EL - ((WL-PIN+3.0D0)/U)**2 - ONE CALL ESTPAC(.FALSE.,MF,ML,EEE,SUM0,QS,WS,DS,DELT,PS,PSS, + TAU,JJL,JJR,SUM,U) EL = EEE WL = SUM ULO = U IF (SUM.EQ.0.0D0) JLOOP = JLOOP + 1 IF (JLOOP.GE.5) THEN C (ESTIMATOR FAILURE) JFLAG = 0 RETURN END IF GO TO 60 END IF C ----------------------- ELSE C INCREASE EEE: EL = EEE WL = SUM ULO = U END IF IF (LIMUP .AND. WU.LT.PIN) THEN EEE = ELIMUP ELSE IF (UUP.EQ.0.0D0) THEN EEE = EMAX + ONE CALL ESTPAC(.FALSE.,MF,ML,EEE,SUM0,QS,WS,DS,DELT,PS,PSS, + TAU,JJL,JJR,SUM,U) EU = EEE WU = SUM UUP = U END IF 70 CONTINUE IF (WU.LE.PIN) THEN C INCREASE EEE: EL = EU WL = WU ULO = UUP EEE = EU + ((PIN-WU+3.0D0)/U)**2 + ONE CALL ESTPAC(.FALSE.,MF,ML,EEE,SUM0,QS,WS,DS,DELT,PS,PSS, + TAU,JJL,JJR,SUM,U) EU = EEE WU = SUM UUP = U GO TO 70 END IF C -------------------------------- 80 CONTINUE C DETERMINE THE INDICES IMIN, IMAX, WHERE THE FUNCTION C QS/WS ATTAINS ITS MINIMUM AND MAXIMUM VALUES OF EMIN, EMAX: IF (ABS(IMAX-IMIN).GE.2 .AND. EU.LE.EMAX) THEN IE = (IMAX+IMIN)/2 EEE = QS(IE)/WS(IE) CALL ESTPAC(.FALSE.,MF,ML,EEE,SUM0,QS,WS,DS,DELT,PS,PSS, + TAU,JJL,JJR,SUM,U) IF (SUM.GT.PIN) THEN IMAX = IE WU = SUM EU = EEE UUP = U ELSE IMIN = IE WL = SUM EL = EEE ULO = U END IF GO TO 80 END IF C --------------------------------- C C IMPROVE APPROXIMATION FOR EIG USING BISECTION OR SECANT METHOD. C SUBSTITUTE 'BALLPARK' ESTIMATE IF APPROXIMATION GROWS TOO LARGE. C DEDW = (EU-EL)/ (WU-WL) FOLD = 0.0D0 IF (INTAB.EQ.1) BALLPK = (PIN/ (A-B))**2 IF (INTAB.EQ.1 .AND. PR) WRITE (T21,FMT=*) ' BALLPK = ',BALLPK C --------------------------------- LOGIC = .TRUE. 90 CONTINUE C NOW TRY TO REFINE THE ESTIMATE EEE: C USE A SECANT METHOD. IF (LOGIC) THEN LOGIC = (WL.LT.PIN-ONE .OR. WU.GT.PIN+ONE) EEE = EL + DEDW* (PIN-WL) FNEW = MIN(PIN-WL,WU-PIN) IF (FNEW.GT.0.4D0*FOLD .OR. FNEW.LE.ONE) EEE = 0.5D0* + (EL+EU) IF (INTAB.EQ.1 .AND. ABS(EEE).GT.1.0D3*BALLPK) THEN EEE = BALLPK GO TO 100 ELSE IF (INTAB.NE.1 .AND. ABS(EEE).GT.1.0D6) THEN EEE = ONE GO TO 100 ELSE FOLD = FNEW CALL ESTPAC(.FALSE.,MF,ML,EEE,SUM0,QS,WS,DS,DELT,PS, + PSS,TAU,JJL,JJR,SUM,U) IF (SUM.LT.PIN) THEN EL = EEE WL = SUM ULO = U ELSE EU = EEE WU = SUM UUP = U END IF DEDW = (EU-EL)/ (WU-WL) GO TO 90 END IF END IF C ----------------------------- END IF 100 CONTINUE RETURN END C SUBROUTINE PRELIM(MF,ML,EEE,BALLPK,XS,QS,WS,DS,DELT,PS,PSS,TAU, + JJL,JJR,LIMUP,ELIMUP,EL,WL,DEDW,GUESS,EIG,AAA, + AA,BBB,BB,NCA,NCB,EIGPI,ALFA,BETA,DTHDEA,DTHDEB, + IMID,TMID) C THIS PROGRAM, USING THE FIRST ESTIMATE FOR THE EIGENVALUE C OBTAINED BY ESTIM, C SETS THE INITIAL INTERVAL (POSSIBLY A SUBINTERVAL OF (A,B)) C TO BE USED FOR THE PROBLEM. IT IS CALLED BY SLEIGN. C IT DETERMINES ALFA, BETA (THE C BOUNDARY VALUES OF THE PRUEFER ANGLE THETA), TRIES TO IMPROVE C THE ESTIMATE OF THE EIGENVALUE (TAKING ALFA, BETA INTO ACCOUNT) C AND SETS THE MIDPOINT, TMID, IN (-1,1) TO BE USED IN THE C COMPUTATIONS. C C IT ALSO CHECKS TO SEE IF AN ENDPOINT IS LP, AND IF SO, WHETHER C THE POINT SPECTRUM MIGHT BE BOUNDED ABOVE. C INPUT QUANTITIES: MF,ML,EEE,BALLPK,XS,QS,WS,DS,DELT,PS,PSS, C TAU,JJL,JJR,LIMUP,ELIMUP,EL,WL,DEDW,GUESS, C EIG,NCA,NCB,EIGPI C OUTPUT QUANTITIES: EEE,PSS,JJL,JJR,TEE,JAY,ZEE,LIMUP,ELIMUP, C AAA,AA,BBB,BB,ALFA,BETA,DTHDEA,DTHDEB, C IMID,TMID C .. Scalar Arguments .. REAL AA,AAA,ALFA,BALLPK,BB,BBB,BETA,DEDW,DTHDEA, + DTHDEB,EEE,EIG,EIGPI,EL,ELIMUP,GUESS,TAU,TMID,WL INTEGER IMID,JJL,JJR,MF,ML,NCA,NCB LOGICAL LIMUP C .. C .. Array Arguments .. REAL DELT(100),DS(100),PS(100),PSS(100),QS(100), + WS(100),XS(100) C .. C .. Scalars in Common .. REAL A,A1,A2,B,B1,B2,HPI,P0ATA,P0ATB,PI,QFATA,QFATB, + TWOPI INTEGER INTAB,MFS,MLS,T21 LOGICAL PR C .. C .. Arrays in Common .. REAL TEE(100),ZEE(100) INTEGER JAY(100) C .. C .. Local Scalars .. REAL BOUNDA,BOUNDB,DERIVL,DERIVR,ELIMA,ELIMB,ONE,PIN, + SUM,SUM0,THOUS,U,TMP INTEGER I,IPA,IPB,JFLAG,KFLAG LOGICAL GESS0,LCOA,LCOB,LIMA,LIMB,SINGA,SINGB C .. C .. External Functions .. REAL TFROMI EXTERNAL TFROMI C .. C .. External Subroutines .. EXTERNAL ALFBET,ENDPT,ESTPAC,SETMID C .. C .. Intrinsic Functions .. INTRINSIC ABS,MAX,MIN,SIGN C .. C .. Common blocks .. COMMON /BCDATA/A1,A2,P0ATA,QFATA,B1,B2,P0ATB,QFATB COMMON /DATADT/A,B,INTAB COMMON /LP/MFS,MLS COMMON /PIE/PI,TWOPI,HPI COMMON /PRIN/PR,T21 COMMON /TEEZ/TEE COMMON /ZEEZ/JAY,ZEE C .. THOUS = 1000.0D0 PIN = EIGPI + PI SINGA = NCA .GE. 3 SINGB = NCB .GE. 3 LCOA = NCA .EQ. 4 LCOB = NCB .EQ. 4 GESS0 = ABS(GUESS) .LE. (1D-7) ONE = 1.0D0 SUM0 = 0.0D0 C LIMUP = .TRUE. MEANS IT APPEARS THE EIGENVALUES ARE BOUNDED ABOVE. C IN THIS CASE, ELIMUP IS THE ESTIMATED UPPER BOUND -- OR, RATHER, C THE LOWER BOUND OF THE CONTINUOUS SPECTRUM. IF (LIMUP .AND. EEE.GE.ELIMUP) EEE = ELIMUP - 1.0D0 C SET-INITIAL-INTERVAL-AND-MATCHPOINT IF (.NOT.GESS0) THEN EEE = EIG CALL ESTPAC(.FALSE.,MF,ML,EEE,SUM0,QS,WS,DS,DELT,PS,PSS,TAU, + JJL,JJR,SUM,U) END IF C C CHOOSE INITIAL INTERVAL AS LARGE AS POSSIBLE THAT AVOIDS OVERFLOW. C JJL AND JJR ARE BOUNDARY INDICES FOR NONNEGATIVITY OF EIG*W-Q. C IF (.NOT.SINGA) THEN AA = -ONE ELSE IF (LCOA) THEN AA = -0.9999D0 AAA = AA ELSE AA = TFROMI(JJL) AAA = MIN(AAA,AA) AA = MAX(AA,AAA) TMP = -0.99D0 AA = MIN(AA,TMP) AA = MAX(AA,AAA) END IF IF (.NOT.SINGB) THEN BB = ONE ELSE IF (LCOB) THEN BB = 0.9999D0 BBB = BB ELSE BB = TFROMI(JJR) BBB = MAX(BBB,BB) BB = MIN(BB,BBB) TMP = 0.99D0 BB = MAX(BB,TMP) BB = MIN(BB,BBB) END IF C C DETERMINE BOUNDARY VALUES ALFA AND BETA FOR THETA AT A AND B. C CALL ALFBET(A,INTAB,AA,A1,A2,EEE,P0ATA,QFATA,SINGA,ALFA,KFLAG, + DERIVL) CALL ALFBET(B,INTAB,BB,B1,B2,EEE,P0ATB,QFATB,SINGB,BETA,JFLAG, + DERIVR) IF (SINGB) BETA = PI - BETA C C TAKE BOUNDARY CONDITIONS INTO ACCOUNT IN ESTIMATION OF EIG. C PIN = EIGPI + BETA - ALFA IF (LCOA) PIN = PIN + ALFA IF (LCOB) PIN = PIN + PI - BETA IF (GESS0) THEN EEE = EL + DEDW* (PIN-WL) IF (.NOT. (LCOA.OR.LCOB) .AND. + ABS(EEE).GT.THOUS) EEE = SIGN(THOUS,EEE) IF (INTAB.EQ.1 .AND. ABS(EEE).GT.1.0D3*BALLPK) EEE = BALLPK END IF CALL ESTPAC(.TRUE.,MF,ML,EEE,SUM0,QS,WS,DS,DELT,PS,PSS,TAU,JJL, + JJR,SUM,U) C C RESET BOUNDARY VALUES ALFA AND BETA, WHICH MAY DEPEND UPON C THE CURRENT VALUE FOR THE EIGENPARAMETER. C CALL ALFBET(A,INTAB,AA,A1,A2,EEE,P0ATA,QFATA,SINGA,ALFA,KFLAG, + DERIVL) CALL ALFBET(B,INTAB,BB,B1,B2,EEE,P0ATB,QFATB,SINGB,BETA,JFLAG, + DERIVR) IF (SINGB) BETA = PI - BETA DTHDEA = DERIVL DTHDEB = -DERIVR IF (PR) WRITE (T21,FMT='(A,E22.14,A,E22.14)') ' ALFA=',ALFA, + ' BETA=',BETA C C CHOOSE INITIAL MATCHING POINT TMID . C IMID = 50 TMID = 0.5D0* (AA+BB) + 0.031D0*PI IF (PR) WRITE (T21,FMT='(A,E15.7,A,F11.8,A,E15.7)') ' ESTIM=',EEE, + ' TMID=',TMID IF (PR) WRITE (T21,FMT='(A,F11.8,A,F11.8,A,F11.8,A,F11.8)') + ' AAA=',AAA,' AA=',AA,' BB=',BB,' BBB=',BBB C END (SET-INITIAL-INTERVAL-AND-MATCHPOINT) C C IF EITHER ENDPOINT IS LP, USE SUBROUTINE ENDPT TO SEE IF THAT C ENDPOINT REQUIRES THE EIGENVALUES TO BE BOUNDED ABOVE, AND IF SO, C WHAT IS AN ESTIMATED UPPER BOUND. LIMA = .FALSE. IF (NCA.EQ.5 .OR. NCA.EQ.6) THEN CALL ENDPT(MF,ML,XS,PS,QS,WS,-ONE,IPA,NCA,BOUNDA) IF (PR) WRITE (T21,FMT=*) ' IPA,NCA,BOUNDA = ',IPA,NCA,BOUNDA IF (BOUNDA.LT. (10000.0D0)) THEN LIMA = .TRUE. ELIMA = BOUNDA END IF END IF LIMB = .FALSE. IF (NCB.EQ.5 .OR. NCB.EQ.6) THEN CALL ENDPT(MF,ML,XS,PS,QS,WS,ONE,IPB,NCB,BOUNDB) IF (PR) WRITE (T21,FMT=*) ' IPB,NCB,BOUNDB = ',IPB,NCB,BOUNDB IF (BOUNDB.LT. (10000.0D0)) THEN LIMB = .TRUE. ELIMB = BOUNDB END IF END IF C C IF BOTH ENDPOINTS CAUSE THE EIGENVALUES TO BE BOUNDED ABOVE, C SET ELIMUP TO BE THE LOWEST UPPER BOUND. C IF (LIMA .OR. LIMB) THEN LIMUP = .TRUE. IF (.NOT.LIMB) THEN ELIMUP = ELIMA ELSE IF (.NOT.LIMA) THEN ELIMUP = ELIMB ELSE ELIMUP = MIN(ELIMA,ELIMB) END IF IF (PR) WRITE (T21,FMT=*) + ' THE CONTINUOUS SPECTRUM HAS A LOWER ' IF (PR) WRITE (T21,FMT=*) ' BOUND, SIGMA0 = ',ELIMUP END IF C IF (EIG.EQ.0.0D0 .AND. LIMUP .AND. EEE.GE.ELIMUP) EEE = ELIMUP - + 0.01D0 CALL SETMID(MF,ML,EEE,QS,WS,IMID,TMID) C SET THE VALUES OF T, TEE(*), CORRESPONDING TO THE PLACES C WHERE THE Z's, ZEE(*), CHANGE WHEN USING SUBROUTINE GERKZ C FOR THE INTEGRATIONS FOR THETA. C THE VALUES OF JAY(*) WERE SET IN SUBROUTINE ESTPAC. DO 85 I = 1,100 TEE(I) = ONE IF (JAY(I).NE.0) TEE(I) = TFROMI(JAY(I)) IF (ZEE(I).GT.100.0D0) ZEE(I) = 100.0D0 IF (JAY(I).NE.0 .AND. PR) WRITE (T21,FMT=*) ' I,T,Z = ',I, + TEE(I),ZEE(I) 85 CONTINUE IF (JAY(3).EQ.0 .AND. ZEE(2).NE.1.0D0) THEN TEE(3) = TEE(2) ZEE(3) = ZEE(2) TEE(2) = 0.0D0 ZEE(2) = 1.0D0 ENDIF RETURN END C SUBROUTINE EIGFCN(A1,A2,B1,B2,AOK,BOK,NCA,NCB,SLFUN,ISLFUN) C ********** C THIS PROGRAM CALCULATES SELECTED EIGENFUNCTION VALUES BY C INTEGRATION OVER THE INTERNAL VARIABLE T. THE USER SELECTED C VALUES OF X IN (A,B), WHICH ARE PRESUMED TO BE IN THE ARRAY C SLFUN, ARE REPLACED BY THE COMPUTED VALUES OF THE EIGENFUNCTION. C IT IS CALLED FROM SLEIGN. C INPUT QUANTITIES: A1,A2,B1,B2,AOK,BOK,NCA,NCB, C SLFUN,ISLFUN C OUTPUT QUANTITIES: SLFUN C C N.B.: IN THIS PROGRAM IT IS ASSUMED THAT THE POINTS T C (CORRESPONDING TO THE GIVEN VALUES OF X IN SLFUN) C ALL LIE WITHIN THE INTERVAL (AA,BB).----- C C .. Scalars in Common .. REAL AA,BB,DTHDAA,DTHDBB,TMID INTEGER MDTHZ,T21 LOGICAL ADDD,PR C .. C .. Local Scalars .. REAL DTHDAT,DTHDBT,DTHDET,EFF,T,THT,TM INTEGER I,IFLAG,J,NC,NMID LOGICAL LCOA,LCOB,OK C .. C .. Local Arrays .. REAL ERL(3),ERR(3),YL(3),YR(3) C .. C .. External Subroutines .. EXTERNAL INTEG C C .. External Functions .. REAL TFROMX EXTERNAL TFROMX C .. C .. Intrinsic Functions .. INTRINSIC EXP,SIN C .. C .. Common blocks .. COMMON /PRIN/PR,T21 COMMON /TDATA/AA,TMID,BB,DTHDAA,DTHDBB,MDTHZ,ADDD C .. C .. Scalar Arguments .. REAL A1,A2,B1,B2,TI,XI INTEGER ISLFUN,NCA,NCB LOGICAL AOK,BOK C .. C .. Array Arguments .. REAL SLFUN(ISLFUN+9) C .. DO 5 I = 1, ISLFUN XI = SLFUN(9+I) TI = TFROMX(XI) SLFUN(9+I) = TI 5 CONTINUE c .. NMID = 0 DO 10 I = 1,ISLFUN IF (SLFUN(9+I).LE.TMID) NMID = I 10 CONTINUE IF (NMID.GT.0) THEN T = AA YL(1) = SLFUN(3) YL(2) = 0.0D0 YL(3) = 0.0D0 LCOA = NCA .EQ. 4 OK = AOK NC = NCA EFF = 0.0D0 DO 20 J = 1,NMID TM = SLFUN(J+9) IF (TM.LT.AA .OR. TM.GT.BB) THEN IF (PR) WRITE (*,FMT=*) ' T.LT.AA .OR. T.GT.BB ' GOTO 20 END IF THT = YL(1) DTHDAT = DTHDAA*EXP(-2.0D0*EFF) DTHDET = YL(2) IF (TM.GT.AA) THEN CALL INTEG(T,THT,DTHDAT,DTHDET,TM,A1,A2,SLFUN(8),YL, + ERL,OK,NC,IFLAG) IF (LCOA) THEN EFF = YL(3) ELSE NC = 1 OK = .TRUE. EFF = EFF + YL(3) END IF END IF SLFUN(J+9) = SIN(YL(1))*EXP(EFF+SLFUN(4)) T = TM IF (T.GT.-1.0D0) NC = 1 IF (T.LT.-0.9D0 .AND. LCOA) THEN NC = 4 T = AA YL(1) = SLFUN(3) YL(2) = 0.0D0 YL(3) = 0.0D0 END IF 20 CONTINUE END IF IF (NMID.LT.ISLFUN) THEN T = BB YR(1) = SLFUN(6) YR(2) = 0.0D0 YR(3) = 0.0D0 LCOB = NCB .EQ. 4 OK = BOK NC = NCB EFF = 0.0D0 DO 30 J = ISLFUN,NMID + 1,-1 TM = SLFUN(J+9) IF (TM.LT.AA .OR. TM.GT.BB) THEN IF (PR) WRITE (*,FMT=*) ' T.LT.AA .OR. T.GT.BB ' GOTO 30 END IF THT = YR(1) DTHDBT = DTHDBB*EXP(-2.0D0*EFF) DTHDET = YR(2) IF (TM.LT.BB) THEN CALL INTEG(T,THT,DTHDBT,DTHDET,TM,B1,B2,SLFUN(8),YR, + ERR,OK,NC,IFLAG) IF (LCOB) THEN EFF = YR(3) ELSE OK = .TRUE. NC = 1 EFF = EFF + YR(3) END IF END IF SLFUN(J+9) = SIN(YR(1))*EXP(EFF+SLFUN(7)) IF (ADDD) SLFUN(J+9) = -SLFUN(J+9) T = TM IF (T.LT.1.0D0) NC = 1 IF (T.GT.0.9D0 .AND. LCOB) THEN NC = 4 T = BB YR(1) = SLFUN(6) YR(2) = 0.0D0 YR(3) = 0.0D0 END IF 30 CONTINUE END IF RETURN END C SUBROUTINE EFDATA(ALFA,DTHDEA,A1,A2,BETA,DTHDEB,B1,B2,EIG,EIGPI, + EPS,AOK,NCA,BOK,NCB,RAY,SLFUN) C C THIS PROGRAM IS USED ONLY AFTER THE WANTED EIGENVALUE HAS BEEN C SUCCESSFULLY COMPUTED. IT IS CALLED FROM SLEIGN. C IT CONVERTS THE T-DATA AA,TMID,BB TO CORRESPONDING X-DATA, C FILLS 7 OF THE FIRST 9 LOCATIONS OF SLFUN, C COMPUTES VALUES OF LOG(RHO(A)), LOG(RHO(B)) SUCH THAT THE C CORRESPONDING SOLUTION Y(X) = RHO(X)*SIN(THETA(X)) IS C CONTINUOUS AT XMID, AND HAS L2-NORM EQUAL TO 1.0 . C THESE VALUES ARE PLACED IN SLFUN(4) AND SLFUN(7). C THE COMPUTED VALUE OF EIG IS FINALLY CHECKED ONE LAST TIME C USING THE JUMPS IN Y AND PY' AT XMID IN A FORM OF C RALEIGH QUOTIENT CORRECTION. (HERE CALLED "RAY".) C C INPUT QUANTITIES: ALFA,DTHDEA,A1,A2,BETA,DTHDEB,B1,B2,EIG, C EIGPI,EPS,AOK,BOK,NCA,NCB,SLFUN C OUTPUT QUANTITIES: XAA,XMID,XBB,EIGSAV,ISAVE,AA,BB,RAY,ADDD, C SLFUN C C .. Scalar Arguments .. REAL A1,A2,ALFA,B1,B2,BETA,DTHDEA,DTHDEB,EIG,EIGPI, + EPS,RAY INTEGER NCA,NCB LOGICAL AOK,BOK C .. C .. Array Arguments .. REAL SLFUN(9) C .. C .. Scalars in Common .. REAL AA,BB,DTHDAA,DTHDBB,EIGSAV,TMID,TSAVEL,TSAVER,Z INTEGER IND,ISAVE,MDTHZ,T21 LOGICAL ADDD,PR C .. C .. Local Scalars .. REAL AAF,BBF,CL,CR,DEN,DUM,E,ONE,PSIL,PSIPL,PSIPR, + PSIR,SL,SQL,SQR,SR,THDAAX,THDBBX,TMP,UL,UR,XAA, + XBB,XMID INTEGER JFLAG C .. C .. Local Arrays .. REAL ERL(3),ERR(3),YL(3),YR(3) C .. C .. External Subroutines .. EXTERNAL AABB,DXDT,INTEG C .. C .. Intrinsic Functions .. INTRINSIC ABS,COS,EXP,LOG,SIN,SQRT C .. C .. Common blocks .. COMMON /DATAF/EIGSAV,IND COMMON /PRIN/PR,T21 COMMON /TDATA/AA,TMID,BB,DTHDAA,DTHDBB,MDTHZ,ADDD COMMON /TSAVE/TSAVEL,TSAVER,ISAVE COMMON /Z1/Z C .. ONE = 1.0D0 Z = 1.0D0 50 CONTINUE CALL DXDT(TMID,TMP,XMID) CALL DXDT(AA,TMP,XAA) SLFUN(1) = XMID SLFUN(2) = XAA SLFUN(3) = ALFA SLFUN(6) = BETA + EIGPI SLFUN(8) = EPS SLFUN(9) = Z C C INTEGRATE FROM BOTH ENDPOINTS TO THE MIDPOINT: C EIGSAV = EIG THDAAX = 0.0D0 YL(1) = 0.0D0 YL(2) = 0.0D0 YL(3) = 0.0D0 ISAVE = 0 CALL INTEG(AA,ALFA,THDAAX,DTHDEA,TMID,A1,A2,EPS,YL,ERL,AOK,NCA, + JFLAG) IF (JFLAG.EQ.5) THEN AAF = AA CALL AABB(AA,-ONE) IF (PR) WRITE (T21,FMT=*) ' AA MOVED FROM ',AAF,'IN TO',AA GO TO 50 END IF 100 CONTINUE CALL DXDT(BB,TMP,XBB) SLFUN(5) = XBB THDBBX = 0.0D0 YR(1) = 0.0D0 YR(2) = 0.0D0 YR(3) = 0.0D0 ISAVE = 0 CALL INTEG(BB,BETA,THDBBX,DTHDEB,TMID,B1,B2,EPS,YR,ERR,BOK,NCB, + JFLAG) IF (JFLAG.EQ.5) THEN BBF = BB CALL AABB(BB,-ONE) IF (PR) WRITE (T21,FMT=*) ' BB MOVED FROM ',BBF,'IN TO',BB GO TO 100 END IF YR(1) = YR(1) + EIGPI SL = SIN(YL(1)) SR = SIN(YR(1)) CL = COS(YL(1)) CR = COS(YR(1)) C UL AND UR ARE THE VALUES OF THE INTEGRAL OF Y**2*W FROM THE C ENDPOINTS A AND B, RESPECTIVELY. UL = (YL(2)-DTHDEA*EXP(-2.0D0*YL(3))) UR = (YR(2)-DTHDEB*EXP(-2.0D0*YR(3))) DEN = 0.5D0*LOG(UL*SR*SR-UR*SL*SL) DUM = 0.5D0*LOG(UL-UR) C COMPUTE SLFUN(4), SLFUN(7) TOWARDS NORMALIZING THE EIGENFUNCTION. SLFUN(4) = -YL(3) - DUM SLFUN(7) = -YR(3) - DUM C DO (CHECK-MATCHING-VALUES-OF-EIGENFUNCTION) C PERFORM FINAL CHECK ON EIG. C DEN = UL*SR*SR - UR*SL*SL E = ABS(SR)/SQRT(DEN) PSIL = E*SL PSIPL = E*CL SQL = E*E*UL E = ABS(SL)/SQRT(DEN) PSIR = E*SR PSIPR = E*CR SQR = E*E*UR ADDD = PSIL*PSIR .LT. 0.0D0 .AND. PSIPL*PSIPR .LT. 0.0D0 RAY = EIG + (PSIL*PSIPL-PSIR*PSIPR)/ (SQL-SQR) IF (PR) WRITE (T21,FMT=*) ' RAY,ADDD = ',RAY,ADDD C END (CHECK-MATCHING-VALUES-OF-EIGENFUNCTION) RETURN END C SUBROUTINE SETMID(MF,ML,EIG,QS,WS,IMID,TMID) C ********** C C THIS PROGRAM SELECTS A POINT IN (-1,1) AS TMID. C IT IS CALLED BY RESET AND PRELIM. C IT TESTS THE INTERVAL SAMPLE POINTS IN THE ORDER C 50,51,49,52,48,...,ETC. FOR THE FIRST ONE WHERE THE EXPRESSION C (LAMBDA*W-Q) IS POSITIVE. THIS T-POINT IS DESIGNATED TMID. C THEORETICALLY IT DOESN'T MATTER WHICH POINT T IN (-1,1) IS USED C FOR TMID, BUT IN PRACTICE IT CAN MAKE A GREAT DEAL OF C DIFFERENCE. THIS SCHEME IS JUST A RULE OF THUMB WHICH SEEMS C TO WORK FAIRLY WELL IN AVOIDING BAD CHOICES FOR TMID. C C INPUT QUANTITIES: MF, ML, EIG, QS, WS C OUTPUT QUANTITIES: IMID, TMID C ********** C .. Local Scalars .. REAL S INTEGER I,J C .. C .. External Functions .. REAL TFROMI EXTERNAL TFROMI C .. C .. Scalar Arguments .. REAL EIG,TMID INTEGER IMID,MF,ML C .. C .. Array Arguments .. REAL QS(*),WS(*) C .. C .. Scalars in Common .. INTEGER T21 LOGICAL PR C .. C .. Common blocks .. COMMON /PRIN/PR,T21 C .. S = -1.0D0 DO 10 J = 1,100 I = 50 + S* (J/2) S = -S IF (I.LT.MF .OR. I.GT.ML) GO TO 20 IF (EIG*WS(I)-QS(I).GT.0.0D0) THEN IMID = I TMID = TFROMI(IMID) GO TO 20 END IF 10 CONTINUE 20 CONTINUE IF (PR) WRITE (T21,FMT=*) ' NEW TMID = ',TMID RETURN END C SUBROUTINE EXTR(MM,F,GQ,GW,XS,PS,QS,WS) C THIS SUBROUTINE IS CALLED FROM SUBROUTINE ENDPT ONLY, C FOR THE PURPOSE OF DETERMINING THE ENDPOINT VALUES C OF F (OR Q, OR W) IN THE INDICIAL EQUATION AT A FINITE C LP ENDPOINT. C IT PUTS VALUES OF F (OR Q OR W), CORRESPONDING TO C EQUALLY SPACED VALUES OF X, IN C AN ARRAY DF, AND EXTRAPOLATES THE VALUES OF F (OR Q, OR W) C TO THE ENDPOINT, XEND, OF THE INTERVAL. C C INPUT QUANTITIES: MM,XS,PS,QS,WS,A,B,EPSMIN C OUTPUT QUANTITIES: F,GQ,GW C .. Scalar Arguments .. REAL F,GQ,GW INTEGER MM C .. C .. Array Arguments .. REAL PS(100),QS(100),WS(100),XS(100) C .. C .. Scalars in Common .. REAL A,B,EPSMIN INTEGER INTAB C .. C .. Local Scalars .. REAL TI,TIE,TMP,XEND,XI,XIE,EPST,ERR INTEGER EP,I,J C .. C .. Local Arrays .. REAL DF(6),DQ(6),DW(6),XX(6) C .. C .. External Functions .. REAL P,TFROMI EXTERNAL P,TFROMI C .. C .. External Subroutines .. EXTERNAL DXDT,INTPOL C .. C .. Common blocks .. COMMON /DATADT/A,B,INTAB COMMON /RNDOFF/EPSMIN C .. IF (MM.LT.50) THEN EP = 1 XEND = A ELSE EP = -1 XEND = B END IF C HERE WE WANT TO FIND THE LIMIT OF C (X-A)*P'(X)/P(X) AT X=A C AND OF C (X-A)**2*Q(X)/P(X) C AND OF C (X-A)**2*W(X)/P(X) C OR THE SAME SORT OF THINGS AT X=B. DO 10 J = 1,5 I = MM + (J-1)*EP TI = TFROMI(I) TIE = TI + 2.0D0*EPSMIN*EP CALL DXDT(TI,TMP,XI) CALL DXDT(TIE,TMP,XIE) DF(J) = (XS(I)-XEND)* (P(XIE)-PS(I))/ ((XIE-XS(I))*PS(I)) DQ(J) = (XS(I)-XEND)**2*QS(I) DW(J) = (XS(I)-XEND)**2*WS(I) XX(J) = XS(I) 10 CONTINUE EPST = 0.00001D0 CALL INTPOL(5,XX,DF,XEND,EPST,3,F,ERR) CALL INTPOL(5,XX,DQ,XEND,EPST,3,GQ,ERR) CALL INTPOL(5,XX,DW,XEND,EPST,3,GW,ERR) IF(ABS(F).LE.ERR) F = 0.0D0 IF(ABS(GQ).LE.ERR) GQ = 0.0D0 IF(GW.LE.ERR) GW = 0.0D0 RETURN END C SUBROUTINE ENDPT(MF,ML,XS,PS,QS,WS,TEND,IPM,NC,BOUND) C THIS PROGRAM IS CALLED FROM PRELIM. C MAINLY IT COMPUTES THE QUANTITIES IN /COMMON/ALBE/ FOR C USE WHEN A FINITE ENDPOINT IS LP. C IF THE ENDPOINT IS SUCH THAT THE QUANTITIES INVOLVED C APPEAR TO NOT HAVE LIMITING VALUES, (AS HAPPENS IN THE C PROBLEM C p(x) = 1, q(x) = cos(x)**2, w(x) = 1 on (0,+Inf) C FOR EXAMPLE), THEN SLEIGN2 CANNOT CONTINUE. C C AT THE SAME TIME, IT TRIES TO DETERMINE WHETHER OR NOT C THE EIGENVALUES HAVE A FINITE UPPER BOUND. C C INPUT QUANTITIES: MF,ML,XS,PS,QS,WS,TEND,INTAB,EPSMIN C OUTPUT QUANTITIES: LPWA,LPQA,LPWB,LPQB,NC,BOUND C C RECALL THAT THE STORED ARRAYS QS AND WS ARE REALLY THE C VALUES OF Q/P AND W/P. C N.B.: INDICIAL EQUATION IS (AT A REGULAR SINGULAR POINT): C S*S + (F-1)*S + G = 0 C WHERE F IS FA OR FB C AND G IS LAMBDA*GWA - GQA C OR LAMBDA*GWB - GQB C (HERE, THE INDEPENDENT VARIABLE IS X IN (A,B) C .. Scalar Arguments .. REAL BOUND,TEND INTEGER IPM,MF,ML,NC C .. C .. Array Arguments .. REAL PS(100),QS(100),WS(100),XS(100) C .. C .. Scalars in Common .. REAL A,B,EPSMIN,LPWA,LPQA,FA,GWA,GQA,LPWB, + LPQB,FB,GWB,GQB INTEGER INTAB,T21 LOGICAL PR C .. C .. Local Scalars .. REAL APQ1,APQ2,APW1,APW2,PQ1,PQ2,PW1,PW2,TMP INTEGER I,IQ,IW,MM,MM6 LOGICAL LOGIC C .. C .. External Subroutines .. EXTERNAL EXTR C .. C .. Intrinsic Functions .. INTRINSIC ABS,MIN C .. C .. Common blocks .. COMMON /ALBE/LPWA,LPQA,FA,GWA,GQA,LPWB,LPQB,FB,GWB,GQB COMMON /DATADT/A,B,INTAB COMMON /PRIN/PR,T21 COMMON /RNDOFF/EPSMIN C .. IF (PR) WRITE (T21,FMT=*) ' IN ENDPT, NC = ',NC IF (TEND.LT.0.0D0) THEN MM = MF ELSE MM = ML END IF BOUND = 1.D+19 IQ = 0 IW = 0 IF (MM.LT.50) THEN MM6 = MM + 6 DO 10 I = MM,MM6 PW1 = PS(I+1)**2*WS(I+1) - PS(I)**2*WS(I) PW2 = PS(I+2)**2*WS(I+2) - PS(I+1)**2*WS(I+1) PQ1 = PS(I+1)**2*QS(I+1) - PS(I)**2*QS(I) PQ2 = PS(I+2)**2*QS(I+2) - PS(I+1)**2*QS(I+1) APW1 = ABS(PW1) APW2 = ABS(PW2) APQ1 = ABS(PQ1) APQ2 = ABS(PQ2) IF ((APW1.LE.EPSMIN.AND.APW2.LE.EPSMIN) .OR. + PW1*PW2.GT.0.D0) IW = IW + 1 IF ((APQ1.LE.EPSMIN.AND.APQ2.LE.EPSMIN) .OR. + PQ1*PQ2.GT.0.D0) IQ = IQ + 1 10 CONTINUE IPM = MIN(IW,IQ) IF (IPM.GE.5) THEN LPWA = PS(MM)**2*WS(MM) LPQA = PS(MM)**2*QS(MM) IF (INTAB.EQ.1 .OR. INTAB.EQ.2) THEN CALL EXTR(MM,FA,GQA,GWA,XS,PS,QS,WS) IF (PR) WRITE (T21,FMT=*) ' FA,GWA,GQA = ',FA,GWA,GQA LOGIC = (ABS(FA)+ABS(GWA)+ABS(GQA)) .LT. 1000.0D0 IF (LOGIC) THEN C INDICIAL EQUATION IS: S*S + (FA-1)*S + (LAMBDA*GWA-GQA)=0 IF (GWA.GT.EPSMIN) THEN C FOR REAL ROOTS, REQUIRE LAMBDA .LE. BOUND, WHERE: TMP = (0.25D0*(FA-1.0D0)**2+GQA)/GWA BOUND = MIN(BOUND,TMP) END IF ELSE C THIS CASE MEANS THAT THE ENDPOINT IS PROBABLY AN "IRREGULAR" C POINT. NC = 6 END IF ELSE C THIS CASE MEANS THAT THE ENDPOINT A IS AT -INF. SO IF (LPWA.GT.EPSMIN) BOUND = MIN(BOUND,LPQA/LPWA) C CHOOSE NC=7 ONLY IF THE USER HAS NOT ALREADY CHOSEN NC=6. IF (NC.NE.6) NC = 7 C TEMPORARILY SET (NC = 7 IS NOT BEING USED YET) NC = 6 END IF ELSE C ARRIVING HERE MEANS THAT IPM.LT.7, WHICH MEANS THAT C P*W AND P*Q DO NOT APPEAR TO BE MONOTONELY TENDING C TO A LIMIT (OR +INF OR -INF). I SUPPOSE THIS MAY C MEAN THAT THERE ARE BANDS & GAPS. I DON'T KNOW. IF (PR) WRITE (T21,FMT=*) + ' P*W & P*Q BEHAVE BADLY NEAR THE END A. ' NC = 8 END IF ELSE MM6 = MM - 6 DO 20 I = MM,MM6,-1 PW1 = PS(I-1)**2*WS(I-1) - PS(I)**2*WS(I) PW2 = PS(I-2)**2*WS(I-2) - PS(I-1)**2*WS(I-1) PQ1 = PS(I-1)**2*QS(I-1) - PS(I)**2*QS(I) PQ2 = PS(I-2)**2*QS(I-2) - PS(I-1)**2*QS(I-1) APW1 = ABS(PW1) APW2 = ABS(PW2) APQ1 = ABS(PQ1) APQ2 = ABS(PQ2) IF ((APW1.LE.EPSMIN.AND.APW2.LE.EPSMIN) .OR. + PW1*PW2.GT.0.D0) IW = IW + 1 IF ((APQ1.LE.EPSMIN.AND.APQ2.LE.EPSMIN) .OR. + PQ1*PQ2.GT.0.D0) IQ = IQ + 1 20 CONTINUE IPM = MIN(IQ,IW) IF (IPM.GE.5) THEN LPWB = PS(MM)**2*WS(MM) LPQB = PS(MM)**2*QS(MM) IF (INTAB.EQ.1 .OR. INTAB.EQ.3) THEN CALL EXTR(MM,FB,GQB,GWB,XS,PS,QS,WS) IF (PR) WRITE (T21,FMT=*) ' FB,GWB,GQB = ',FB,GWB,GQB LOGIC = (ABS(FB)+ABS(GWB)+ABS(GQB)) .LT. 1000.0D0 IF (LOGIC) THEN C INDICIAL EQUATION IS: S*S + (FB-1)*S + (LAMBDA*GWB-GQB)=0 IF (GWB.GT.EPSMIN) THEN C FOR REAL ROOTS, REQUIRE LAMBDA .LE. BOUND, WHERE: TMP = (0.25D0*(FB-1.0D0)**2+GQB)/GWB BOUND = MIN(BOUND, TMP) END IF ELSE C THIS CASE MEANS THAT THE ENDPOINT IS PROBABLY AN "IRREGULAR" C POINT. NC = 6 END IF ELSE C THIS CASE MEANS THAT THE ENDPOINT B IS AT +INF. SO IF (LPWB.GT.EPSMIN) BOUND = MIN(BOUND,LPQB/LPWB) C CHOOSE NC=7 ONLY IF THE USER HAS NOT ALREADY CHOSEN NC=6. IF (NC.NE.6) NC = 7 C TEMPORARILY SET (NC = 7 IS NOT BEING USED YET) NC = 6 END IF ELSE C ARRIVING HERE MEANS THAT IPM.LT.7, WHICH MEANS THAT C P*W AND P*Q DO NOT APPEAR TO BE MONOTONELY TENDING C TO A LIMIT (OR +INF OR -INF). I SUPPOSE THIS MAY C MEAN THAT THERE ARE BANDS & GAPS. I DON'T KNOW. IF (PR) WRITE (T21,FMT=*) + ' P*W & P*Q BEHAVE BADLY NEAR THE END B. ' NC = 8 END IF END IF RETURN END C SUBROUTINE THZ2TH(U,ERU,Z,Y,ERY) C ********** C THIS PROGRAM IS CALLED FROM GERKZ. THERE THE PRUEFER C ANGLE HAS A CONSTANT Z IN ITS DEFINITION, AND IS HERE C DENOTED BY THZ. THE USUAL THETA (EQUIVALENT TO Z = 1) C IS DENOTED BY TH. C THIS PROGRAM CONVERTS FROM THZ TO TH WHERE C TAN(TH)=Y/(P*Y') C AND C TAN(THZ)=Z*Y/(P*Y'), C OR C TAN(THZ)=Z*TAN(TH) . C SO WE HAVE C DTHZ=Z*(COS(THZ)/COS(TH))**2 * DTH , C OR C DTH=(1/Z)*(COS(TH)/COS(THZ))**2 * DTHZ . C C INPUT QUANTITIES: U,ERU,Z C OUTPUT QUANTITIES: Y,ERY C ********** C .. Scalars in Common .. REAL HPI,PI,TWOPI C .. C .. Local Scalars .. REAL DTH,DTHZ,DUM,FAC,PIK,REMTHZ,TH,THZ INTEGER K C .. C .. Intrinsic Functions .. INTRINSIC ABS,ATAN,COS,LOG,SIN,TAN C .. C .. Common blocks .. COMMON /PIE/PI,TWOPI,HPI C .. C .. Scalar Arguments .. REAL Z C .. C .. Array Arguments .. REAL ERU(3),ERY(3),U(3),Y(3) C .. THZ = U(1) DTHZ = U(2) K = THZ/PI IF (THZ.LT.0.0D0) K = K - 1 PIK = K*PI REMTHZ = THZ - PIK IF (4.0D0*REMTHZ.LE.PI) THEN TH = ATAN(TAN(REMTHZ)/Z) + PIK ELSE IF (4.0D0*REMTHZ.GE.3.0D0*PI) THEN TH = ATAN(TAN(REMTHZ)/Z) + PIK + PI ELSE TH = ATAN(Z*TAN(REMTHZ-HPI)) + PIK + HPI END IF C THE TWO DIFFERENT APPEARING FORMULAS BELOW FOR FAC ARE C IN FACT EQUIVALENT. WE USE WHICHEVER ONE AVOIDS C DIVIDING BY A SMALL NUMBER. DUM = ABS(COS(TH)) IF (DUM.GE.0.5D0) THEN FAC = (DUM/COS(THZ))**2/Z ELSE FAC = Z* (SIN(TH)/SIN(THZ))**2 END IF DTH = FAC*DTHZ Y(1) = TH Y(2) = DTH Y(3) = U(3) + 0.5D0*LOG(Z/FAC) C ALSO CONVERT THE ESTIMATED ERRORS IN THE THZ COMPUTATIONS C TO CORRESPONDING ERRORS IN TH. ERY(1) = FAC*ERU(1) ERY(2) = FAC*ERU(2) ERY(3) = FAC*ERU(3) RETURN END C SUBROUTINE TH2THZ(Y,Z,U) C ********** C THIS PROGRAM IS CALLED FROM GERKZ. THERE THE PRUEFER C ANGLE HAS A CONSTANT Z IN ITS DEFINITION, AND IS HERE C DENOTED BY THZ. THE USUAL THETA (EQUIVALENT TO Z = 1) C IS DENOTED BY TH. C THIS PROGRAM CONVERTS FROM TH TO THZ WHERE C TAN(TH)=Y/(P*Y') C AND C TAN(THZ)=Z*Y/(P*Y'), C OR C TAN(THZ)=Z*TAN(TH) . C SO WE HAVE C DTHZ=Z*(COS(THZ)/COS(TH))**2 * DTH , C OR C DTH=(1/Z)*(COS(TH)/COS(THZ))**2 * DTHZ . C C INPUT QUANTITIES: Y,Z C OUTPUT QUANTITIES: U C ********** C .. Scalars in Common .. REAL HPI,PI,TWOPI C .. C .. Local Scalars .. REAL DTH,DTHZ,DUM,FAC,PIK,REMTH,TH,THZ INTEGER K C .. C .. Intrinsic Functions .. INTRINSIC ABS,ATAN,COS,LOG,SIN,TAN C .. C .. Common blocks .. COMMON /PIE/PI,TWOPI,HPI C .. C .. Scalar Arguments .. REAL Z C .. C .. Array Arguments .. REAL U(3),Y(3) C .. TH = Y(1) DTH = Y(2) K = TH/PI IF (TH.LT.0.0D0) K = K - 1 PIK = K*PI REMTH = TH - PIK IF (4.0D0*REMTH.LE.PI) THEN THZ = ATAN(Z*TAN(REMTH)) + PIK ELSE IF (4.0D0*REMTH.GE.3.0D0*PI) THEN THZ = ATAN(Z*TAN(REMTH)) + PIK + PI ELSE THZ = ATAN(TAN(REMTH-HPI)/Z) + PIK + HPI END IF C THE TWO DIFFERENT APPEARING FORMULAS BELOW FOR FAC ARE C IN FACT EQUIVALENT. WE USE WHICHEVER ONE AVOIDS C DIVIDING BY A SMALL NUMBER. DUM = ABS(COS(THZ)) IF (DUM.GE.0.5D0) THEN FAC = Z* (DUM/COS(TH))**2 ELSE FAC = (SIN(THZ)/SIN(TH))**2/Z END IF DTHZ = FAC*DTH U(1) = THZ U(2) = DTHZ U(3) = Y(3) - 0.5D0*LOG(Z*FAC) RETURN END C SUBROUTINE PQEXT(F) C THIS SUBROUTINE IS CALLED ONLY BY SUBROUTINE WR. C IT IS USED TO EXTRAPOLATE THE VALUES F(I), C I=2,5, TO F(1) . IT ASSUMES THE INDEPENDENT C VARIABLE VALUES ARE EQUALLY SPACED. C C INPUT QUANTITIES: F(2),F(3),F(4),F(5) C OUTPUT QUANTITIES: F(1) C .. Array Arguments .. REAL F(6) C .. C .. Local Scalars .. REAL D2F1,D2F2,D2F3,D3F2 INTEGER I C .. C .. Local Arrays .. REAL DF(6) C .. DO 10 I = 2,4 DF(I) = F(I+1) - F(I) 10 CONTINUE D2F3 = DF(4) - DF(3) D2F2 = DF(3) - DF(2) D3F2 = D2F3 - D2F2 D2F1 = D2F2 + D3F2 DF(1) = DF(2) + D2F1 F(1) = F(2) + DF(1) RETURN END C SUBROUTINE FIT(TH1,TH,TH2) C ********** C C THIS PROGRAM IS CALLED ONLY FROM SUBROUTINE UVPHI, WHICH C IS CALLED BY SUBROUTINE LCO. C IT CONVERTS TH INTO AN 'EQUIVALENT' ANGLE BETWEEN C TH1 AND TH2. WE ASSUME TH1.LT.TH2 AND PI.LE.(TH2-TH1). C C INPUT QUANTITIES: TH1,TH2 C OUTPUT QUANTITIES: TH C C ********** C .. Scalars in Common .. REAL HPI,PI,TWOPI C .. C .. Intrinsic Functions .. INTRINSIC AINT C .. C .. Common blocks .. COMMON /PIE/PI,TWOPI,HPI C .. C .. Scalar Arguments .. REAL TH,TH1,TH2 C .. IF (TH.LT.TH1) TH = TH + AINT((TH1-TH+PI)/PI)*PI IF (TH.GT.TH2) TH = TH - AINT((TH-TH2+PI)/PI)*PI RETURN END C SUBROUTINE F(U,Y,YP) C ********** C C THIS SUBROUTINE EVALUATES THE DERIVATIVE FUNCTIONS FOR USE WITH C INTEGRATOR GERK IN SOLVING THE SYSTEM OF DIFFERENTIAL C EQUATIONS FOR THETA, RHO (OR, RATHER, LOG(RHO)), AND C DTHDE = D(THETA)/D(LAMBDA), WHERE C Y = RHO*SIN(THETA) C AND C PY' = Z*RHO*COS(THETA) . C THE DIFFERENTIAL EQUATIONS ARE: C THETA' = (Z/P)*COS(THETA)**2 + (EIG*W - Q)*SIN(THETA)**2/Z C LOG(RHO)' = (Z/P - (EIG*W - Q)/Z)*SIN(THETA)*COS(THETA) C DTHDE' = -2*(Z/P - (EIG*W - Q)/Z)*DTHDE + W*SIN(THETA)**2/Z C C EXCEPT WHEN CALLED FROM GERKZ, Z IS ALWAYS = 1.0 . C C INPUT QUANTITIES: U,Y,EIG,IND,Z C OUTPUT QUANTITIES: YP C ********** C .. Scalars in Common .. REAL EIG,Z INTEGER IND C .. C .. Local Scalars .. REAL C,C2,DT,QX,S,S2,T,TH,V,WW,WX,X,XP C .. C .. External Functions .. REAL P,Q,W EXTERNAL P,Q,W C .. C .. External Subroutines .. EXTERNAL DXDT C .. C .. Intrinsic Functions .. INTRINSIC COS,MOD,SIN C .. C .. Common blocks .. COMMON /DATAF/EIG,IND COMMON /Z1/Z C .. C .. Scalar Arguments .. REAL U C .. C .. Array Arguments .. REAL Y(2),YP(3) C .. IF (MOD(IND,2).EQ.1) THEN T = U TH = Y(1) ELSE T = Y(1) TH = U END IF CALL DXDT(T,DT,X) XP = Z/P(X) QX = Q(X)/Z WX = W(X)/Z V = EIG*WX - QX S = SIN(TH) C = COS(TH) S2 = S*S C2 = C*C YP(1) = DT* (XP*C2+V*S2) IF (IND.EQ.1) THEN WW = (XP-V)*S*C YP(2) = DT* (-2.0D0*WW*Y(2)+WX*S2) YP(3) = DT*WW ELSE IF (IND.EQ.2) THEN C IN THIS CASE THE INDEPENDENT AND DEPENDENT VARIABLES C (t and dtheta, etc.) HAVE BEEN INTERCHANGED. YP(2) = YP(2)/YP(1) YP(3) = YP(3)/YP(1) YP(1) = 1.0D0/YP(1) ELSE IF (IND.EQ.3) THEN ELSE YP(1) = 1.0D0/YP(1) END IF RETURN END C SUBROUTINE UVPHI(U,PUP,V,PVP,THU,THV,PHI,TH) C ********** C C THIS PROGRAM IS CALLED BY SUBROUTINE LCO. C IT FINDS TH (= THETA) APPROPRIATE TO THU, THV, AND PHI, WHERE C THU IS THE PHASE ANGLE FOR U, AND THV IS THE PHASE ANGLE FOR V. C C INPUT QUANTITIES: U,PUP,V,PVP,THU,THV C OUTPUT QUANTITIES: TH C C ********** C .. Scalars in Common .. REAL HPI,PI,TWOPI C .. C .. Local Scalars .. REAL C,D,PYP,S,Y C .. C .. External Subroutines .. EXTERNAL FIT C .. C .. Intrinsic Functions .. INTRINSIC ATAN2,COS,SIN C .. C .. Common blocks .. COMMON /PIE/PI,TWOPI,HPI C .. C .. Scalar Arguments .. REAL PHI,PUP,PVP,TH,THU,THV,U,V C .. TH = THU IF (PHI.EQ.0.0D0) RETURN IF (THV-THU.LT.PI) THEN TH = THV IF (PHI.EQ.-HPI) RETURN TH = THV - PI IF (PHI.EQ.HPI) RETURN ELSE TH = THV - PI IF (PHI.EQ.-HPI) RETURN TH = THV - TWOPI IF (PHI.EQ.HPI) RETURN END IF C = COS(PHI) S = SIN(PHI) Y = U*C + V*S PYP = PUP*C + PVP*S TH = ATAN2(Y,PYP) IF (Y.LT.0.0D0) TH = TH + TWOPI D = U*PVP - V*PUP IF (D*PHI.GT.0.0D0) THEN CALL FIT(THU-PI,TH,THU) ELSE CALL FIT(THU,TH,THU+PI) END IF RETURN END SUBROUTINE WR(FG,NC,TSTAR,YSTAR,THEND,DTHDE,TOUT,Y,TT,YY,IFLAG, + ERR,WORK,IWORK) C ********** C C THIS PROGRAM IS CALLED BY WREG AND LCNO. C IT INTEGRATES Y' = F(T,Y) FROM T = +/-1.0 TO THEND C WHEN F CANNOT BE EVALUATED AT T. (T*,Y*) IS CHOSEN AS A C NEARBY POINT, AND THE EQUATION IS INTEGRATED FROM THERE AND C CHECKED FOR CONSISTENCY WITH HAVING INTEGRATED FROM T. IF NOT, C A DIFFERENT (T*,Y*) IS CHOSEN UNTIL CONSISTENCY IS ACHIEVED. C C INPUT QUANTITIES: NC,TSTAR,YSTAR,THEND,DTHDE C OUTPUT QUANTITIES: IFLAG,TT,YY,IND,Y,TOUT,ERR,WORK,IWORK,TSTAR C C The criterion to be used for "consistency" here is a so-called C "correction formula", obtained by integrating Newton's backward C difference formula - which can be written as C Dx0 = h(q4 - (7/2)Dq3 + (53/12)DDq2 - (55/24)DDDq1 + C + (251/720)DDDDq0 ....) C (Here Dx0 = x1-x0, Dq1 = q2-q1, etc., and the qi terms are the C derivatives of the function x(*) at the points yi.) When the C numerical integration for x has reached y4 and there is some C uncertainty about the earlier value of x1 at y1, this formula C can be used to effect a correction. C C THE BASIC IDEA IS THIS: C SUPPOSE THE PROBLEM IS TO INTEGRATE C y'(x) = f(x, y(x)), y(a) = A C FROM a TO c, BUT f(a,A) is +Infinity. C INTERCHANGE INDEPENDENT AND DEPENDENT VARIABLES SO THAT THE C PROBLEM BECOMES C x'(y) = 1/f(x,y), x(A) = a. C Let yi, i=1,2,3,4, be equally spaced points, distance h apart. C Let x1 be an initial "guess" for the value of the solution at y1, C and integrate the d.e. for the variable x from y1 to y4, obtaining C values xi for i=1,2,3,4. Let qi be the values for the derivative C function for x at the points yi. Then, according to the above C correction formula, a "corrected" value of the initial guess x1 is C x1 = x0 + h(q4 - (7/2)Dq3 + (53/12)DDq2 - (55/24)DDDq1 + C (251/720)DDDDq0 + ... ) C This process can be repeated until, hopefully, x1 no longer C changes. C IN THE EVENT THAT THE EQUATION y'(x) = f(x, y(x)) IS NOT C +Infinity, BUT SIMPLY CANNOT BE EVALUATED AT x=a, THE ONLY C DIFFERENCE IS THAT THERE IS NO NEED TO INTERCHANGE INDEPENDENT C AND DEPENDENT VARIABLES. C C ********** C .. Scalars in Common .. REAL EIG,EPSMIN INTEGER IND,T21 LOGICAL PR C .. C .. Local Scalars .. REAL CHM,CHNG,D2F,D2G,D3F,D3G,D4F,D4G,EPST,HT,HU, + OLDSS2,OLDYY2,ONE,RR,SLO,SOUT,SUMM,SUP,T,TEN5, + TIN,U,UOUT,USTAR,YLO,YOUT,YUP INTEGER I,K,KFLAG,NN3 C .. C .. Local Arrays .. REAL DF(4),DG(4),FF(6),GG(5),S(3),SS(6,3),UU(6) C .. C .. External Subroutines .. EXTERNAL GERK,PQEXT C .. C .. Intrinsic Functions .. INTRINSIC ABS,MAX,SIGN C .. C .. Scalar Arguments .. REAL DTHDE,THEND,TOUT,TSTAR,YSTAR INTEGER IFLAG,NC C .. C .. Array Arguments .. REAL ERR(3),TT(7),WORK(27),Y(3),YY(7,3) INTEGER IWORK(5) C .. C .. Subroutine Arguments .. EXTERNAL FG C .. C .. Common blocks .. COMMON /DATAF/EIG,IND COMMON /PRIN/PR,T21 COMMON /RNDOFF/EPSMIN C .. IFLAG = 1 ONE = 1.0D0 TEN5 = 100000.0D0 EPST = 1.D-7 C C INTEGRATE Y' = F(T,Y,YP), STARTING AT T = TSTAR. C TIN = SIGN(ONE,TSTAR) TT(1) = TIN YY(1,1) = THEND YY(1,2) = DTHDE YY(1,3) = 0.0D0 C THE NEXT THREE VALUES FOR YY(2,*) ARE JUST INITIAL GUESSES. YY(2,1) = YSTAR YY(2,2) = DTHDE YY(2,3) = 0.0D0 C YLO = -TEN5 YUP = TEN5 CHM = MAX(0.01D0*ONE,0.1D0*ABS(THEND)) C C NORMALLY IND = 1 OR 3; IND IS SET TO 2 OR 4 WHEN Y IS TO BE USED C AS THE INDEPENDENT VARIABLE IN FG(T,Y,YP), INSTEAD OF THE USUAL T. C BEFORE LEAVING THIS SUBROUTINE, IND IS RESET TO 1. C 10 CONTINUE T = TSTAR HT = TSTAR - TIN TT(2) = TSTAR Y(1) = YY(2,1) Y(2) = YY(2,2) Y(3) = YY(2,3) KFLAG = 1 IND = 1 DO 30 K = 3,6 TOUT = T + HT NN3 = 0 20 CONTINUE CALL GERK(FG,3,Y,T,TOUT,EPST,EPST,KFLAG,ERR,WORK,IWORK) IF (KFLAG.GT.3 .AND. ABS(TSTAR).GE.0.89D0) THEN TSTAR = TIN + 5.0D0* (TSTAR-TIN) IF (PR) WRITE (T21,FMT=*) ' KFLAG,NEW TSTAR = ',KFLAG, + TSTAR GO TO 10 END IF IF (KFLAG.EQ.3) THEN NN3 = NN3 + 1 IF (NN3.GE.6) THEN IFLAG = 0 RETURN END IF GO TO 20 END IF YOUT = Y(1) TT(K) = T YY(K,1) = YOUT YY(K,2) = Y(2) YY(K,3) = Y(3) 30 CONTINUE IND = 3 DO 40 I = 2,6 CALL FG(TT(I),YY(I,1),FF(I)) 40 CONTINUE C SET FF(1) = 0., JUST TO HAVE IT DEFINED BEFORE GOING INTO PQEXT() : C SUBROUTINE PQEXT, IF USED, TRIES TO DETERMINE WHETHER THE VALUES C OF FF(I) EXTRAPOLATE TO A FINITE OR INFINITE VALUE AT FF(1). FF(1) = 0.0D0 IF (NC.LE.4 .OR. NC.EQ.7) THEN CALL PQEXT(FF(1)) ELSE CALL FG(TT(1),YY(1,1),FF(1)) END IF WRITE (T21,FMT=*) ' in wr, ff1 = ',FF(1) IF (ABS(FF(1)).LE.50.0D0) THEN C C NOW WE WANT TO APPLY THE ABOVE CONSISTENCY CRITERION, C TO SEE IF THESE RESULTS ARE CONSISTENT WITH HAVING C INTEGRATED FROM (TT(1),YY(1,1). C WRITE (T21,FMT=*) ' in wr, finite ' DO 50 I = 1,4 DF(I) = FF(I+1) - FF(I) 50 CONTINUE D2F = DF(4) - DF(3) D3F = DF(4) - 2.0D0*DF(3) + DF(2) D4F = DF(4) - 3.0D0*DF(3) + 3.0D0*DF(2) - DF(1) SUMM = HT* (FF(5)-3.5D0*DF(4)+53.0D0*D2F/12.0D0- + 55.0D0*D3F/24.0D0+251.0D0*D4F/720.0D0) C C PRESUMABLY, YY(2,1) SHOULD BE YY(1,1) + SUMM. C OLDYY2 = YY(2,1) C TO COUNTER THE TENDENCY TO OVERCORRECT, USE A FACTOR C OF RR < 1.0 : RR = 0.95D0 YY(2,1) = YY(1,1) + RR*SUMM C C ALSO IMPROVE THE VALUE OF Y(2) AT TSTAR. C YY(2,2) = 0.5D0* (YY(1,2)+YY(3,2)) YY(2,3) = 0.5D0* (YY(1,3)+YY(3,3)) CHNG = YY(2,1) - OLDYY2 IF (CHNG.GE.0.0D0 .AND. OLDYY2.GT.YLO) YLO = OLDYY2 IF (CHNG.LE.0.0D0 .AND. OLDYY2.LT.YUP) YUP = OLDYY2 IF ((YY(2,1).GE.YUP.AND.YLO.GT.-TEN5) .OR. + (YY(2,1).LE.YLO.AND.YUP.LT.TEN5)) YY(2,1) = 0.5D0* + (YLO+YUP) CHNG = YY(2,1) - OLDYY2 IF (ABS(CHNG).GT.CHM) CHNG = SIGN(CHM,CHNG) YY(2,1) = OLDYY2 + CHNG c IF(PR)WRITE(T21,*) ' YY2, CHNG = ',YY(2,1),CHNG IF (ABS(YY(2,1)-OLDYY2).GT.EPST) GO TO 10 TOUT = TT(6) ELSE C C HERE, Y' IS ASSUMED INFINITE AT T = TIN. IN THIS CASE, C IT CANNOT BE EXPECTED TO APPROXIMATE Y WITH A POLYNOMIAL, C SO THE INDEPENDENT AND DEPENDENT VARIABLES ARE INTERCHANGED. C THE POINTS ARE ASSUMED EQUALLY SPACED. C WRITE (T21,FMT=*) ' in wr, infinite ' HU = (YY(6,1)-YY(1,1))/5.0D0 UU(1) = YY(1,1) SS(1,1) = TT(1) SS(1,2) = DTHDE SS(1,3) = 0.0D0 UU(2) = UU(1) + HU SS(2,1) = TSTAR SS(2,2) = DTHDE SS(2,3) = 0.0D0 USTAR = UU(2) SLO = -TEN5 SUP = TEN5 60 CONTINUE U = USTAR S(1) = SS(2,1) S(2) = SS(2,2) S(3) = SS(2,3) KFLAG = 1 IND = 2 DO 80 K = 3,6 UOUT = U + HU NN3 = 0 70 CONTINUE CALL GERK(FG,3,S,U,UOUT,EPST,EPST,KFLAG,ERR,WORK,IWORK) IF (KFLAG.GT.3 .AND. ABS(TSTAR).GE.0.89D0) THEN TSTAR = TIN + 5.0D0* (TSTAR-TIN) GO TO 10 END IF IF (KFLAG.EQ.3) THEN NN3 = NN3 + 1 IF (NN3.GE.6) THEN IFLAG = 0 RETURN END IF GO TO 70 END IF SOUT = S(1) UU(K) = U SS(K,1) = SOUT SS(K,2) = S(2) SS(K,3) = S(3) 80 CONTINUE IND = 4 DO 90 I = 2,5 CALL FG(UU(I),SS(I,1),GG(I)) 90 CONTINUE GG(1) = 0.0D0 DO 100 I = 1,4 DG(I) = GG(I+1) - GG(I) 100 CONTINUE D2G = DG(4) - DG(3) D3G = DG(4) - 2.0D0*DG(3) + DG(2) D4G = DG(4) - 3.0D0*DG(3) + 3.0D0*DG(2) - DG(1) SUMM = HU* (GG(5)-3.5D0*DG(4)+53.0D0*D2G/12.0D0- + 55.0D0*D3G/24.0D0+251.0D0*D4G/720.0D0) C C PRESUMABLY, SS(2,1) SHOULD BE SS(1,1) + SUMM. C OLDSS2 = SS(2,1) SS(2,1) = SS(1,1) + SUMM IF (SS(2,1).LE.-1.0D0) SS(2,1) = -1.0D0 + EPSMIN IF (SS(2,1).GE.1.0D0) SS(2,1) = 1.0D0 - EPSMIN C C ALSO IMPROVE THE VALUE OF Y(2) AT TSTAR. C SS(2,2) = 0.5D0* (SS(1,2)+SS(3,2)) SS(2,3) = 0.5D0* (SS(1,3)+SS(3,3)) CHNG = SS(2,1) - OLDSS2 C IF(PR)WRITE(T21,*) ' CHNG = ',CHNG IF (CHNG.GE.0.0D0 .AND. OLDSS2.GT.SLO) SLO = OLDSS2 IF (CHNG.LE.0.0D0 .AND. OLDSS2.LT.SUP) SUP = OLDSS2 IF ((SS(2,1).GE.SUP.AND.SLO.GT.-TEN5) .OR. + (SS(2,1).LE.SLO.AND.SUP.LT.TEN5)) SS(2,1) = 0.5D0* + (SLO+SUP) IF (ABS(SS(2,1)-OLDSS2).GT.EPST) GO TO 60 END IF IF (IND.EQ.4) THEN C NOW INTEGRATE AGAIN BUT WITH T AS THE INDEPENDENT VARIABLE: C THIS IS USEFUL FOR OBTAINING THE USUAL GLOBAL ERROR C ESTIMATES FOR THE QUANTITIES Y(1),Y(2),Y(3). 120 CONTINUE TT(6) = SS(6,1) YY(6,1) = UU(6) YY(6,2) = SS(6,2) YY(6,3) = SS(6,3) T = TT(6) HT = T - TIN Y(1) = YY(6,1) Y(2) = YY(6,2) Y(3) = YY(6,3) KFLAG = 1 IND = 1 DO 130 K = 7,12 TOUT = T + HT CALL GERK(FG,3,Y,T,TOUT,EPST,EPST,KFLAG,ERR,WORK,IWORK) IF (KFLAG.EQ.5) THEN EPST = 5.0D0*EPST WRITE (T21,FMT=*) ' in wr, epst = ',EPST GO TO 120 END IF 130 CONTINUE TOUT = T END IF C TOUT = TT(6) IND = 1 RETURN END C SUBROUTINE SETTHU(X,THU) C ********** C C THIS SUBROUTINE IS CALLED ONLY FROM SUBROUTINE LCO. C IT ESTABLISHES A DEFINITE VALUE FOR THU, C THE PHASE ANGLE FOR THE FUNCTION U, INCLUDING AN C APPROPRIATE INTEGER MULTIPLE OF PI C IT NEEDS THE NUMBERS MMW(*) FOUND IN THUM C C INPUT QUANTITIES: X,THU,YS,MMW,MMWD C OUTPUT QUANTITIES: THU C C ********** C .. Scalars in Common .. REAL HPI,PI,TWOPI INTEGER MMWD C .. C .. Arrays in Common .. REAL YS(200) INTEGER MMW(100) C .. C .. Local Scalars .. INTEGER I C .. C .. Common blocks .. COMMON /PASS/YS,MMW,MMWD COMMON /PIE/PI,TWOPI,HPI C .. C .. Scalar Arguments .. REAL THU,X C .. DO 10 I = 1,MMWD IF (X.GE.YS(MMW(I)) .AND. X.LE.YS(MMW(I)+1)) THEN IF (THU.GT.PI) THEN THU = THU + (I-1)*TWOPI RETURN