PROGRAM MAIN INTEGER I,IFLAG,INTAB,ISLFUN,JPAIRS,JSUM,NUMEIG,NWRITE DOUBLE PRECISION A,A1,A2,B,B1,B2,EIG,P0ATA,P0ATB,QFATA,QFATB,TOL DOUBLE PRECISION PAIRS(2),SLFUN(9) EXTERNAL PARAMS,SLEIGN,ZCOUNT C DATA NWRITE/6/ C CALL PARAMS(A,B,INTAB,P0ATA,QFATA,P0ATB,QFATB,A1,A2,B1,B2, 1 JPAIRS,PAIRS) C C COUNT NUMBER OF ZEROS IN (A,B) OF EIGENFUNCTION ASSOCIATED C WITH SMALLEST EIGENVALUE. C CALL ZCOUNT(A,B,A1,A2,B1,B2,JPAIRS,PAIRS,JSUM) WRITE (NWRITE,10) ' ZERO COUNT =',JSUM 10 FORMAT (A,I6) C C CALCULATE SMALLEST EIGENVALUE AND ASSOCIATED EIGENFUNCTION DATA. C NUMEIG = JSUM + 1 EIG = 0.0 TOL = 1.0D-6 ISLFUN = 0 CALL SLEIGN(A,B,INTAB,P0ATA,QFATA,P0ATB,QFATB,A1,A2,B1,B2, 1 NUMEIG,EIG,TOL,IFLAG,ISLFUN,SLFUN) WRITE (NWRITE,20) ' IFLAG,NUMEIG,EIG,TOL = ',IFLAG,NUMEIG,EIG,TOL 20 FORMAT (A,2I6,F12.5,1PE12.2) WRITE (NWRITE,30) ' SLFUN(1-9) =',(SLFUN(I),I=1,9) 30 FORMAT (A,F12.5/(13X,3F12.5)) STOP END C SUBROUTINE PARAMS(A,B,INTAB,P0ATA,QFATA,P0ATB,QFATB,A1,A2,B1,B2, 1 JPAIRS,PAIRS) INTEGER INTAB,JPAIRS DOUBLE PRECISION A,A1,A2,B,B1,B2,ONE,PI,P0ATA,P0ATB,QFATA,QFATB DOUBLE PRECISION PAIRS(*) INTRINSIC ATAN COMMON PI ONE = 1.0 PI = 4.0*ATAN(ONE) A = -PI B = PI INTAB = 1 P0ATA = -1.0 QFATA = 1.0 P0ATB = -1.0 QFATB = 1.0 A1 = 1.0 A2 = 0.0 B1 = 1.0 B2 = 0.0 JPAIRS = 1 PAIRS(1) = -PI/2.0 PAIRS(2) = PI/2.0 RETURN END C FUNCTION P(X) DOUBLE PRECISION P,X P = 1.0 RETURN END C FUNCTION Q(X) DOUBLE PRECISION Q,X Q = 4.0 RETURN END C FUNCTION R(X) DOUBLE PRECISION PI,R,X COMMON PI R = 0.0 IF (X.LE.-PI/2.0 .OR. X.GT.PI/2.0) R = 1.0 RETURN END ************************************************************************ SUBROUTINE ZCOUNT(A,B,A1,A2,B1,B2,JPAIRS,PAIRS,JSUM) INTEGER JPAIRS,JSUM DOUBLE PRECISION A,B,A1,A2,B1,B2 DOUBLE PRECISION PAIRS(2*JPAIRS) C ********** C C THIS SUBROUTINE COUNTS THE ZEROS, OVER SPECIFIED SUBINTERVALS, OF C THE SOLUTIONS OF A SECOND ORDER DIFFERENTIAL EQUATION IN THE FORM C C (P(X)*Y'(X))' + Q(X)*Y(X) = 0 ON (A,B) C C FOR USER-SUPPLIED COEFFICIENT FUNCTIONS P AND Q. THIS COUNT IN C TURN CORRESPONDS TO THE NUMBER OF ZEROS, IN THE INTERIOR OF (A,B), C OF THE FIRST EIGENFUNCTION OF THE RELATED STURM-LIOUVILLE PROBLEM C WHOSE WEIGHT FUNCTION VANISHES IDENTICALLY IN THE SUBINTERVALS. C C THE APPLICABLE INITIAL CONDITION DEPENDS UPON THREE CASES. C C CASE 1 -- ON A SUBINTERVAL WITH LEFT ENDPOINT A, C A1*Y(A) + A2*Y'(A)*P(A) = 0. C C CASE 2 -- ON A SUBINTERVAL WITH RIGHT ENDPOINT B, C B1*Y(B) + B2*Y'(B)*P(B) = 0. C C CASE 3 -- ON A SUBINTERVAL WITH NEITHER A NOR B AS ENDPOINT, C Y(XAA) = 0, WHERE XAA IS THE LEFT ENDPOINT. C C THE SUBROUTINE STATEMENT IS C C SUBROUTINE ZCOUNT(A,B,A1,A2,B1,B2,JPAIRS,PAIRS,JSUM) C C WHERE C C A AND B ARE INPUT VARIABLES DEFINING THE FULL INTERVAL. C A MUST BE LESS THAN B. C C A1 AND A2 ARE INPUT VARIABLES SET TO PRESCRIBE THE INITIAL C CONDITION AT A (CASE 1). C C B1 AND B2 ARE INPUT VARIABLES SET TO PRESCRIBE THE INITIAL C CONDITION AT B (CASE 2). C C JPAIRS IS AN INTEGER INPUT VARIABLE SET TO THE NUMBER OF C SPECIFIED SUBINTERVALS OF (A,B). C C PAIRS IS AN INPUT ARRAY OF LENGTH 2*JPAIRS WHOSE SUCCESSIVE C ORDERED ELEMENT PAIRS SPECIFY THE SUBINTERVALS. C C JSUM IS AN INTEGER OUTPUT VARIABLE SET TO THE TOTAL ZERO COUNT. C C SUBPROGRAMS CALLED C C SLEIGN-SUPPLIED ... G,GERK C C USER-SUPPLIED ..... P,Q C C THIS VERSION DATED 8/11/88. C BURTON S. GARBOW, ARGONNE NATIONAL LABORATORY C C ********** C .. SCALARS IN COMMON .. DOUBLE PRECISION Z C .. C .. LOCAL SCALARS .. INTEGER I,J,LFLAG,MF,ML DOUBLE PRECISION EPS,EPSMIN,H,ONE,PI,PSUM,PX,QX,RT,RTSAV, 1 T,U,V,W,WSAV,X,X50,XAA,XBB,XSAV,ZAV C .. C .. LOCAL ARRAYS .. INTEGER IWORK(5) DOUBLE PRECISION DS(98),GERROR(1),PS(99),QS(99),WORK(11),Y(1) C .. C .. EXTERNAL FUNCTIONS .. DOUBLE PRECISION P,Q EXTERNAL P,Q C .. C .. EXTERNAL SUBROUTINES .. EXTERNAL G,GERK C .. C .. INTRINSIC FUNCTIONS .. INTRINSIC ABS,ATAN,INT,SIGN,SQRT C .. C .. COMMON BLOCKS .. COMMON /ZEE/Z C .. C C SET CONSTANTS EPSMIN AND PI. (DATA ELEMENT U IN GERK C INTEGRATOR PACKAGE MUST ALSO BE SET TO THE VALUE OF EPSMIN.) C C ********** WARNING - THE FOLLOWING CONSTANT IS MACHINE DEPENDENT C (ACTIVATED VALUE IS FOR IEEE MACHINES): C EPSMIN - THE COMPUTER UNIT ROUNDOFF ERROR. C C CONSTANTS FOR SOME MACHINES (ACTIVATE BY REMOVING 'C' FROM COLUMN 1) C --------------------------------------------------------------------- C VAX 11/780 C DATA EPSMIN /1.4D-17/ C IEEE (ALLIANT FX/8, ENCORE MULTIMAX, SEQUENT BALANCE, SUN, ETC.) DATA EPSMIN /2.2D-16/ C IBM 3033 C DATA EPSMIN /2.2D-16/ C CRAY C DATA EPSMIN /2.5D-29/ C --------------------------------------------------------------------- C PI (VARIABLE ONE SET TO 1.0 EASES PRECISION CONVERSION.) ONE = 1.0 PI = 4.0*ATAN(ONE) C C SET RELATIVE AND ABSOLUTE ERROR TOLERANCES FOR GERK. C EPS = SQRT(EPSMIN) C JSUM = 0 DO 70 J = 1,JPAIRS XAA = PAIRS(2*J-1) XBB = PAIRS(2*J) C DO (CALCULATE-MODIFIED-PRUFER-TRANSFORMATION-CONSTANT) C C EVALUATE P,Q TO OBTAIN PRELIMINARY INFORMATION ABOUT THE C DIFFERENTIAL EQUATION. C C DO (SAMPLE-COEFFICIENTS) X50 = 0.5*((XBB+XAA)+(XBB-XAA)*EPSMIN) PX = P(X50) QX = Q(X50) PS(50) = PX QS(50) = QX/PX C C MF AND ML ARE THE LEAST AND GREATEST INDEX VALUES, RESPECTIVELY. C XSAV = X50 H = 0.9/40.0 DO 10 I=49,1,-1 IF (I.GE.10) T = H*(I-50) IF (I.LT.10) T = T - 0.75*(1.0+T) X = 0.5*((XBB+XAA)+(XBB-XAA)*T) PX = P(X) QX = Q(X) PS(I) = PX QS(I) = QX/PX DS(I) = XSAV - X XSAV = X MF = I 10 CONTINUE XSAV = X50 DO 20 I=51,99 IF (I.LE.90) T = H*(I-50) IF (I.GT.90) T = T + 0.75*(1.0-T) X = 0.5*((XBB+XAA)+(XBB-XAA)*T) PX = P(X) QX = Q(X) PS(I) = PX QS(I) = QX/PX DS(I-1) = X - XSAV XSAV = X ML = I - 1 20 CONTINUE C END (SAMPLE-COEFFICIENTS) C DO (ESTIMATE-PHASE-ANGLE-CHANGE) C C U MEASURES THE TOTAL LENGTH OF SUBINTERVALS WHERE Q/P .GT. 0.0. C ZAV IS THE AVERAGE VALUE OF SQRT(Q*P) OVER THOSE SUBINTERVALS. C U = 0.0 ZAV = 0.0 WSAV = QS(MF) IF (WSAV.GT.0.0) THEN RTSAV = SQRT(WSAV) ELSE RTSAV = 0.0 END IF DO 30 I=MF+1,ML W = QS(I) IF (W.GT.0.0) THEN U = U + DS(I-1) RT = SQRT(W) ELSE RT = 0.0 END IF IF (W.EQ.0.0 .OR. WSAV.EQ.0.0 .OR. 1 W.EQ.SIGN(W,WSAV)) THEN V = RT + RTSAV ELSE V = (W*RT+WSAV*RTSAV)/ABS(W-WSAV) END IF WSAV = W RTSAV = RT PSUM = DS(I-1)*V IF (U.GT.0.0) ZAV = ZAV + PSUM*(PS(I)+PS(I-1)) 30 CONTINUE ZAV = 0.25*ZAV C END (ESTIMATE-PHASE-ANGLE-CHANGE) Z = 1.0 IF (U.GT.0.0) Z = ZAV/U C END (CALCULATE-MODIFIED-PRUFER-TRANSFORMATION-CONSTANT) LFLAG = 1 IF (XAA.EQ.A) THEN C C CASE 1 ---------- C Y(1) = PI/2.0 IF (A1.NE.0.0) Y(1) = ATAN(-Z*A2/A1) IF (Y(1).LT.0.0) Y(1) = Y(1) + PI 40 CONTINUE CALL GERK(G,1,Y,XAA,XBB,EPS,EPS,LFLAG,GERROR,WORK,IWORK) IF (LFLAG.EQ.3) GO TO 40 JSUM = JSUM + INT((Y(1)+ABS(EPS))/PI) ELSE IF (XBB.EQ.B) THEN C C CASE 2 ---------- C Y(1) = PI/2.0 IF (B1.NE.0.0) Y(1) = ATAN(-Z*B2/B1) IF (Y(1).GT.0.0) Y(1) = Y(1) - PI 50 CONTINUE CALL GERK(G,1,Y,XBB,XAA,EPS,EPS,LFLAG,GERROR,WORK,IWORK) IF (LFLAG.EQ.3) GO TO 50 JSUM = JSUM - INT((Y(1)-ABS(EPS))/PI) ELSE C C CASE 3 ---------- C Y(1) = 0.0 60 CONTINUE CALL GERK(G,1,Y,XAA,XBB,EPS,EPS,LFLAG,GERROR,WORK,IWORK) IF (LFLAG.EQ.3) GO TO 60 JSUM = JSUM + INT((Y(1)+ABS(EPS))/PI) END IF 70 CONTINUE RETURN END SUBROUTINE G(X,Y,YP) DOUBLE PRECISION X DOUBLE PRECISION Y(1),YP(1) C ********** C C THIS SUBROUTINE EVALUATES THE DERIVATIVE FUNCTION FOR USE WITH C INTEGRATOR GERK IN SOLVING A DIFFERENTIAL EQUATION IN THE FORM C C (P(X)*Y'(X))' + Q(X)*Y(X) = 0 ON (A,B) C C FOR USER-SUPPLIED COEFFICIENT FUNCTIONS P AND Q. C C SUBPROGRAMS CALLED C C USER-SUPPLIED ..... P,Q C C ********** C .. SCALARS IN COMMON .. DOUBLE PRECISION Z C .. C .. LOCAL SCALARS .. DOUBLE PRECISION C,C2,QX,S,S2,XP C .. C .. EXTERNAL FUNCTIONS .. DOUBLE PRECISION P,Q EXTERNAL P,Q C .. C .. INTRINSIC FUNCTIONS .. INTRINSIC COS,SIN C .. C .. COMMON BLOCKS .. COMMON /ZEE/Z C .. QX = Q(X)/Z XP = Z/P(X) S = SIN(Y(1)) C = COS(Y(1)) S2 = S*S C2 = C*C YP(1) = XP*C2+QX*S2 RETURN END SUBROUTINE SLEIGN(A,B,INTAB,P0ATA,QFATA,P0ATB,QFATB,A1,A2,B1,B2, 1 NUMEIG,EIG,TOL,IFLAG,ISLFUN,SLFUN) INTEGER INTAB,NUMEIG,IFLAG,ISLFUN DOUBLE PRECISION A,B,P0ATA,QFATA,P0ATB,QFATB,A1,A2,B1,B2,EIG,TOL DOUBLE PRECISION SLFUN(9) C ********** C C THIS SUBROUTINE IS DESIGNED FOR THE CALCULATION OF A SPECIFIED C EIGENVALUE, EIG, OF A STURM-LIOUVILLE PROBLEM IN THE FORM C C (P(X)*Y'(X))' + (Q(X) + EIG*R(X))*Y(X) = 0 ON (A,B) C C FOR USER-SUPPLIED COEFFICIENT FUNCTIONS P, Q, AND R. C THE PROBLEM MAY BE EITHER NONSINGULAR OR SINGULAR. IN THE C NONSINGULAR CASE, BOUNDARY CONDITIONS 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 ARE PRESCRIBED BY SPECIFYING THE NUMBERS A1, A2, B1, AND B2. 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 THE SUBROUTINE STATEMENT IS C C SUBROUTINE SLEIGN(A,B,INTAB,P0ATA,QFATA,P0ATB,QFATB,A1,A2,B1,B2, C NUMEIG,EIG,TOL,IFLAG,ISLFUN,SLFUN) 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 CAN BE ARBITRARY. C C P0ATA, QFATA, P0ATB, AND QFATB ARE INPUT VARIABLES SET TO C 1.0 OR -1.0 AS THE FOLLOWING PROPERTIES OF P AND Q AT THE C INTERVAL ENDPOINTS ARE TRUE OR FALSE, RESPECTIVELY. C C P0ATA - P(A) IS ZERO. (IF TRUE, A IS SINGULAR.) C QFATA - Q(A) IS FINITE. (IF FALSE, A IS SINGULAR.) C P0ATB - P(B) IS ZERO. (IF TRUE, B IS SINGULAR.) C QFATB - Q(B) IS FINITE. (IF FALSE, B IS SINGULAR.) C C A1 AND A2 ARE INPUT VARIABLES SET TO PRESCRIBE THE BOUNDARY C CONDITION AT A. THEY CAN BE ARBITRARY IF A IS SINGULAR. C C B1 AND B2 ARE INPUT VARIABLES SET TO PRESCRIBE THE BOUNDARY C CONDITION AT B. THEY CAN BE ARBITRARY IF B IS SINGULAR. C C NUMEIG IS AN INTEGER VARIABLE. ON INPUT, IT SHOULD BE SET TO C THE INDEX OF THE DESIRED EIGENVALUE (INCREASING SEQUENCE). C ON OUTPUT, IT IS UNCHANGED UNLESS THE PROBLEM (APPARENTLY) C LACKS NUMEIG EIGENVALUES, IN WHICH CASE IT IS RESET TO THE C INDEX OF THE LARGEST EIGENVALUE. 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, SLEIGN 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 = 1 - SUCCESSFUL PROBLEM SOLUTION. C IFLAG = 2 - IMPROPER INPUT PARAMETERS. C IFLAG = 3 - NUMEIG EXCEEDS ACTUAL NUMBER OF EIGENVALUES. C IFLAG = 4 - SOME UNCERTAINTY ABOUT ACCURACY ESTIMATE TOL. C IFLAG = 5 - CONVERGENCE TOO SLOW, BEST RESULTS RETURNED. C IFLAG = 6 - FAILURE, INTEGRATOR COULD NOT COMPLETE. 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 ZERO OR NEGATIVE. 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 R) 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 SUBPROGRAMS CALLED C C USER-SUPPLIED ..... P,Q,R C C SLEIGN-SUPPLIED ... ALFBET,DXDT,ESTEIG,ESTPAC,F,GERK C C THIS VERSION DATED 8/11/88. C PAUL B. BAILEY, SANDIA LABORATORIES, ALBUQUERQUE C BURTON S. GARBOW, ARGONNE NATIONAL LABORATORY C C ********** C .. SCALARS IN COMMON .. INTEGER INTSAV LOGICAL EIGF DOUBLE PRECISION ASAV,BSAV,C1,C2,EIGSAV,Z C .. C .. LOCAL SCALARS .. INTEGER I,IA,IB,IMAX,IMIN,IOUT,IP,J,JFLAG,KFLAG,LFLAG,MF,ML, 1 NEIG,NMID LOGICAL AOK,BOK,BRACKT,CHNGAB,CHNGM,CONVRG,FYNYT,FYNYT1, 1 LIMIT,LOGIC,NEWTON,ONEDIG,PRIN,THEGT0,THELT0 DOUBLE PRECISION AA,AAA,ALFA,ASL,ASR,BB,BBB,BETA, 1 C,CHNG,CL,CR,DAV,DE,DEDW,DEN,DERIVL,DERIVR,DIST, 2 DPSIL,DPSIPL,DPSIPR,DPSIR,DT,DTHDA,DTHDAA,DTHDB, 3 DTHDBB,DTHDE,DTHDEA,DTHDEB,DTHETA,DTHOLD,E,EEE, 4 EIGLO,EIGLT,EIGRT,EIGUP,EL,ELIM,EMAX,EMIN,EOLD,EPS,EPSMIN, 5 ER1,ER2,ESTERR,FLO,FMAX,FUP,GMAX,GUESS,H,ONE,PI,PIN, 6 PSIL,PSIPL,PSIPR,PSIR,PX,QAV,QX,RATIO,RAV,RAY,RX, 7 SL,SQL,SQR,SR,T,T1,T2,T3,TAU,TEMP,THRESH,TMAX,TMID,TMP, 8 U,UL,UR,V,WL,X,X50,XAA,XBB,XBC,XMID,XSAV,ZAV C .. C .. LOCAL ARRAYS .. INTEGER IWORK(5) DOUBLE PRECISION DS(98),ERL(3),ERR(3),PS(99),QS(99),RS(99), 1 WORK(27),YL(3),YR(3) C .. C .. EXTERNAL FUNCTIONS .. DOUBLE PRECISION P,Q,R EXTERNAL P,Q,R C .. C .. EXTERNAL SUBROUTINES .. EXTERNAL ALFBET,DXDT,ESTEIG,ESTPAC,F,GERK C .. C .. INTRINSIC FUNCTIONS .. INTRINSIC ABS,ATAN,COS,EXP,INT,LOG,MAX,MIN,SIGN,SIN,TAN C .. C .. COMMON BLOCKS .. COMMON /DATADT/ASAV,BSAV,C1,C2,INTSAV COMMON /DATAF/EIGSAV,EIGF COMMON /ZEE/Z C .. C C SET CONSTANTS EPSMIN AND PI. (DATA ELEMENT U IN GERK C INTEGRATOR PACKAGE MUST ALSO BE SET TO THE VALUE OF EPSMIN.) C C ********** WARNING - THE FOLLOWING CONSTANT IS MACHINE DEPENDENT C (ACTIVATED VALUE IS FOR IEEE MACHINES): C EPSMIN - THE COMPUTER UNIT ROUNDOFF ERROR. C C CONSTANTS FOR SOME MACHINES (ACTIVATE BY REMOVING 'C' FROM COLUMN 1) C --------------------------------------------------------------------- C VAX 11/780 C DATA EPSMIN /1.4D-17/ C IEEE (ALLIANT FX/8, ENCORE MULTIMAX, SEQUENT BALANCE, SUN, ETC.) DATA EPSMIN /2.2D-16/ C IBM 3033 C DATA EPSMIN /2.2D-16/ C CRAY C DATA EPSMIN /2.5D-29/ C --------------------------------------------------------------------- C PI (VARIABLE ONE SET TO 1.0 EASES PRECISION CONVERSION.) ONE = 1.0 PI = 4.0*ATAN(ONE) C C SET OUTPUT DEVICE NUMBER. C IOUT = 6 C C CHECK INPUT PARAMETERS FOR ERRORS. IF ERRORS, RETURN IFLAG=2. C LOGIC = TOL.NE.0.0 .AND. 1.LE.INTAB .AND. INTAB.LE.4 .AND. 1 P0ATA*QFATA*P0ATB*QFATB*NUMEIG.NE.0.0 IF (INTAB.EQ.1) LOGIC = LOGIC .AND. A.LT.B IF (.NOT.LOGIC) THEN IFLAG = 2 RETURN END IF C C SET PRIN = .TRUE. TO TRIGGER TRACE PRINTOUT OF SUCCESSIVE STEPS. C PRIN = .FALSE. IF (TOL.LT.0.0) PRIN = .TRUE. C C SET EPS TO THE (INITIAL) INTEGRATION ACCURACY. C EPS = 0.001 C C AOK (BOK) SIGNALS, IF TRUE, THAT ENDPOINT A (B) IS NOT SINGULAR. C AOK = INTAB.LT.3 .AND. P0ATA.LT.0.0 .AND. QFATA.GT.0.0 BOK = (INTAB.EQ.1 .OR. INTAB.EQ.3) .AND. 1 P0ATB.LT.0.0 .AND. QFATB.GT.0.0 NEIG = ABS(NUMEIG) - 1 C C INITIAL C1 AND C2, USED IN THE MAPPING BETWEEN X AND T INTERVALS. C C1 = 1.0 C2 = 0.0 C DO (SAVE-INPUT-DATA) ASAV = A BSAV = B INTSAV = INTAB TAU = ABS(TOL) C END (SAVE-INPUT-DATA) C C EVALUATE P,Q,R TO OBTAIN PRELIMINARY INFORMATION ABOUT THE C DIFFERENTIAL EQUATION. C C DO (SAMPLE-COEFFICIENTS) THRESH = 1.0E+17 10 CONTINUE CALL DXDT(EPSMIN,TEMP,X50) PX = P(X50) QX = Q(X50) RX = R(X50) PS(50) = PX QS(50) = QX/PX RS(50) = RX/PX C C DAV,QAV,RAV ARE USED IN SPECIAL CASE ESTIMATION WHEN NUMEIG = 1,2. C EMIN = MIN(-Q(X)/R(X)), ACHIEVED AT INDEX VALUE IMIN. C EMAX = MAX(-Q(X)/R(X)), ACHIEVED AT INDEX VALUE IMAX. C MF AND ML ARE THE LEAST AND GREATEST INDEX VALUES, RESPECTIVELY. C DAV = 0.0 QAV = 0.0 RAV = 0.0 XSAV = X50 EMIN = THRESH EMAX = -THRESH IF (RX.NE.0.0) THEN EMIN = -QX/RX EMAX = EMIN IMIN = 50 IMAX = 50 END IF H = 0.9/40.0 DO 20 I=49,1,-1 IF (I.GE.10) T = H*(I-50) IF (I.LT.10) T = T - 0.75*(1.0+T) CALL DXDT(T,TEMP,X) PX = P(X) QX = Q(X) RX = R(X) PS(I) = PX QS(I) = QX/PX RS(I) = RX/PX DS(I) = XSAV - X DAV = DAV + DS(I) QAV = QAV + DS(I)*(0.5*(QS(I)+QS(I+1))-QAV)/DAV RAV = RAV + DS(I)*(0.5*(RS(I)+RS(I+1))-RAV)/DAV XSAV = X C C TRY TO AVOID OVERFLOW BY STOPPING WHEN FUNCTIONS ARE LARGE NEAR A. C FYNYT = (ABS(RX)+ABS(QX)+1.0/PX).LE.THRESH IF (RX.NE.0.0) THEN IF (-QX/RX.LT.EMIN) THEN EMIN = -QX/RX IMIN = I END IF IF (-QX/RX.GT.EMAX) THEN EMAX = -QX/RX IMAX = I END IF END IF MF = I IF (.NOT.FYNYT) GO TO 30 20 CONTINUE 30 CONTINUE AAA = T IF (AOK) AAA = -1.0 XSAV = X50 DO 40 I=51,99 IF (I.LE.90) T = H*(I-50) IF (I.GT.90) T = T + 0.75*(1.0-T) CALL DXDT(T,TEMP,X) PX = P(X) QX = Q(X) RX = R(X) PS(I) = PX QS(I) = QX/PX RS(I) = RX/PX DS(I-1) = X - XSAV DAV = DAV + DS(I-1) QAV = QAV + DS(I-1)*(0.5*(QS(I-1)+QS(I))-QAV)/DAV RAV = RAV + DS(I-1)*(0.5*(RS(I-1)+RS(I))-RAV)/DAV XSAV = X C C TRY TO AVOID OVERFLOW BY STOPPING WHEN FUNCTIONS ARE LARGE NEAR B. C FYNYT1 = (ABS(QX)+ABS(RX)+1.0/PX).LE.THRESH IF (RX.NE.0.0) THEN IF (-QX/RX.LT.EMIN) THEN EMIN = -QX/RX IMIN = I END IF IF (-QX/RX.GT.EMAX) THEN EMAX = -QX/RX IMAX = I END IF END IF ML = I - 1 IF (.NOT.FYNYT1) GO TO 50 40 CONTINUE 50 CONTINUE BBB = T IF (BOK) BBB = 1.0 LOGIC = C1.EQ.1.0 .AND. (.NOT.FYNYT .OR. .NOT.FYNYT1) C C MODIFY (T,X) TRANSFORMATION CORRESPONDING TO TRUNCATED INTERVAL. C IF (LOGIC) THEN C1 = 0.5*(BBB-AAA) C2 = 0.5*(AAA+BBB) GO TO 10 END IF C C ESTIMATE UPPER BOUND ELIM FOR EIG SUCH THAT BOUNDARY CONDITIONS C CAN BE SATISFIED. C ELIM = EMAX + 1.0 IF (INTAB.EQ.3 .OR. (P0ATA.GT.0.0 .AND. QFATA.LT.0.0)) THEN IF (-QS(MF)/RS(MF).LE.ELIM) THEN ELIM = -QS(MF)/RS(MF) IMAX = MF END IF END IF IF (INTAB.EQ.2 .OR. (P0ATB.GT.0.0 .AND. QFATB.LT.0.0)) THEN IF (-QS(ML)/RS(ML).LE.ELIM) THEN ELIM = -QS(ML)/RS(ML) IMAX = ML END IF END IF IF (INTAB.EQ.4) THEN ELIM = MIN(ELIM,-QS(MF)/RS(MF),-QS(ML)/RS(ML)) IF (-QS(MF)/RS(MF).EQ.ELIM) IMAX = MF IF (-QS(ML)/RS(ML).EQ.ELIM) IMAX = ML END IF ELIM = ELIM - EPSMIN IF (ELIM.EQ.0.0) ELIM = -EPSMIN LIMIT = ELIM.LE.EMAX C END (SAMPLE-COEFFICIENTS) PIN = (NEIG+1)*PI IF (EIG.EQ.0.0) THEN C DO (ESTIMATE-EIG) CALL ESTEIG(MF,ML,LIMIT,ELIM,EMAX,EMIN,PIN,QS,RS,DS,PS, 1 IMAX,IMIN,EEE,EIG,IA,IB,EL,WL,DEDW) C END (ESTIMATE-EIG) END IF GUESS = EIG C DO (SET-INITIAL-INTERVAL-AND-MATCHPOINT) IF (GUESS.NE.0.0) THEN C C REDUCE OVERLY LARGE GUESS FOR EIG TO UPPER BOUND IF NECESSARY. C IF (LIMIT .AND. EIG.GT.ELIM) EIG = ELIM EEE = EIG C DO (ESTIMATE-PHASE-ANGLE-CHANGE) CALL ESTPAC(MF,ML,EEE,PIN,QS,RS,DS,PS, 1 IA,IB,IP,TEMP,U,TMP) C END (ESTIMATE-PHASE-ANGLE-CHANGE) END IF C C CHOOSE INITIAL INTERVAL AS LARGE AS POSSIBLE THAT AVOIDS OVERFLOW. C IA AND IB ARE BOUNDARY INDICES FOR NONNEGATIVITY OF EIG*R+Q. C AA = -1.0 IF (.NOT.AOK) THEN AA = H*(IA-50) IF (IA.LT.10) AA = -1.0 + 0.1/4.0**(10-IA) END IF BB = 1.0 IF (.NOT.BOK) THEN BB = H*(IB-50) IF (IB.GT.90) BB = 1.0 - 0.1/4.0**(IB-90) END IF AA = AA + 0.6*(AAA-AA) BB = BB + 0.6*(BBB-BB) C C DETERMINE BOUNDARY VALUES ALFA AND BETA FOR THETA AT A AND B. C Z = 1.0 CALL ALFBET(A,INTAB,AA,A1,A2,EEE,P0ATA,QFATA,AOK, 1 ALFA,KFLAG,DERIVL) CALL ALFBET(B,INTAB,BB,B1,B2,EEE,P0ATB,QFATB,BOK, 1 BETA,JFLAG,DERIVR) IF (.NOT.BOK) BETA = PI - BETA C C TAKE BOUNDARY CONDITIONS INTO ACCOUNT IN ESTIMATION OF EIG. C PIN = PIN + BETA - ALFA - PI IF (GUESS.EQ.0.0) EEE = EL + DEDW*(PIN-WL) C C SUBROUTINE ESTPAC MUST BE CALLED AGAIN BECAUSE PIN HAS CHANGED. C C DO (ESTIMATE-PHASE-ANGLE-CHANGE) CALL ESTPAC(MF,ML,EEE,PIN,QS,RS,DS,PS,IA,IB,IP,TEMP,U,ZAV) C END (ESTIMATE-PHASE-ANGLE-CHANGE) C C CHOOSE THE CONSTANT Z. C IF (U.GT.0.0) Z = ZAV/U C C RESET BOUNDARY VALUES ALFA AND BETA. C CALL ALFBET(A,INTAB,AA,A1,A2,EEE,P0ATA,QFATA,AOK, 1 ALFA,KFLAG,DERIVL) CALL ALFBET(B,INTAB,BB,B1,B2,EEE,P0ATB,QFATB,BOK, 1 BETA,JFLAG,DERIVR) IF (.NOT.BOK) BETA = PI - BETA IF (PRIN) WRITE(IOUT,'(A,E22.14,A,E22.14)') 1 ' ALFA=',ALFA,' BETA=',BETA C C SPECIAL CASE FORMULA FOR ESTIMATION OF EIG WHEN NUMEIG = 1,2. C IF (U.EQ.0.0 .AND. NEIG.LE.1 .AND. (BETA+NEIG*PI).LT.ALFA) THEN XBC = MAX(-1.0/TAN(ALFA),1.0/TAN(BETA)) EEE = -(XBC*XBC-QAV)/RAV DEDW = XBC*(1.0+XBC*XBC)/RAV END IF C C CHOOSE INITIAL MATCHING POINT TMID. C TMID = H*(IP-50) IF (TMID.LT.-0.8) TMID = -0.4 IF (TMID.GT.0.8) TMID = 0.4 IF (PRIN) WRITE(IOUT,'(A,E15.7,A,F11.8,A,E15.7)') 1 ' ESTIM=',EEE,' TMID=',TMID,' Z=',Z IF (PRIN) WRITE(IOUT,'(A,F11.8,A,F11.8,A,F11.8,A,F11.8)') 1 ' AAA=',AAA,' AA=',AA,' BB=',BB,' BBB=',BBB C END (SET-INITIAL-INTERVAL-AND-MATCHPOINT) 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 CHNGM - MATCHING POINT TMID SHOULD BE CHANGED. C CHNGAB - ONE OR BOTH ENDPOINTS SHOULD BE MOVED FARTHER OUT. 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 EIGF - DERIVATIVE ARGUMENT IS IN ORIGINAL COORDINATE SYSTEM. C EIG = EEE EIGF = .FALSE. CHNGM = .FALSE. CHNGAB = .TRUE. 60 CONTINUE IF (CHNGAB) THEN C DO (INITIAL-IZE) CHNGAB = .FALSE. BRACKT = .FALSE. CONVRG = .FALSE. THELT0 = .FALSE. THEGT0 = .FALSE. EIGLO = 0.0 EIGLT = 0.0 EIGRT = 0.0 EIGUP = 0.0 DTHOLD = 0.0 C END (INITIAL-IZE) END IF C C RECOMPUTE BOUNDARY CONDITIONS AT SINGULAR ENDPOINT(S). C C DO (RESET-BOUNDARY-CONDITIONS) DERIVL = 0.0 IF (.NOT.AOK) CALL ALFBET(A,INTAB,AA,A1,A2,EIG, 1 P0ATA,QFATA,.FALSE.,ALFA,KFLAG,DERIVL) DERIVR = 0.0 IF (.NOT.BOK) THEN CALL ALFBET(B,INTAB,BB,B1,B2,EIG,P0ATB,QFATB,.FALSE., 1 BETA,JFLAG,DERIVR) BETA = PI - BETA END IF IF (PRIN) WRITE(IOUT,'(A,E22.14,A,E22.14)') 1 ' ALFA=',ALFA,' BETA=',BETA C END (RESET-BOUNDARY-CONDITIONS) 70 CONTINUE IF (EIG.NE.GUESS .AND. .NOT.BRACKT) THEN C C IF INITIAL GUESS WAS SUPPLIED, CHECK THAT BOUNDARY CONDITIONS C CAN BE SATISFIED AT SINGULAR ENDPOINTS. IF NOT, TRY FOR C SLIGHTLY LOWER EIG CONSISTENT WITH BOUNDARY CONDITIONS. C 80 CONTINUE IF (.NOT.AOK) CALL ALFBET(A,INTAB,AA,A1,A2,EIG, 1 P0ATA,QFATA,.FALSE.,TMP,KFLAG,TEMP) IF (.NOT.BOK) CALL ALFBET(B,INTAB,BB,B1,B2,EIG, 1 P0ATB,QFATB,.FALSE.,TMP,JFLAG,TEMP) IF (KFLAG*JFLAG.NE.1) THEN IF (THEGT0) EIG = 0.6*EIG + 0.4*EIGUP IF (THELT0) EIG = 0.6*EIG + 0.4*EIGLO IF (THELT0 .AND. EIGLO.LT.ELIM) EIGUP = ELIM GO TO 80 END IF END IF C DO (OBTAIN-DTHETA-WITH-ONE-CORRECT-DIGIT) 90 CONTINUE IF (PRIN) WRITE(IOUT,'(/A,E22.14,A,E10.3,A,E10.3)') 1 ' GUESS=',EIG,' EPS=',EPS,' TMID=',TMID C DO (INTEGRATE-FOR-DTHETA) C DO (SET-INITIAL-CONDITIONS) DTHDEA = DERIVL DTHDAA = 0.0 IF (.NOT.AOK) THEN CALL DXDT(AA,DT,X) PX = P(X)/Z QX = Q(X)/Z RX = R(X)/Z C = EIG*RX + QX DTHDAA = -(COS(ALFA)**2/PX + 1 C*SIN(ALFA)**2)*DT C C TWO SPECIAL CASES FOR DTHDAA. C IF (C.GE.0.0 .AND. P0ATA.LT.0.0 .AND. 1 QFATA.LT.0.0) DTHDAA = DTHDAA + 1 ALFA*DT/(X-A) IF (C.GE.0.0 .AND. P0ATA.GT.0.0 .AND. 1 QFATA.GT.0.0) DTHDAA = DTHDAA + 2 (ALFA-0.5*PI)*DT/(X-A) END IF DTHDEB = -DERIVR DTHDBB = 0.0 IF (.NOT.BOK) THEN CALL DXDT(BB,DT,X) PX = P(X)/Z QX = Q(X)/Z RX = R(X)/Z C = EIG*RX + QX DTHDBB = -(COS(BETA)**2/PX + 1 C*SIN(BETA)**2)*DT C C TWO SPECIAL CASES FOR DTHDBB. C IF (C.GE.0.0 .AND. P0ATB.LT.0.0 .AND. 1 QFATB.LT.0.0) DTHDBB = DTHDBB + 2 (PI-BETA)*DT/(B-X) IF (C.GE.0.0 .AND. P0ATB.GT.0.0 .AND. 1 QFATB.GT.0.0) DTHDBB = DTHDBB + 2 (0.5*PI-BETA)*DT/(B-X) END IF TMAX = TMID GMAX = ABS(DTHDEA) EIGSAV = EIG C END (SET-INITIAL-CONDITIONS) C T C YL = (THETA,D(THETA)/D(EIG),D(THETA)/DA) C T = AA YL(1) = ALFA YL(2) = DTHDEA YL(3) = 0.0 C C USE INTEGRATOR IN ONE-STEP MODE TOWARDS CHANGE TO DIFFERENT TMID. C LFLAG = 1 IF (CHNGM) LFLAG = -1 100 CONTINUE CALL GERK(F,3,YL,T,TMID,EPS,EPS,LFLAG,ERL, 1 WORK,IWORK) IF (LFLAG.EQ.-2 .AND. T.GT.-0.8 .AND. 1 ABS(YL(2)).GT.GMAX) THEN TMAX = T GMAX = ABS(YL(2)) END IF IF (LFLAG.EQ.3 .OR. LFLAG.EQ.-2) GO TO 100 IF (LFLAG.GT.3) THEN IFLAG = 6 RETURN END IF DTHDA = DTHDAA*EXP(-2.0*YL(3)) C T C YR = (THETA,D(THETA)/D(EIG),D(THETA)/DB) C T = BB YR(1) = BETA + PI*NEIG YR(2) = DTHDEB YR(3) = 0.0 C C USE INTEGRATOR IN ONE-STEP MODE TOWARDS CHANGE TO DIFFERENT TMID. C LFLAG = 1 IF (CHNGM) LFLAG = -1 110 CONTINUE CALL GERK(F,3,YR,T,TMID,EPS,EPS,LFLAG,ERR, 1 WORK,IWORK) IF (LFLAG.EQ.-2 .AND. T.LT.0.8 .AND. 1 ABS(YR(2)).GT.GMAX) THEN TMAX = T GMAX = ABS(YR(2)) END IF IF (LFLAG.EQ.3 .OR. LFLAG.EQ.-2) GO TO 110 IF (LFLAG.GT.3) THEN IFLAG = 6 RETURN END IF DTHDB = DTHDBB*EXP(-2.0*YR(3)) C C DTHETA MEASURES THETA DIFFERENCE FROM LEFT AND RIGHT INTEGRATIONS. C DTHETA = YL(1) - YR(1) DTHDE = YL(2) - YR(2) ER1 = ERL(1) - ERR(1) ER2 = ERL(2) - ERR(2) TMID = TMAX C END (INTEGRATE-FOR-DTHETA) C C DEFINE ONEDIG TO TRY TO BE SURE OF ONE CORRECT DIGIT IN DTHETA. C REDO INTEGRATIONS WITH TIGHTER TOLERANCE UNTIL ONEDIG IS TRUE. C ONEDIG = ABS(ER1).LE.0.1*ABS(DTHETA) .AND. 1 ABS(ER2).LE.0.1*ABS(DTHDE) NEWTON = ABS(DTHETA).LT.0.06 IF (NEWTON) THEN C DO (COMPUTE-CONVRG) C C MEASURE CONVERGENCE AFTER ADDING SEPARATE CONTRIBUTIONS TO ERROR. C T1 = ABS(DTHETA)+50.0*ABS(ER1) T2 = (1.0+AA)*ABS(DTHDA) T3 = (1.0-BB)*ABS(DTHDB) ESTERR = (T1+T2+T3)/ABS(DTHDE)/MAX(ONE,ABS(EIG)) CONVRG = ESTERR.LE.TAU IF (PRIN) WRITE(IOUT,'(A,L2)') 1 ' CONVERGE=',CONVRG IF (PRIN .AND. .NOT.CONVRG) 1 WRITE(IOUT,'(A,E15.7)') 2 ' ESTIM. ACC.=',ESTERR C END (COMPUTE-CONVRG) END IF IF (.NOT.ONEDIG .OR. 1 ABS(ER1).GT.0.01*ABS(DTHETA)) THEN C C REDUCE LOCAL ERROR CRITERION, BUT RETURN IFLAG=5 IF TOO SMALL. C EPS = 0.05*EPS IF (EPS.LE.EPSMIN) THEN IFLAG = 5 RETURN END IF END IF IF (.NOT.(ONEDIG .OR. CONVRG)) GO TO 90 IF (ABS(DTHETA).LT.0.1) CHNGM = .FALSE. IF (PRIN) WRITE(IOUT,'(A,E15.7,A,E15.7)') 1 ' DTHETA=',DTHETA,' DTHDE=',DTHDE IF (PRIN) WRITE(IOUT,'(/A,E15.7,A,E15.7)') 1 ' THETAL=',YL(1),' THETAR=',YR(1) C END (OBTAIN-DTHETA-WITH-ONE-CORRECT-DIGIT) C DO (SET-BRACKET-AND-FORM-NEWTON-ITERATES) C C EIG IS BRACKETED WHEN BOTH THEGT0=.TRUE. AND THELT0=.TRUE. C IF (DTHETA*DTHDE.GT.0.0) 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 BRACKT = THELT0 .AND. THEGT0 IF (PRIN) WRITE(IOUT,'(A,E22.14,A,E22.14)') 1 ' EIGRT=',EIGRT,' EIGUP=',EIGUP IF (PRIN) WRITE(IOUT,'(A,E22.14,A,E22.14)') 1 ' EIGLT=',EIGLT,' EIGLO=',EIGLO C END (SET-BRACKET-AND-FORM-NEWTON-ITERATES) IF (NEWTON) THEN CHNGM = .TRUE. IF (CONVRG) THEN CHNG = DTHDA*(-1.0-AA) - DTHDB*(1.0-BB) TEMP = (DTHETA+CHNG)/DTHDE EIG = EIG - TEMP TOL = ABS(TEMP)/MAX(ONE,ABS(EIG)) ELSE CHNGAB = T1.LT.0.5*(T2+T3) C C MOVE ENDPOINT(S) OUT OR TAKE NEWTON STEP, ACCORDING TO CHNGAB. C IF (CHNGAB) THEN C DO (MOVE-ENDPOINTS) IF (T2.GT.T1 .AND. AA.GT.AAA) 1 AA = AA + 0.8*(-1.0-AA) IF (T3.GE.T1 .AND. BB.LT.BBB) 1 BB = BB + 0.8*(1.0-BB) AA = MAX(AA,AAA) BB = MIN(BB,BBB) IF ((AAA-AA).EQ.(BBB-BB)) THEN C C CANNOT MOVE ENDPOINT(S) AGAIN. STORE ESTIMATES AND RETURN IFLAG=5. C CHNG = (DTHDA-DTHDB)*(AAA-AA) TEMP = (DTHETA+CHNG)/DTHDE EIG = EIG - TEMP TOL = ABS(TEMP)/MAX(ONE,ABS(EIG)) IFLAG = 5 RETURN END IF EEE = EIG IF (PRIN) WRITE(IOUT,'(A,2E15.8)') 1 ' NEW ENDPOINTS ',AA,BB C END (MOVE-ENDPOINTS) ELSE EIG = EIG - DTHETA/DTHDE END IF END IF ELSE IF (BRACKT) THEN C C OBTAIN NEXT ESTIMATE OF EIG BY BISECTION OR LINEAR INTERPOLATION. C FMAX = MAX(-FLO,FUP) EOLD = EIG RATIO = DTHETA/DTHOLD EIG = 0.5*(EIGLO+EIGUP) IF (FMAX.LE.1.5) THEN U = -FLO/(FUP-FLO) DIST = EIGUP - EIGLO EIG = EIGLO + U*DIST V = MIN(EIGLT,EIGRT) IF (EIG.LE.V) EIG = 0.5*(EIG+V) V = MAX(EIGLT,EIGRT) IF (EIG.GE.V) EIG = 0.5*(EIG+V) DE = EIG - EOLD CHNGAB = RATIO.GE.0.4 .AND. .NOT.(AOK.AND.BOK) IF (ABS(DE).LT.EPSMIN) THEN TOL = ABS(DE)/MAX(ONE,ABS(EIG)) IFLAG = 5 RETURN END IF END IF CHNGM = .NOT.CHNGM .AND. RATIO.GE.0.4 ELSE C DO (TRY-FOR-BRACKET) C C TAKE TWICE THE NEWTON STEP IN TRYING FOR A BRACKET. C IF (EIG.EQ.EEE) THEN IF (GUESS.NE.0.0) DEDW = 1.0/DTHDE CHNG = -(DEDW+1.0/DTHDE)*DTHETA IF (ABS(CHNG).GT.0.1*ABS(EIG)) 1 CHNG = -0.1*SIGN(EIG,DTHETA) ELSE CHNG = -2.0*DTHETA/DTHDE END IF LOGIC = EIG.NE.EEE .AND. DTHOLD.LT.0.0 .AND. 1 LIMIT .AND. CHNG.GT.(ELIM-EIG) IF (LOGIC) THEN CHNG = 0.99*(ELIM-EIG+EPSMIN) IF (CHNG.LT.EPSMIN) THEN C C IF CHANGE IS TOO SMALL, EIG IS PRESUMED NOT TO EXIST (IFLAG=3). C NUMEIG = NEIG - INT(-DTHETA/PI) IFLAG = 3 RETURN END IF C C LIMIT CHANGE IN EIG TO A FACTOR OF 10. C ELSE IF (ABS(CHNG).GT.(1.0+10.0*ABS(EIG))) THEN CHNG = SIGN(1.0+10.0*ABS(EIG),CHNG) ELSE IF (ABS(EIG).GE.1.0 .AND. 1 ABS(CHNG).LT.0.1*ABS(EIG)) THEN CHNG = 0.1*SIGN(EIG,CHNG) END IF EOLD = EIG EIG = EIG + CHNG C END (TRY-FOR-BRACKET) END IF END IF DTHOLD = DTHETA IF (.NOT.(CONVRG .OR. CHNGAB .OR. NEWTON)) GO TO 70 IF (.NOT.CONVRG) GO TO 60 IFLAG = 1 IF (PRIN) WRITE(IOUT,'(A,I7,A,E22.14,A,E10.3)') 1 ' NUMEIG=',NUMEIG,' EIG=',EIG,' TOL=',TOL C DO (COMPUTE-EIGENFUNCTION-DATA) C C CONVERT FROM T TO X VALUES, FILL 7 OF FIRST 9 LOCATIONS OF SLFUN. C CALL DXDT(TMID,TEMP,XMID) CALL DXDT(AA,TEMP,XAA) CALL DXDT(BB,TEMP,XBB) SLFUN(1) = XMID SLFUN(2) = XAA SLFUN(3) = ALFA SLFUN(5) = XBB SLFUN(6) = BETA + PI*NEIG SLFUN(8) = EPS SLFUN(9) = Z C C COMPUTE SLFUN(4), SLFUN(7) TOWARDS NORMALIZING THE EIGENFUNCTION. C EIGSAV = EIG Z = -Z T = AA YL(1) = ALFA YL(2) = DTHDEA YL(3) = 0.0 LFLAG = 1 120 CONTINUE CALL GERK(F,3,YL,T,TMID,EPS,EPS,LFLAG,ERL,WORK,IWORK) IF (LFLAG.EQ.3) GO TO 120 T = BB YR(1) = BETA + PI*NEIG YR(2) = DTHDEB YR(3) = 0.0 LFLAG = 1 130 CONTINUE CALL GERK(F,3,YR,T,TMID,EPS,EPS,LFLAG,ERR,WORK,IWORK) IF (LFLAG.EQ.3) GO TO 130 Z = -Z SL = SIN(YL(1)) SR = SIN(YR(1)) CL = COS(YL(1)) CR = COS(YR(1)) UL = (YL(2)-DTHDEA*EXP(-2.0*YL(3)))*Z UR = (YR(2)-DTHDEB*EXP(-2.0*YR(3)))*Z ASL = ABS(SL) ASR = ABS(SR) DEN = 0.5*LOG(UL*ASR*ASR-UR*ASL*ASL) SLFUN(4) = LOG(ASR) - YL(3) - DEN SLFUN(7) = LOG(ASL) - YR(3) - DEN C END (COMPUTE-EIGENFUNCTION-DATA) C DO (CHECK-MATCHING-VALUES-OF-EIGENFUNCTION) C C PERFORM FINAL CHECK ON EIG. RETURN IFLAG=4 IF NOT ACCURATE ENOUGH. C E = ASR*EXP(-DEN) PSIL = E*SL PSIPL = E*CL SQL = E*E*UL DPSIL = PSIL*ERL(3) + PSIPL*ERL(1) DPSIPL = PSIL*ERL(1) + PSIPL*ERL(3) PSIPL = PSIPL*Z E = ASL*EXP(-DEN) PSIR = E*SR PSIPR = E*CR SQR = E*E*UR DPSIR = PSIR*ERR(3) + PSIPR*ERR(1) DPSIPR = PSIR*ERR(1) + PSIPR*ERR(3) PSIPR = PSIPR*Z RAY = EIG + (PSIL*PSIPL-PSIR*PSIPR)/(SQL-SQR) IF (PRIN) THEN WRITE(IOUT,'(A,E22.14)') ' RAY=',RAY WRITE(IOUT,'(A,E22.14,A,E22.14)') 1 ' PSIL=',PSIL,' PSIR=',PSIR WRITE(IOUT,'(A,E22.14,A,E22.14)') 1 ' PSIPL=',PSIPL,' PSIPR=',PSIPR WRITE(IOUT,'(A,E22.14,A,E22.14)') 1 ' SQL=',SQL,' SQR=',SQR WRITE(IOUT,'(A,E22.14,A,E22.14)') 1 ' DPSIL=',DPSIL,' DPSIR=',DPSIR WRITE(IOUT,'(A,E22.14,A,E22.14)') 1 ' DPSIPL=',DPSIPL,' DPSIPR=',DPSIPR END IF C END (CHECK-MATCHING-VALUES-OF-EIGENFUNCTION) IF (ABS(RAY-EIG).GT.TAU*MAX(ONE,ABS(EIG))) IFLAG = 4 IF (ISLFUN.GT.0) THEN C C CALCULATE SELECTED EIGENFUNCTION VALUES BY INTEGRATION. C C DO (GENERATE-EIGENFUNCTION-VALUES) EIGF = .TRUE. NMID = 0 DO 140 I=1,ISLFUN IF (SLFUN(9+I).LE.SLFUN(1)) NMID = I 140 CONTINUE IF (NMID.GT.0) THEN X = SLFUN(2) YL(1) = SLFUN(3) YL(2) = 0.0 YL(3) = SLFUN(4) LFLAG = 1 DO 160 J=1,NMID 150 CONTINUE CALL GERK(F,3,YL,X,SLFUN(J+9),SLFUN(8),SLFUN(8), 1 LFLAG,ERL,WORK,IWORK) IF (LFLAG.EQ.3) GO TO 150 IF (LFLAG.EQ.6) LFLAG = 2 SLFUN(J+9) = EXP(YL(3))*SIN(YL(1)) 160 CONTINUE END IF IF (NMID.LT.ISLFUN) THEN X = SLFUN(5) YR(1) = SLFUN(6) YR(2) = 0.0 YR(3) = SLFUN(7) LFLAG = 1 DO 180 J=ISLFUN,NMID+1,-1 170 CONTINUE CALL GERK(F,3,YR,X,SLFUN(J+9),SLFUN(8),SLFUN(8), 1 LFLAG,ERR,WORK,IWORK) IF (LFLAG.EQ.3) GO TO 170 IF (LFLAG.EQ.6) LFLAG = 2 SLFUN(J+9) = EXP(YR(3))*SIN(YR(1)) 180 CONTINUE END IF C END (GENERATE-EIGENFUNCTION-VALUES) END IF RETURN END SUBROUTINE ESTEIG(MF,ML,LIMIT,ELIM,EMAX,EMIN,PIN,QS,RS,DS,PS, 1 IMAX,IMIN,EEE,EIG,IA,IB,EL,WL,DEDW) INTEGER MF,ML,IMAX,IMIN,IA,IB LOGICAL LIMIT DOUBLE PRECISION ELIM,EMAX,EMIN,PIN,EEE,EIG,EL,WL,DEDW DOUBLE PRECISION QS(ML),RS(ML),DS(ML),PS(ML) C ********** C C THIS SUBROUTINE GENERATES AN INITIAL GUESS FOR A SPECIFIED C EIGENVALUE OF A STURM-LIOUVILLE PROBLEM IN THE FORM C C (P(X)*Y'(X))' + (Q(X) + EIG*R(X))*Y(X) = 0 ON (A,B) C C FOR USER-SUPPLIED COEFFICIENT FUNCTIONS P, Q, AND R. IT IS C CALLED FROM SLEIGN WHEN NO INITIAL GUESS IS PROVIDED BY THE USER. C C THE METHOD USED IS TO APPROXIMATELY SOLVE THE EQUATION FOR EIG C C INTEGRAL (SQRT((EIG*R+Q)/P)) = NUMEIG*PI C C WHERE THE INTEGRAL IS TAKEN OVER THOSE X IN (A,B) FOR WHICH C C (EIG*R+Q)/P .GT. 0 C C AND NUMEIG IS THE INDEX OF THE DESIRED EIGENVALUE (PIN=NUMEIG*PI). C C SUBPROGRAMS CALLED C C SLEIGN-SUPPLIED ... ESTPAC C C ********** C .. SCALARS IN COMMON .. INTEGER INTAB DOUBLE PRECISION A,B,C1,C2 C .. C .. LOCAL SCALARS .. INTEGER IE,IP LOGICAL LOGIC DOUBLE PRECISION BALLPK,EU,FNEW,FOLD,SUM,TEMP,U,WU C .. C .. EXTERNAL SUBROUTINES .. EXTERNAL ESTPAC C .. C .. INTRINSIC FUNCTIONS .. INTRINSIC ABS,MIN C .. C .. COMMON BLOCKS .. COMMON /DATADT/A,B,C1,C2,INTAB C .. EEE = MIN(ELIM,EMAX) C DO (ESTIMATE-PHASE-ANGLE-CHANGE) CALL ESTPAC(MF,ML,EEE,PIN,QS,RS,DS,PS,IA,IB,IP,SUM,U,TEMP) C END (ESTIMATE-PHASE-ANGLE-CHANGE) C C CHOOSE BOUNDS FOR EIG AND ASSOCIATE FUNCTION (INTEGRAL) VALUES. C EL = EMIN WL = 0.0 EU = EEE WU = SUM IF (LIMIT .AND. WU.LT.PIN) THEN EIG = ELIM ELSE IF (U.EQ.0.0) THEN EL = EMAX EEE = EMAX + 1.0 C DO (ESTIMATE-PHASE-ANGLE-CHANGE) CALL ESTPAC(MF,ML,EEE,PIN,QS,RS,DS,PS, 1 IA,IB,IP,SUM,U,TEMP) C END (ESTIMATE-PHASE-ANGLE-CHANGE) EU = EEE WU = SUM END IF 10 CONTINUE IF (WU.LE.PIN) THEN C C INCREASE TRIAL VALUE IF INTEGRAL IS STILL TOO SMALL. C EL = EU WL = WU EEE = EU + ((PIN-WU+3.0)/U)**2 C DO (ESTIMATE-PHASE-ANGLE-CHANGE) CALL ESTPAC(MF,ML,EEE,PIN,QS,RS,DS,PS, 1 IA,IB,IP,SUM,U,TEMP) C END (ESTIMATE-PHASE-ANGLE-CHANGE) EU = EEE WU = SUM GO TO 10 END IF C C EIG IS BRACKETED. NOW TRY TO REDUCE THE SIZE OF THE BRACKET C BY SEARCHING AMONG THE SAVED VALUES OF -QS()/RS(). C 20 CONTINUE IF (ABS(IMAX-IMIN).GE.2 .AND. EU.LE.EMAX) THEN IE = (IMAX+IMIN)/2 IF (RS(IE).NE.0.0) THEN EEE = -QS(IE)/RS(IE) C DO (ESTIMATE-PHASE-ANGLE-CHANGE) CALL ESTPAC(MF,ML,EEE,PIN,QS,RS,DS,PS, 1 IA,IB,IP,SUM,U,TEMP) C END (ESTIMATE-PHASE-ANGLE-CHANGE) IF (SUM.GT.PIN) THEN IMAX = IE WU = SUM EU = EEE ELSE IMIN = IE WL = SUM EL = EEE END IF GO TO 20 END IF END IF 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.0 BALLPK = (PIN/(A-B))**2 LOGIC = .TRUE. 30 CONTINUE IF (LOGIC) THEN LOGIC = (WL.LT.(PIN-1.0) .OR. WU.GT.(PIN+1.0)) EEE = EL + DEDW*(PIN-WL) FNEW = MIN(PIN-WL,WU-PIN) IF (FNEW.GT.0.4*FOLD .OR. FNEW.LE.1.0) 1 EEE = 0.5*(EL+EU) IF (INTAB.EQ.1 .AND. ABS(EEE).GT.1.0E+3*BALLPK) THEN EIG = BALLPK RETURN ELSE IF (INTAB.NE.1 .AND. ABS(EEE).GT.1.0E+6) THEN EIG = 1.0 RETURN ELSE FOLD = FNEW C DO (ESTIMATE-PHASE-ANGLE-CHANGE) CALL ESTPAC(MF,ML,EEE,PIN,QS,RS,DS,PS, 1 IA,IB,IP,SUM,U,TEMP) C END (ESTIMATE-PHASE-ANGLE-CHANGE) IF (SUM.LT.PIN) THEN EL = EEE WL = SUM ELSE EU = EEE WU = SUM END IF DEDW = (EU-EL)/(WU-WL) GO TO 30 END IF END IF END IF RETURN END SUBROUTINE ESTPAC(MF,ML,EEE,PIN,QS,RS,DS,PS,IA,IB,IP,SUM,U,ZAV) INTEGER MF,ML,IA,IB,IP DOUBLE PRECISION EEE,PIN,SUM,U,ZAV DOUBLE PRECISION QS(ML),RS(ML),DS(ML),PS(ML) 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) + EIG*R(X))*Y(X) = 0 ON (A,B) C C FOR USER-SUPPLIED COEFFICIENT FUNCTIONS P, Q, AND R. IT IS C CALLED FROM SLEIGN, AND ALSO FROM ESTEIG WHEN NO INITIAL GUESS C IS PROVIDED BY THE USER. C C THE SUBROUTINE APPROXIMATES THE (TRAPEZOIDAL RULE) INTEGRAL OF C C SQRT((EIG*R+Q)/P) C C WHERE THE INTEGRAL IS TAKEN OVER THOSE X IN (A,B) FOR WHICH C C (EIG*R+Q)/P .GT. 0 C C ********** C .. LOCAL SCALARS .. INTEGER J DOUBLE PRECISION PSUM,RT,RTSAV,V,W,WSAV C .. C .. INTRINSIC FUNCTIONS .. INTRINSIC ABS,SIGN,SQRT C .. IA = MF IB = 80 IP = MF C C SUM ACCUMULATES THE INTEGRAL APPROXIMATION. U MEASURES THE TOTAL C LENGTH OF SUBINTERVALS WHERE (EIG*R+Q)/P .GT. 0.0. ZAV IS THE C AVERAGE VALUE OF SQRT((EIG*R+Q)*P) OVER THOSE SUBINTERVALS. C SUM = 0.0 U = 0.0 ZAV = 0.0 WSAV = EEE*RS(MF) + QS(MF) IF (WSAV.GT.0.0) THEN RTSAV = SQRT(WSAV) ELSE RTSAV = 0.0 END IF DO 10 J=MF+1,ML W = EEE*RS(J) + QS(J) IF (W.GT.0.0) THEN IF (J.GE.80) IB = J U = U + DS(J-1) RT = SQRT(W) ELSE RT = 0.0 IF (U.EQ.0.0 .AND. RTSAV.EQ.0.0 .AND. IA.LE.19) IA = IA + 1 END IF IF (W.EQ.0.0 .OR. WSAV.EQ.0.0 .OR. W.EQ.SIGN(W,WSAV)) THEN V = RT + RTSAV ELSE V = (W*RT+WSAV*RTSAV)/ABS(W-WSAV) END IF WSAV = W RTSAV = RT PSUM = DS(J-1)*V IF (PSUM.LT.(PIN-SUM)) IP = J SUM = SUM + PSUM IF (U.GT.0.0) ZAV = ZAV + PSUM*(PS(J)+PS(J-1)) 10 CONTINUE SUM = 0.5*SUM ZAV = 0.25*ZAV RETURN END SUBROUTINE ALFBET(XEND,INTAB,TT,COEF1,COEF2,EIG,P0,QF,OK, 1 VALUE,IFLAG,DERIV) INTEGER INTAB,IFLAG LOGICAL OK DOUBLE PRECISION XEND,TT,COEF1,COEF2,EIG,P0,QF,VALUE,DERIV C ********** C C THIS SUBROUTINE COMPUTES A BOUNDARY VALUE FOR A SPECIFIED ENDPOINT C OF THE INTERVAL FOR A STURM-LIOUVILLE PROBLEM IN THE FORM C C (P(X)*Y'(X))' + (Q(X) + EIG*R(X))*Y(X) = 0 ON (A,B) C C FOR USER-SUPPLIED COEFFICIENT FUNCTIONS P, Q, AND R. IT IS CALLED C FROM SLEIGN. BOTH NONSINGULAR AND SINGULAR ENDPOINTS ARE TREATED. C C SUBPROGRAMS CALLED C C USER-SUPPLIED ..... P,Q,R C C SLEIGN-SUPPLIED ... DXDT,EXTRAP C C ********** C .. SCALARS IN COMMON .. DOUBLE PRECISION Z C .. C .. LOCAL SCALARS .. LOGICAL LOGIC DOUBLE PRECISION C,CD,D,HH,ONE,PI,PX,QX,RX,T,TEMP,X C .. C .. EXTERNAL FUNCTIONS .. DOUBLE PRECISION P,Q,R EXTERNAL P,Q,R C .. C .. EXTERNAL SUBROUTINES .. EXTERNAL DXDT,EXTRAP C .. C .. INTRINSIC FUNCTIONS .. INTRINSIC ABS,ATAN,SIGN,SQRT C .. C .. COMMON BLOCKS .. COMMON /ZEE/Z C .. C SET MACHINE DEPENDENT CONSTANT. C C PI (VARIABLE ONE SET TO 1.0 EASES PRECISION CONVERSION). ONE = 1.0 PI = 4.0*ATAN(ONE) C IFLAG = 1 DERIV = 0.0 IF (OK) THEN VALUE = 0.5*PI IF (COEF1.NE.0.0) VALUE = ATAN(-Z*COEF2/COEF1) LOGIC = (TT.LT.0.0 .AND. VALUE.LT.0.0) .OR. 1 (TT.GT.0.0 .AND. VALUE.LE.0.0) IF (LOGIC) VALUE = VALUE + PI ELSE LOGIC = (INTAB.EQ.2 .AND. TT.GT.0.0) .OR. 1 (INTAB.EQ.3 .AND. TT.LT.0.0) .OR. 2 INTAB.EQ.4 .OR. (P0.GT.0.0 .AND. QF.LT.0.0) IF (LOGIC) THEN T = SIGN(ONE,TT) CALL EXTRAP(T,TT,EIG,VALUE,DERIV,IFLAG) ELSE CALL DXDT(TT,TEMP,X) PX = P(X)/Z QX = Q(X)/Z RX = R(X)/Z C = 2.0*(EIG*RX+QX) IF (C.LT.0.0) THEN VALUE = 0.0 IF (P0.GT.0.0) VALUE = 0.5*PI ELSE HH = ABS(XEND-X) D = 2.0*HH/PX CD = C*D*HH IF (P0.GT.0.0) THEN VALUE = C*HH IF (CD.LT.1.0) VALUE = VALUE/(1.0+SQRT(1.0-CD)) VALUE = VALUE + 0.5*PI ELSE VALUE = D IF (CD.LT.1.0) VALUE = VALUE/(1.0+SQRT(1.0-CD)) END IF END IF END IF END IF RETURN END SUBROUTINE F(T,Y,YP) DOUBLE PRECISION T DOUBLE PRECISION Y(2),YP(3) C ********** C C THIS SUBROUTINE EVALUATES THE DERIVATIVE FUNCTIONS FOR USE WITH C INTEGRATOR GERK IN SOLVING A STURM-LIOUVILLE PROBLEM IN THE FORM C C (P(X)*Y'(X))' + (Q(X) + EIG*R(X))*Y(X) = 0 ON (A,B) C C FOR USER-SUPPLIED COEFFICIENT FUNCTIONS P, Q, AND R. C C SUBPROGRAMS CALLED C C USER-SUPPLIED ..... P,Q,R C C SLEIGN-SUPPLIED ... DXDT C C ********** C .. SCALARS IN COMMON .. LOGICAL EIGF DOUBLE PRECISION EIG,Z C .. C .. LOCAL SCALARS .. DOUBLE PRECISION C,C2,DT,QX,RX,S,S2,V,W,X,XP,ZP C .. C .. EXTERNAL FUNCTIONS .. DOUBLE PRECISION P,Q,R EXTERNAL P,Q,R C .. C .. EXTERNAL SUBROUTINES .. EXTERNAL DXDT C .. C .. INTRINSIC FUNCTIONS .. INTRINSIC ABS,COS,SIN C .. C .. COMMON BLOCKS .. COMMON /DATAF/EIG,EIGF COMMON /ZEE/Z C .. DT = 1.0 X = T IF (.NOT.EIGF) CALL DXDT(T,DT,X) ZP = ABS(Z) QX = Q(X)/ZP RX = R(X)/ZP XP = ZP/P(X) V = EIG*RX + QX S = SIN(Y(1)) C = COS(Y(1)) S2 = S*S C2 = C*C YP(1) = DT*(XP*C2+V*S2) W = (XP-V)*S*C IF (Z.LT.0.0) RX = ABS(RX) YP(2) = DT*(-2.0*W*Y(2)+RX*S2) YP(3) = DT*W RETURN END SUBROUTINE DXDT(T,DT,X) DOUBLE PRECISION T,DT,X C ********** C C THIS SUBROUTINE TRANSFORMS COORDINATES FROM T ON (-1,1) TO C X ON (A,B) IN THE SOLUTION OF A STURM-LIOUVILLE PROBLEM. C IT IS CALLED FROM SUBROUTINES SLEIGN, ALFBET, F, AND EXTRAP. C C ********** C .. SCALARS IN COMMON .. INTEGER INTAB DOUBLE PRECISION A,B,C1,C2 C .. C .. LOCAL SCALARS .. DOUBLE PRECISION U C .. C .. INTRINSIC FUNCTIONS .. INTRINSIC ABS C .. C .. COMMON BLOCKS .. COMMON /DATADT/A,B,C1,C2,INTAB C .. U = C1*T + C2 GO TO (10,20,30,40), INTAB 10 CONTINUE DT = C1*0.5*(B-A) X = 0.5*((B+A)+(B-A)*U) RETURN 20 CONTINUE DT = C1*2.0/(1.0-U)**2 X = A + (1.0+U)/(1.0-U) RETURN 30 CONTINUE DT = C1*2.0/(1.0+U)**2 X = B - (1.0-U)/(1.0+U) RETURN 40 CONTINUE DT = C1/(1.0-ABS(U))**2 X = U/(1.0-ABS(U)) RETURN END SUBROUTINE EXTRAP(T,TT,EIG,VALUE,DERIV,IFLAG) INTEGER IFLAG DOUBLE PRECISION T,TT,EIG,VALUE,DERIV 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) + EIG*R(X))*Y(X) = 0 ON (A,B) C C FOR USER-SUPPLIED COEFFICIENT FUNCTIONS P, Q, AND R. C C EXTRAP, WHICH IN TURN CALLS INTPOL, EXTRAPOLATES THE FUNCTION C C ARCTAN(1.0/SQRT(-P*(EIG*R+Q))) C C FROM ITS VALUES FOR T WITHIN (-1,1) TO AN ENDPOINT. C C SUBPROGRAMS CALLED C C USER-SUPPLIED ..... P,Q,R C C SLEIGN-SUPPLIED ... DXDT,INTPOL C C ********** C .. SCALARS IN COMMON .. DOUBLE PRECISION Z C .. C .. LOCAL SCALARS .. INTEGER KGOOD DOUBLE PRECISION ANS,CTN,ERROR,PROD,PX,QX,RX,T1,TEMP,X C .. C .. LOCAL ARRAYS .. DOUBLE PRECISION FN1(5),XN(5) C .. C .. EXTERNAL FUNCTIONS .. DOUBLE PRECISION P,Q,R EXTERNAL P,Q,R C .. C .. EXTERNAL SUBROUTINES .. EXTERNAL DXDT,INTPOL C .. C .. INTRINSIC FUNCTIONS .. INTRINSIC ABS,ATAN,SQRT,TAN C .. C .. COMMON BLOCKS .. COMMON /ZEE/Z C .. IFLAG = 1 KGOOD = 0 T1 = TT 10 CONTINUE CALL DXDT(T1,TEMP,X) PX = P(X)/Z QX = Q(X)/Z RX = R(X)/Z PROD = -PX*(EIG*RX+QX) IF (PROD.LE.0.0) THEN T1 = 0.5*(T1+T) IF ((1.0+(T1-T)**2).GT.1.0) GO TO 10 IFLAG = 5 RETURN ELSE KGOOD = KGOOD + 1 XN(KGOOD) = T1 FN1(KGOOD) = ATAN(1.0/SQRT(PROD)) T1 = 0.5*(T+T1) IF (KGOOD.LT.5) GO TO 10 END IF T1 = 0.01 CALL INTPOL(5,XN,FN1,T,T1,3,ANS,ERROR) VALUE = ABS(ANS) CTN = 1.0/TAN(VALUE) DERIV = 0.5*PX*RX/CTN/(1.0+CTN**2) TT = XN(1) RETURN END SUBROUTINE INTPOL(N,XN,FN,X,ABSERR,MAXDEG,ANS,ERROR) INTEGER N,MAXDEG DOUBLE PRECISION X,ABSERR,ANS,ERROR DOUBLE PRECISION XN(N),FN(N) C ********** C C THIS SUBROUTINE FORMS AN INTERPOLATING POLYNOMIAL FOR DATA PAIRS. C IT IS CALLED FROM EXTRAP IN SOLVING A STURM-LIOUVILLE PROBLEM. C C ********** C .. LOCAL SCALARS .. INTEGER I,I1,II,IJ,IK,IKM1,J,K,L,LIMIT DOUBLE PRECISION PROD C .. C .. LOCAL ARRAYS .. INTEGER INDEX(10) DOUBLE PRECISION V(10,10) C .. C .. INTRINSIC FUNCTIONS .. INTRINSIC ABS,MIN 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.0 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 ************************************************************************ SUBROUTINE GERK(F, NEQN, Y, T, TOUT, RELERR, ABSERR, IFLAG, * GERROR, WORK, IWORK) C C FEHLBERG FOURTH(FIFTH) ORDER RUNGE-KUTTA METHOD WITH C GLOBAL ERROR ASSESSMENT C C WRITTEN BY H.A.WATTS AND L.F.SHAMPINE C SANDIA LABORATORIES C C GERK IS DESIGNED TO SOLVE SYSTEMS OF DIFFERENTIAL EQUATIONS C WHEN IT IS IMPORTANT TO HAVE A READILY AVAILABLE GLOBAL ERROR C ESTIMATE. PARALLEL INTEGRATION IS PERFORMED TO YIELD TWO C SOLUTIONS ON DIFFERENT MESH SPACINGS AND GLOBAL EXTRAPOLATION C IS APPLIED TO PROVIDE AN ESTIMATE OF THE GLOBAL ERROR IN THE C MORE ACCURATE SOLUTION. C C FOR IBM SYSTEM 360 AND 370 AND OTHER MACHINES OF SIMILAR C ARITHMETIC CHARACTERISTICS, THIS CODE SHOULD BE CONVERTED TO C DOUBLE PRECISION. C C******************************************************************* C ABSTRACT C******************************************************************* C C SUBROUTINE GERK INTEGRATES A SYSTEM OF NEQN FIRST ORDER C ORDINARY DIFFERENTIAL EQUATIONS OF THE FORM C DY(I)/DT = F(T,Y(1),Y(2),...,Y(NEQN)) C WHERE THE Y(I) ARE GIVEN AT T . C TYPICALLY THE SUBROUTINE IS USED TO INTEGRATE FROM T TO TOUT C BUT IT CAN BE USED AS A ONE-STEP INTEGRATOR TO ADVANCE THE C SOLUTION A SINGLE STEP IN THE DIRECTION OF TOUT. ON RETURN, AN C ESTIMATE OF THE GLOBAL ERROR IN THE SOLUTION AT T IS PROVIDED C AND THE PARAMETERS IN THE CALL LIST ARE SET FOR CONTINUING THE C INTEGRATION. THE USER HAS ONLY TO CALL GERK AGAIN (AND PERHAPS C DEFINE A NEW VALUE FOR TOUT). ACTUALLY, GERK IS MERELY AN C INTERFACING ROUTINE WHICH ALLOCATES VIRTUAL STORAGE IN THE C ARRAYS WORK, IWORK AND CALLS SUBROUTINE GERKS FOR THE SOLUTION. C GERKS IN TURN CALLS SUBROUTINE FEHL WHICH COMPUTES AN APPROX- C IMATE SOLUTION OVER ONE STEP. C C GERK USES THE RUNGE-KUTTA-FEHLBERG (4,5) METHOD DESCRIBED C IN THE REFERENCE C E.FEHLBERG , LOW-ORDER CLASSICAL RUNGE-KUTTA FORMULAS WITH C STEPSIZE CONTROL , NASA TR R-315 C C C THE PARAMETERS REPRESENT- C F -- SUBROUTINE F(T,Y,YP) TO EVALUATE DERIVATIVES C YP(I)=DY(I)/DT C NEQN -- NUMBER OF EQUATIONS TO BE INTEGRATED C Y(*) -- SOLUTION VECTOR AT T C T -- INDEPENDENT VARIABLE C TOUT -- OUTPUT POINT AT WHICH SOLUTION IS DESIRED C RELERR,ABSERR -- RELATIVE AND ABSOLUTE ERROR TOLERANCES FOR C LOCAL ERROR TEST. AT EACH STEP THE CODE REQUIRES THAT C ABS(LOCAL ERROR) .LE. RELERR*ABS(Y) + ABSERR C FOR EACH COMPONENT OF THE LOCAL ERROR AND SOLUTION C VECTORS. C IFLAG -- INDICATOR FOR STATUS OF INTEGRATION C GERROR(*) -- VECTOR WHICH ESTIMATES THE GLOBAL ERROR AT T. C THAT IS, GERROR(I) APPROXIMATES Y(I)-TRUE C SOLUTION(I). C WORK(*) -- ARRAY TO HOLD INFORMATION INTERNAL TO GERK WHICH C IS NECESSARY FOR SUBSEQUENT CALLS. MUST BE DIMENSIONED C AT LEAST 3+8*NEQN C IWORK(*) -- INTEGER ARRAY USED TO HOLD INFORMATION INTERNAL C TO GERK WHICH IS NECESSARY FOR SUBSEQUENT CALLS. MUST C BE DIMENSIONED AT LEAST 5 C C C******************************************************************* C FIRST CALL TO GERK C******************************************************************* C C THE USER MUST PROVIDE STORAGE IN HIS CALLING PROGRAM FOR THE C ARRAYS IN THE CALL LIST - Y(NEQN), WORK(3+8*NEQN), IWORK(5), C DECLARE F IN AN EXTERNAL STATEMENT, SUPPLY SUBROUTINE F(T,Y,YP) C AND INITIALIZE THE FOLLOWING PARAMETERS- C C NEQN -- NUMBER OF EQUATIONS TO BE INTEGRATED. (NEQN .GE. 1) C Y(*) -- VECTOR OF INITIAL CONDITIONS C T -- STARTING POINT OF INTEGRATION , MUST BE A VARIABLE C TOUT -- OUTPUT POINT AT WHICH SOLUTION IS DESIRED. C T=TOUT IS ALLOWED ON THE FIRST CALL ONLY,IN WHICH CASE C GERK RETURNS WITH IFLAG=2 IF CONTINUATION IS POSSIBLE. C RELERR,ABSERR -- RELATIVE AND ABSOLUTE LOCAL ERROR TOLERANCES C WHICH MUST BE NON-NEGATIVE BUT MAY BE CONSTANTS. WE CAN C USUALLY EXPECT THE GLOBAL ERRORS TO BE SOMEWHAT SMALLER C THAN THE REQUESTED LOCAL ERROR TOLERANCES. TO AVOID C LIMITING PRECISION DIFFICULTIES THE CODE ALWAYS USES C THE LARGER OF RELERR AND AN INTERNAL RELATIVE ERROR C PARAMETER WHICH IS MACHINE DEPENDENT. C IFLAG -- +1,-1 INDICATOR TO INITIALIZE THE CODE FOR EACH NEW C PROBLEM. NORMAL INPUT IS +1. THE USER SHOULD SET IFLAG= C -1 ONLY WHEN ONE-STEP INTEGRATOR CONTROL IS ESSENTIAL. C IN THIS CASE, GERK ATTEMPTS TO ADVANCE THE SOLUTION A C SINGLE STEP IN THE DIRECTION OF TOUT EACH TIME IT IS C CALLED. SINCE THIS MODE OF OPERATION RESULTS IN EXTRA C COMPUTING OVERHEAD, IT SHOULD BE AVOIDED UNLESS NEEDED. C C C******************************************************************* C OUTPUT FROM GERK C******************************************************************* C C Y(*) -- SOLUTION AT T C T -- LAST POINT REACHED IN INTEGRATION. C IFLAG = 2 -- INTEGRATION REACHED TOUT. INDICATES SUCCESSFUL C RETURN AND IS THE NORMAL MODE FOR CONTINUING C INTEGRATION. C =-2 -- A SINGLE SUCCESSFUL STEP IN THE DIRECTION OF C TOUT HAS BEEN TAKEN. NORMAL MODE FOR CONTINUING C INTEGRATION ONE STEP AT A TIME. C = 3 -- INTEGRATION WAS NOT COMPLETED BECAUSE MORE THAN C 9000 DERIVATIVE EVALUATIONS WERE NEEDED. THIS C IS APPROXIMATELY 500 STEPS. C = 4 -- INTEGRATION WAS NOT COMPLETED BECAUSE SOLUTION C VANISHED MAKING A PURE RELATIVE ERROR TEST C IMPOSSIBLE. MUST USE NON-ZERO ABSERR TO CONTINUE. C USING THE ONE-STEP INTEGRATION MODE FOR ONE STEP C IS A GOOD WAY TO PROCEED. C = 5 -- INTEGRATION WAS NOT COMPLETED BECAUSE REQUESTED C ACCURACY COULD NOT BE ACHIEVED USING SMALLEST C ALLOWABLE STEPSIZE. USER MUST INCREASE THE ERROR C TOLERANCE BEFORE CONTINUED INTEGRATION CAN BE C ATTEMPTED. C = 6 -- GERK IS BEING USED INEFFICIENTLY IN SOLVING C THIS PROBLEM. TOO MUCH OUTPUT IS RESTRICTING THE C NATURAL STEPSIZE CHOICE. USE THE ONE-STEP C INTEGRATOR MODE. C = 7 -- INVALID INPUT PARAMETERS C THIS INDICATOR OCCURS IF ANY OF THE FOLLOWING IS C SATISFIED - NEQN .LE. 0 C T=TOUT AND IFLAG .NE. +1 OR -1 C RELERR OR ABSERR .LT. 0. C IFLAG .EQ. 0 OR .LT. -2 OR .GT. 7 C GERROR(*) -- ESTIMATE OF THE GLOBAL ERROR IN THE SOLUTION AT T C WORK(*),IWORK(*) -- INFORMATION WHICH IS USUALLY OF NO C INTEREST TO THE USER BUT NECESSARY FOR SUBSEQUENT C CALLS. WORK(1),...,WORK(NEQN) CONTAIN THE FIRST C DERIVATIVES OF THE SOLUTION VECTOR Y AT T. C WORK(NEQN+1) CONTAINS THE STEPSIZE H TO BE C ATTEMPTED ON THE NEXT STEP. IWORK(1) CONTAINS C THE DERIVATIVE EVALUATION COUNTER. C C C******************************************************************* C SUBSEQUENT CALLS TO GERK C******************************************************************* C C SUBROUTINE GERK RETURNS WITH ALL INFORMATION NEEDED TO CONTINUE C THE INTEGRATION. IF THE INTEGRATION REACHED TOUT, THE USER NEED C ONLY DEFINE A NEW TOUT AND CALL GERK AGAIN. IN THE ONE-STEP C INTEGRATOR MODE (IFLAG=-2) THE USER MUST KEEP IN MIND THAT EACH C STEP TAKEN IS IN THE DIRECTION OF THE CURRENT TOUT. UPON C REACHING TOUT (INDICATED BY CHANGING IFLAG TO 2), THE USER MUST C THEN DEFINE A NEW TOUT AND RESET IFLAG TO -2 TO CONTINUE IN THE C ONE-STEP INTEGRATOR MODE. C C IF THE INTEGRATION WAS NOT COMPLETED BUT THE USER STILL WANTS C TO CONTINUE (IFLAG=3 CASE), HE JUST CALLS GERK AGAIN. THE C FUNCTION COUNTER IS THEN RESET TO 0 AND ANOTHER 9000 FUNCTION C EVALUATIONS ARE ALLOWED. C C HOWEVER, IN THE CASE IFLAG=4, THE USER MUST FIRST ALTER THE C ERROR CRITERION TO USE A POSITIVE VALUE OF ABSERR BEFORE C INTEGRATION CAN PROCEED. IF HE DOES NOT,EXECUTION IS TERMINATED. C C ALSO, IN THE CASE IFLAG=5, IT IS NECESSARY FOR THE USER TO C RESET IFLAG TO 2 (OR -2 WHEN THE ONE-STEP INTEGRATION MODE IS C BEING USED) AS WELL AS INCREASING EITHER ABSERR,RELERR OR BOTH C BEFORE THE INTEGRATION CAN BE CONTINUED. IF THIS IS NOT DONE, C EXECUTION WILL BE TERMINATED. THE OCCURRENCE OF IFLAG=5 C INDICATES A TROUBLE SPOT (SOLUTION IS CHANGING RAPIDLY, C SINGULARITY MAY BE PRESENT) AND IT OFTEN IS INADVISABLE TO C CONTINUE. C C IF IFLAG=6 IS ENCOUNTERED, THE USER SHOULD USE THE ONE-STEP C INTEGRATION MODE WITH THE STEPSIZE DETERMINED BY THE CODE. IF C THE USER INSISTS UPON CONTINUING THE INTEGRATION WITH GERK IN C THE INTERVAL MODE, HE MUST RESET IFLAG TO 2 BEFORE CALLING GERK C AGAIN. OTHERWISE,EXECUTION WILL BE TERMINATED. C C IF IFLAG=7 IS OBTAINED, INTEGRATION CAN NOT BE CONTINUED UNLESS C THE INVALID INPUT PARAMETERS ARE CORRECTED. C C IT SHOULD BE NOTED THAT THE ARRAYS WORK,IWORK CONTAIN C INFORMATION REQUIRED FOR SUBSEQUENT INTEGRATION. ACCORDINGLY, C WORK AND IWORK SHOULD NOT BE ALTERED. C C******************************************************************* C C .. SCALAR ARGUMENTS .. INTEGER IFLAG,NEQN DOUBLE PRECISION ABSERR,RELERR,T,TOUT C .. C .. ARRAY ARGUMENTS .. INTEGER IWORK(5) DOUBLE PRECISION GERROR(NEQN),WORK(3+8*NEQN),Y(NEQN) C .. C .. SUBROUTINE ARGUMENTS .. EXTERNAL F C .. C .. LOCAL SCALARS .. INTEGER K1,K1M,K2,K3,K4,K5,K6,K7,K8 C .. C .. EXTERNAL SUBROUTINES .. EXTERNAL GERKS C .. C COMPUTE INDICES FOR THE SPLITTING OF THE WORK ARRAY K1M = NEQN + 1 K1 = K1M + 1 K2 = K1 + NEQN K3 = K2 + NEQN K4 = K3 + NEQN K5 = K4 + NEQN K6 = K5 + NEQN K7 = K6 + NEQN K8 = K7 + NEQN C ******************************************************************* C THIS INTERFACING ROUTINE MERELY RELIEVES THE USER OF A LONG C CALLING LIST VIA THE SPLITTING APART OF TWO WORKING STORAGE C ARRAYS. IF THIS IS NOT COMPATIBLE WITH THE USERS COMPILER, C HE MUST USE GERKS DIRECTLY. C ******************************************************************* CALL GERKS(F, NEQN, Y, T, TOUT, RELERR, ABSERR, IFLAG, * GERROR, WORK(1), WORK(K1M), WORK(K1), WORK(K2), WORK(K3), * WORK(K4), WORK(K5), WORK(K6), WORK(K7), WORK(K8), * WORK(K8+1), IWORK(1), IWORK(2), IWORK(3), IWORK(4), IWORK(5)) RETURN END SUBROUTINE GERKS(F, NEQN, Y, T, TOUT, RELERR, ABSERR, IFLAG, * GERROR, YP, H, F1, F2, F3, F4, F5, YG, YGP, SAVRE, SAVAE, * NFE, KOP, INIT, JFLAG, KFLAG) C FEHLBERG FOURTH(FIFTH) ORDER RUNGE-KUTTA METHOD WITH C GLOBAL ERROR ASSESSMENT C ******************************************************************* C GERKS INTEGRATES A SYSTEM OF FIRST ORDER ORDINARY DIFFERENTIAL C EQUATIONS AS DESCRIBED IN THE COMMENTS FOR GERK. THE ARRAYS C YP,F1,F2,F3,F4,F5,YG AND YGP (OF DIMENSION AT LEAST NEQN) AND C THE VARIABLES H,SAVRE,SAVAE,NFE,KOP,INIT,JFLAG,AND KFLAG ARE C USED INTERNALLY BY THE CODE AND APPEAR IN THE CALL LIST TO C ELIMINATE LOCAL RETENTION OF VARIABLES BETWEEN CALLS. C ACCORDINGLY, THEY SHOULD NOT BE ALTERED. ITEMS OF POSSIBLE C INTEREST ARE C YP - DERIVATIVE OF SOLUTION VECTOR AT T C H - AN APPROPRIATE STEPSIZE TO BE USED FOR THE NEXT STEP C NFE- COUNTER ON THE NUMBER OF DERIVATIVE FUNCTION C EVALUATIONS. C ******************************************************************* C .. SCALAR ARGUMENTS .. INTEGER IFLAG,INIT,JFLAG,KFLAG,KOP,NEQN,NFE DOUBLE PRECISION ABSERR,H,RELERR,SAVAE,SAVRE,T,TOUT C .. C .. ARRAY ARGUMENTS .. DOUBLE PRECISION F1(NEQN),F2(NEQN),F3(NEQN),F4(NEQN),F5(NEQN), 1 GERROR(NEQN),Y(NEQN),YG(NEQN),YGP(NEQN),YP(NEQN) C .. C .. SUBROUTINE ARGUMENTS .. EXTERNAL F C .. C .. LOCAL SCALARS .. INTEGER K,MAXNFE,MFLAG LOGICAL HFAILD,OUTPUT DOUBLE PRECISION A,AE,DT,EE,EEOET,ESTTOL,ET,HH,HMIN,REMIN,RER,S, 1 SCALE,TOL,TOLN,TS,U,U26,YPK C .. C .. EXTERNAL SUBROUTINES .. EXTERNAL FEHL C .. C .. INTRINSIC FUNCTIONS .. INTRINSIC ABS,MAX,SIGN C .. C ******************************************************************* C ********** WARNING - THE FOLLOWING CONSTANT IS MACHINE DEPENDENT C (ACTIVATED VALUE IS FOR IEEE MACHINES): C U - THE COMPUTER UNIT ROUNDOFF ERROR U IS THE SMALLEST POSITIVE C VALUE REPRESENTABLE IN THE MACHINE SUCH THAT 1.+ U .GT. 1. C C CONSTANTS FOR SOME MACHINES (ACTIVATE BY REMOVING 'C' FROM COLUMN 1) C --------------------------------------------------------------------- C VAX 11/780 C DATA U /1.4D-17/ C IEEE (ALLIANT FX/8, ENCORE MULTIMAX, SEQUENT BALANCE, SUN, ETC.) DATA U /2.2D-16/ C IBM 3033 C DATA U /2.2D-16/ C CRAY C DATA U /2.5D-29/ C --------------------------------------------------------------------- C ******************************************************************* C REMIN IS A TOLERANCE THRESHOLD WHICH IS ALSO DETERMINED BY THE C INTEGRATION METHOD. IN PARTICULAR, A FIFTH ORDER METHOD WILL C GENERALLY NOT BE CAPABLE OF DELIVERING ACCURACIES NEAR LIMITING C PRECISION ON COMPUTERS WITH LONG WORDLENGTHS. DATA REMIN /3.D-11/ C ******************************************************************* C THE EXPENSE IS CONTROLLED BY RESTRICTING THE NUMBER C OF FUNCTION EVALUATIONS TO BE APPROXIMATELY MAXNFE. C AS SET,THIS CORRESPONDS TO ABOUT 500 STEPS. DATA MAXNFE /9000/ C ******************************************************************* C CHECK INPUT PARAMETERS IF (NEQN.LT.1) GO TO 10 IF ((RELERR.LT.0.) .OR. (ABSERR.LT.0.)) GO TO 10 MFLAG = ABS(IFLAG) IF ((MFLAG.GE.1) .AND. (MFLAG.LE.7)) GO TO 20 C INVALID INPUT 10 IFLAG = 7 RETURN C IS THIS THE FIRST CALL 20 IF (MFLAG.EQ.1) GO TO 70 C CHECK CONTINUATION POSSIBILITIES IF (T.EQ.TOUT) GO TO 10 IF (MFLAG.NE.2) GO TO 30 C IFLAG = +2 OR -2 IF (INIT.EQ.0) GO TO 60 IF (KFLAG.EQ.3) GO TO 50 IF ((KFLAG.EQ.4) .AND. (ABSERR.EQ.0.)) GO TO 40 IF ((KFLAG.EQ.5) .AND. (RELERR.LE.SAVRE) .AND. * (ABSERR.LE.SAVAE)) GO TO 40 GO TO 70 C IFLAG = 3,4,5,6, OR 7 30 IF (IFLAG.EQ.3) GO TO 50 IF ((IFLAG.EQ.4) .AND. (ABSERR.GT.0.)) GO TO 60 C INTEGRATION CANNOT BE CONTINUED SINCE USER DID NOT RESPOND TO C THE INSTRUCTIONS PERTAINING TO IFLAG=4,5,6 OR 7 40 STOP C ******************************************************************* C RESET FUNCTION EVALUATION COUNTER 50 NFE = 0 IF (MFLAG.EQ.2) GO TO 70 C RESET FLAG VALUE FROM PREVIOUS CALL 60 IFLAG = JFLAG C SAVE INPUT IFLAG AND SET CONTINUATION FLAG VALUE FOR SUBSEQUENT C INPUT CHECKING 70 JFLAG = IFLAG KFLAG = 0 C SAVE RELERR AND ABSERR FOR CHECKING INPUT ON SUBSEQUENT CALLS SAVRE = RELERR SAVAE = ABSERR C RESTRICT RELATIVE ERROR TOLERANCE TO BE AT LEAST AS LARGE AS C 32U+REMIN TO AVOID LIMITING PRECISION DIFFICULTIES ARISING C FROM IMPOSSIBLE ACCURACY REQUESTS RER = MAX(RELERR,32.*U+REMIN) U26 = 26.*U DT = TOUT - T IF (MFLAG.EQ.1) GO TO 80 IF (INIT.EQ.0) GO TO 90 GO TO 110 C ******************************************************************* C INITIALIZATION -- C SET INITIALIZATION COMPLETION INDICATOR,INIT C SET INDICATOR FOR TOO MANY OUTPUT POINTS,KOP C EVALUATE INITIAL DERIVATIVES C COPY INITIAL VALUES AND DERIVATIVES FOR THE C PARALLEL SOLUTION C SET COUNTER FOR FUNCTION EVALUATIONS,NFE C ESTIMATE STARTING STEPSIZE 80 INIT = 0 KOP = 0 A = T CALL F(A, Y, YP) NFE = 1 IF (T.NE.TOUT) GO TO 90 IFLAG = 2 RETURN 90 INIT = 1 H = ABS(DT) TOLN = 0. DO 100 K=1,NEQN YG(K) = Y(K) YGP(K) = YP(K) TOL = RER*ABS(Y(K)) + ABSERR IF (TOL.LE.0.) GO TO 100 TOLN = TOL YPK = ABS(YP(K)) IF (YPK*H**5.GT.TOL) H = (TOL/YPK)**0.2 100 CONTINUE IF (TOLN.LE.0.) H = 0. H = MAX(H,U26*MAX(ABS(T),ABS(DT))) C ******************************************************************* C SET STEPSIZE FOR INTEGRATION IN THE DIRECTION FROM T TO TOUT 110 H = SIGN(H,DT) C TEST TO SEE IF GERK IS BEING SEVERELY IMPACTED BY TOO MANY C OUTPUT POINTS IF (ABS(H).GT.2.*ABS(DT)) KOP = KOP + 1 IF (KOP.NE.100) GO TO 120 KOP = 0 IFLAG = 6 RETURN 120 IF (ABS(DT).GT.U26*ABS(T)) GO TO 140 C IF TOO CLOSE TO OUTPUT POINT,EXTRAPOLATE AND RETURN DO 130 K=1,NEQN YG(K) = YG(K) + DT*YGP(K) Y(K) = Y(K) + DT*YP(K) 130 CONTINUE A = TOUT CALL F(A, YG, YGP) CALL F(A, Y, YP) NFE = NFE + 2 GO TO 230 C INITIALIZE OUTPUT POINT INDICATOR 140 OUTPUT = .FALSE. C TO AVOID PREMATURE UNDERFLOW IN THE ERROR TOLERANCE FUNCTION, C SCALE THE ERROR TOLERANCES SCALE = 2./RER AE = SCALE*ABSERR C ******************************************************************* C ******************************************************************* C STEP BY STEP INTEGRATION 150 HFAILD = .FALSE. C SET SMALLEST ALLOWABLE STEPSIZE HMIN = U26*ABS(T) C ADJUST STEPSIZE IF NECESSARY TO HIT THE OUTPUT POINT. C LOOK AHEAD TWO STEPS TO AVOID DRASTIC CHANGES IN THE STEPSIZE C AND THUS LESSEN THE IMPACT OF OUTPUT POINTS ON THE CODE. DT = TOUT - T IF (ABS(DT).GE.2.*ABS(H)) GO TO 170 IF (ABS(DT).GT.ABS(H)) GO TO 160 C THE NEXT SUCCESSFUL STEP WILL COMPLETE THE INTEGRATION TO THE C OUTPUT POINT OUTPUT = .TRUE. H = DT GO TO 170 160 H = 0.5*DT C ******************************************************************* C CORE INTEGRATOR FOR TAKING A SINGLE STEP C ******************************************************************* C THE TOLERANCES HAVE BEEN SCALED TO AVOID PREMATURE UNDERFLOW C IN COMPUTING THE ERROR TOLERANCE FUNCTION ET. C TO AVOID PROBLEMS WITH ZERO CROSSINGS, RELATIVE ERROR IS C MEASURED USING THE AVERAGE OF THE MAGNITUDES OF THE SOLUTION C AT THE BEGINNING AND END OF A STEP. C THE ERROR ESTIMATE FORMULA HAS BEEN GROUPED TO CONTROL LOSS OF C SIGNIFICANCE. C TO DISTINGUISH THE VARIOUS ARGUMENTS, H IS NOT PERMITTED C TO BECOME SMALLER THAN 26 UNITS OF ROUNDOFF IN T. C PRACTICAL LIMITS ON THE CHANGE IN THE STEPSIZE ARE ENFORCED TO C SMOOTH THE STEPSIZE SELECTION PROCESS AND TO AVOID EXCESSIVE C CHATTERING ON PROBLEMS HAVING DISCONTINUITIES. C TO PREVENT UNNECESSARY FAILURES, THE CODE USES 9/10 THE C STEPSIZE IT ESTIMATES WILL SUCCEED. C AFTER A STEP FAILURE, THE STEPSIZE IS NOT ALLOWED TO INCREASE C FOR THE NEXT ATTEMPTED STEP. THIS MAKES THE CODE MORE C EFFICIENT ON PROBLEMS HAVING DISCONTINUITIES AND MORE C EFFECTIVE IN GENERAL SINCE LOCAL EXTRAPOLATION IS BEING USED C AND THE ERROR ESTIMATE MAY BE UNRELIABLE OR UNACCEPTABLE WHEN C A STEP FAILS. C ******************************************************************* C TEST NUMBER OF DERIVATIVE FUNCTION EVALUATIONS. C IF OKAY,TRY TO ADVANCE THE INTEGRATION FROM T TO T+H 170 IF (NFE.LE.MAXNFE) GO TO 180 C TOO MUCH WORK IFLAG = 3 KFLAG = 3 RETURN C ADVANCE AN APPROXIMATE SOLUTION OVER ONE STEP OF LENGTH H 180 CALL FEHL(F, NEQN, YG, T, H, YGP, F1, F2, F3, F4, F5, F1) NFE = NFE + 5 C COMPUTE AND TEST ALLOWABLE TOLERANCES VERSUS LOCAL ERROR C ESTIMATES AND REMOVE SCALING OF TOLERANCES. NOTE THAT RELATIVE C ERROR IS MEASURED WITH RESPECT TO THE AVERAGE MAGNITUDES OF THE C OF THE SOLUTION AT THE BEGINNING AND END OF THE STEP. EEOET = 0. DO 200 K=1,NEQN ET = ABS(YG(K)) + ABS(F1(K)) + AE IF (ET.GT.0.) GO TO 190 C INAPPROPRIATE ERROR TOLERANCE IFLAG = 4 KFLAG = 4 RETURN 190 EE = ABS((-2090.*YGP(K)+(21970.*F3(K)-15048.*F4(K))) * +(22528.*F2(K)-27360.*F5(K))) EEOET = MAX(EEOET,EE/ET) 200 CONTINUE ESTTOL = ABS(H)*EEOET*SCALE/752400. IF (ESTTOL.LE.1.) GO TO 210 C UNSUCCESSFUL STEP C REDUCE THE STEPSIZE , TRY AGAIN C THE DECREASE IS LIMITED TO A FACTOR OF 1/10 HFAILD = .TRUE. OUTPUT = .FALSE. S = 0.1 IF (ESTTOL.LT.59049.) S = 0.9/ESTTOL**0.2 H = S*H IF (ABS(H).GT.HMIN) GO TO 170 C REQUESTED ERROR UNATTAINABLE AT SMALLEST ALLOWABLE STEPSIZE IFLAG = 5 KFLAG = 5 RETURN C SUCCESSFUL STEP C STORE ONE-STEP SOLUTION YG AT T+H C AND EVALUATE DERIVATIVES THERE 210 TS = T T = T + H DO 220 K=1,NEQN YG(K) = F1(K) 220 CONTINUE A = T CALL F(A, YG, YGP) NFE = NFE + 1 C NOW ADVANCE THE Y SOLUTION OVER TWO STEPS OF C LENGTH H/2 AND EVALUATE DERIVATIVES THERE HH = 0.5*H CALL FEHL(F, NEQN, Y, TS, HH, YP, F1, F2, F3, F4, F5, Y) TS = TS + HH A = TS CALL F(A, Y, YP) CALL FEHL(F, NEQN, Y, TS, HH, YP, F1, F2, F3, F4, F5, Y) A = T CALL F(A, Y, YP) NFE = NFE + 12 C CHOOSE NEXT STEPSIZE C THE INCREASE IS LIMITED TO A FACT-6R OF 5 C IF STEP FAILURE HAS JUST OCCURRED, NEXT C STEPSIZE IS NOT ALLOWED TO INCREASE S = 5. IF (ESTTOL.GT.1.889568E-4) S = 0.9/ESTTOL**0.2 IF (HFAILD .AND. (S.GT.1.)) S=1. H = SIGN(MAX(S*ABS(H),HMIN),H) C ******************************************************************* C END OF CORE INTEGRATOR C ******************************************************************* C SHOULD WE TAKE ANOTHER STEP IF (OUTPUT) GO TO 230 IF (IFLAG.GT.0) GO TO 150 C ******************************************************************* C ******************************************************************* C INTEGRATION SUCCESSFULLY COMPLETED C ONE-STEP MODE IFLAG = -2 GO TO 240 C INTERVAL MODE 230 T = TOUT IFLAG = 2 240 DO 250 K=1,NEQN GERROR(K) = (YG(K)-Y(K))/31. 250 CONTINUE RETURN END SUBROUTINE FEHL(F, NEQN, Y, T, H, YP, F1, F2, F3, F4, F5, S) C FEHLBERG FOURTH-FIFTH ORDER RUNGE-KUTTA METHOD C ******************************************************************* C FEHL INTEGRATES A SYSTEM OF NEQN FIRST ORDER C ORDINARY DIFFERENTIAL EQUATIONS OF THE FORM C DY(I)/DT=F(T,Y(1),---,Y(NEQN)) C WHERE THE INITIAL VALUES Y(I) AND THE INITIAL DERIVATIVES C YP(I) ARE SPECIFIED AT THE STARTING POINT T. FEHL ADVANCES C THE SOLUTION OVER THE FIXED STEP H AND RETURNS C THE FIFTH ORDER (SIXTH ORDER ACCURATE LOCALLY) SOLUTION C APPROXIMATION AT T+H IN ARRAY S(I). C F1,---,F5 ARE ARRAYS OF DIMENSION NEQN WHICH ARE NEEDED C FOR INTERNAL STORAGE. C THE FORMULAS HAVE BEEN GROUPED TO CONTROL LOSS OF SIGNIFICANCE. C FEHL SHOULD BE CALLED WITH AN H NOT SMALLER THAN 13 UNITS OF C ROUNDOFF IN T SO THAT THE VARIOUS INDEPENDENT ARGUMENTS CAN BE C DISTINGUISHED. C ******************************************************************* C .. SCALAR ARGUMENTS .. INTEGER NEQN DOUBLE PRECISION H,T C .. C .. ARRAY ARGUMENTS .. DOUBLE PRECISION F1(NEQN),F2(NEQN),F3(NEQN),F4(NEQN),F5(NEQN), 1 S(NEQN),Y(NEQN),YP(NEQN) C .. C .. SUBROUTINE ARGUMENTS .. EXTERNAL F C .. C .. LOCAL SCALARS .. INTEGER K DOUBLE PRECISION CH C .. CH = 0.25*H DO 10 K=1,NEQN F5(K) = Y(K) + CH*YP(K) 10 CONTINUE CALL F(T+0.25*H, F5, F1) CH = 0.09375*H DO 20 K=1,NEQN F5(K) = Y(K) + CH*(YP(K)+3.*F1(K)) 20 CONTINUE CALL F(T+0.375*H, F5, F2) CH = H/2197. DO 30 K=1,NEQN F5(K) = Y(K) + CH*(1932.*YP(K)+(7296.*F2(K)-7200.*F1(K))) 30 CONTINUE CALL F(T+12./13.*H, F5, F3) CH = H/4104. DO 40 K=1,NEQN F5(K) = Y(K) + CH*((8341.*YP(K)-845.*F3(K))+(29440.*F2(K) * -32832.*F1(K))) 40 CONTINUE CALL F(T+H, F5, F4) CH = H/20520. DO 50 K=1,NEQN F1(K) = Y(K) + CH*((-6080.*YP(K)+(9295.*F3(K)-5643.*F4(K))) * +(41040.*F1(K)-28352.*F2(K))) 50 CONTINUE CALL F(T+0.5*H, F1, F5) C COMPUTE APPROXIMATE SOLUTION AT T+H CH = H/7618050. DO 60 K=1,NEQN S(K) = Y(K) + CH*((902880.*YP(K)+(3855735.*F3(K)-1371249.* * F4(K)))+(3953664.*F2(K)+277020.*F5(K))) 60 CONTINUE RETURN END ************************************************************************