SUBROUTINE DMRODE(F, A, B, HMIN, H, HMAX, N, X0, WK, EPS, DMR 10 * JSTART, AA, TT, XX, ZZ, FINIT, ALPHA, WKA, IER, IQ) C C F -AUTOMATIC STEP CHANGE MERSON DIFFERENTIAL EQUATION C U - SOLVER MODIFIED TO HANDLE FUNCTIONAL DIFFERENTIAL C N - EQUATIONS OF THE FORM C C - (1) DX/DT =F(T,X,Z,I) C T - WHERE F (DESCRIBED BELOW) IS A VECTOR VALUED C I - FUNCTION OF DIMENSION N, X AND Z ARE C O - N VECTORS SUCH THAT THE ITH COMPONENT OF Z IS C N - THE ITH COMPONENT OF X EVALUATED AT ALPHA(X,T,I) C - ALPHA IS A VECTOR VALUED FUNCTION (DESCRIBED C - BELOW). ALPHA IS A USER SUPPLIED STATE DEPENDENT C - RETARDING FUNCTION, AND (1) IS CALLED A C - RETARDED ORDINARY DIFF. EQ. (RODE) WITH STATE C - DEPENDENT LAG (SDL). NOTE -AS WRITTEN ABOVE C - EACH COMPONENT OF DX/DT MAY DEPEND ON A DIFFERENT C - RETARDING FUNCTION. IT IS ALSO POSSIBLE TO ALLOW C - A GIVEN COMPONENT TO DEPEND ON MORE THAN ONE C - RETARDING FUNCTION BY AUGMENTING THE ORDER OF C - THE SYSTEM (SEE DISCRIPTION FOR DETAILS). C P F-USER SUPPLIED EXTERNAL FUNCTION. THE DERIVATIVE C A - OF THE ITH COMPONENT OF X IS F(T, X, Z,I), WHERE C R - Z AND X ARE ARRAYS OF DIMENSION N. X(I) C A - IS THE SOLUTION OF THE ITH COMPONENET AT T AND C M - Z (I) IS THE ITH COMPONENT OF X AT ALPHA(X,T,I) C E - Z(I) IS CALCULATED INTERNALLY WITHIN DMRODE AND C T - PASSED TO F (SEE ALPHA BELOW). T IS THE CURRENT C E - VALUE OF THE INDEPENDENT VARIABLE. C R A-STARTING POINT OF INTEGRATION C S B-FINAL POINT OF INTEGRATION (NOT NEC. .GT. A) C HMIN-THE ABSOLUTE VALUE OF MINIMUM ALLOWABLE STEP SIZE C - THE STEP SIZE WILL NEVER BE SMALLER THAN HMIN C - IN ABSOLUTE VALUE EXCEPT AT THE END OF THE C -INTERVAL WHERE H IS ALWAYS .GT. HMIN/2 IN ABS.. C H-INITIAL GUESS OF THE INTEGRATION STEP SIZE. THE C - ADJUSTED STEPSIZE IS RETURNED IN H. C HMAX-THE ABSOLUTE VALUE OF MAXIMUM ALLOWABLE STEP SIZE C N-NUMBER OF SIMULTANEOUS EQUATIONS C X0-VECTOR OF LENGTH N CONTAINING THE VALUE OF X(A). C - X(B) IS RETURNED IN X0. C WK-WORK AREA OF DIMENSION .GE. 4*N C EPS-INPUT REQUEST OF RELATIVE ERROR IN LARGEST COMP. C JSTART-A PARAMETER USED TO INTIALIZE Z ARRAY INTERNALLY C - IN DMRODE. JSTART MUST BE SET EQUAL TO ZERO ON C - FIRST CALL TO DMRODE. TO INTEGRATE BEYOND B AFTER C - FIRST CALL CERTAIN ARRAYS MUST BE SAVED. THIS C - IS DONE AUTOMATICALLY BY LETTING JSTART=1. C - IF JSTART=-1, THEN THE INITIAL CONDITIONS AT C - THE PREVIOUS CALL ARE REINSTATED AND THE INTE- C - GRATION FROM A TO POSSIBLY A NEW CHOICE FOR B C - CAN BE EXECUTED. THIS ALLOWS INTERATIVE USE OF C - DMRODE IN ORDER TO CALCULATE DERIVATIVE JUMP C - POINTS IN THE SOLUTION(SEE DESCRIPTION) C AA-A PARAMETER THAT MUST BE SET EQUAL TO THE VALUE OF C - THE STARTING POINT A WHEN FIRST CALL TO DMRODE C - WAS MADE (I.E. WHEN JSTART=0) (SEE FINIT) C TT-STORAGE ARRAY FOR SAVING TIME STEP VALUES C - TT IS OF DIMENSION IQ WHERE IQ IS .GT.LENGTH OF C - ENTIRE INTERVAL/HMIN C XX,ZZ-STORAGE ARRAYS OF DIMENSION N BY IQ USED FOR C - SAVING VALUES OF X AND Z DESCRIBED IN F. C FINIT-USER SUPPLIED EXTERNAL FUNCTION, FINIT(T,I). FINIT C - RETURNS THE VALUE OF THE ITH COMPONENT OF C - THE INITIAL FUNCTION AT T WHENEVER ((T. LE.AA) C - .AND.(H.GE.0)).OR.((T.GE.AA).AND.(H.LE.0)). C ALPHA-USER SUPPLIED EXTERNAL FUNCTION OF THE C - FORM ALPHA(X,T,I) WHERE X IS ASSUMED TO BE C - THE CURRENT SOLUTION (SUPPLIED BY DMRODE) AND C - T THE CURRENT TIME. IN GENERAL ALPHA WILL BE C - DIFFERENT FOR EACH COMPONENT X(I). IN THE EVENT C - A RODE IS POSED WITH MORE THAN ONE LAG IN ANY C - GIVEN COMPONENT THE SYSTEM MAY BE AUGMENTED TO C - HANDLE THIS SITUATION (SEE DISCRIPTION). C WKA-WORK AREA OF DIMENSION .GE.11*N UNLESS C - IT IS KNOWN THE LAG T-ALPHA(T,X(T),I) IN ALL C - COMPONENTS WILL ALWAYS EXCEED HMAX OR VANISH C - OVER THE ENTIRE INTERVAL IN A GIVEN COMPONENT, C - IN WHICH CASE A DIMENSION OF 3*N WILL SUFFICE C IER-OUTPUT PARAMETER IF IER=0 AND HMAX.NE.HMIN C - NORMAL ERROR CHECKING TOOK PLACE. IF IER=1 C - HMIN WAS ATTAINED BY THE STEP CHANGING C - PROCEDURE, HENCE REQUESTED ACCURACY IS NOT C - GUARANTEED C IQ-THE DIMENSION OF TT IN THE CALLING PROGRAM (ALSO C - THE DIMENSION OF THE 2ND ARGUMENT OF THE XX AND C - ZZ ARRAYS). C PRECISION - SINGLE C AUTHOR/IMPLEMENTOR- KENNETH W. NEVES C REQUIRED ROUTINES C - DZETA, AN INTERPOLATION SCHEME OF ORDER COMP- C - ARABLE TO 4TH ORDER MERSON METHOD. DZETA IN TURN C - WILL CALL STARTER METHOD WHEN LAG IS SMALLER C - THAN CURRENT STEP SIZE. C DIMENSION WK(1), X0(1) DIMENSION WKA(1) DIMENSION TT(1), XX(N,1), ZZ(N,1) EXTERNAL F, FINIT, ALPHA INTEGER SW LOGICAL BE, BH, BR, BX, BT DATA ZERO, P5, OP5, THREE, FOUR /0.,.5,1.5,3.,4./ E5 = .5*EPS IER = 0 IF (A.EQ.B) GO TO 290 IB1 = N + N IB2 = IB1 + N C ************************************ C CHECK FOR THE PROPER SIGN OF H C ************************************ H = SIGN(ABS(H),B-A) C ************************************ C CHECK FOR 1ST CALL TO DMRODE. IF C YES, THEN INITIALIZE NEC. ARRAYS C ************************************ IF (JSTART.NE.0) GO TO 20 JRAY = 1 JSAVE = JRAY TT(1) = A DO 10 J=1,N XX(J,1) = X0(J) ZZ(J,1)=DZETA(TT(1),X0,TT,XX,ZZ,AA,F,FINIT,ALPHA, * N,J,WKA(N+1),H,1,IQ) 10 CONTINUE 20 CONTINUE BH = .TRUE. BR = .TRUE. BX = .TRUE. BT = .TRUE. X = A XS = A IF (ABS(H).GE.HMAX) H = SIGN(1.,H)*HMAX DO 30 J=1,N IJK0 = N + J WK(IJK0) = X0(J) 30 CONTINUE C ************************************ C CHECK JSTART AND EXECUTE C PROPER INITIALIZATION C ************************************ IF (JSTART) 40, 90, 50 40 JRAY = JSAVE 50 DO 60 J=1,N IJK0 = N + J WK(IJK0) = XX(J,JRAY) 60 CONTINUE JSAVE = JRAY GO TO 90 70 XS = X C ************************************ C AFTER EACH SUCCESSFUL TIME STEP C UPDATE TT,XX,ZZ ARRAYS. CALL TO C DZETA PERFORMS NECESSARY INTERPO- C LATION AT LAG POINT. C ************************************ TT(JRAY+1) = X DO 80 J=1,N IJK0 = N + J WK(IJK0) = X0(J) XX(J,JRAY+1) = X0(J) HH = X - TT(JRAY) ZZ(J,JRAY+1) = DZETA(TT(JRAY+1),X0,TT,XX,ZZ,AA,F,FINIT, * ALPHA,N,J,WKA(N+1),HH,JRAY,IQ) 80 CONTINUE C ************************************ C INSERT CALL TO PRINT ROUTINE HERE C IF DESIRED. C ************************************ JRAY = JRAY + 1 IF (BR) GO TO 90 H = HSS RETURN 90 HSS = H C ************************************ C TEST FOR END OF PROBLEM C ************************************ Q = X + H - B BE = .TRUE. IF (.NOT.((H.GT.ZERO .AND. Q/H.GE.-.5) .OR. (H.LT.ZERO .AND. * Q/(-H).LE..5))) GO TO 110 IF (BT) HSS = H H = .5*(B-X) IF (.NOT.BT) H = B - X IF (BT) GO TO 100 BR = .FALSE. 100 BT = .FALSE. 110 CONTINUE H3 = H/THREE C ************************************ C CALCULATE SOLN. AT X+H C ************************************ DO 270 SW=1,5 C ************************************ C DZETA CALCULATES LAG AND SOLUTION C AT PAST POINT (NEC. FOR F ) C ************************************ DO 120 JS=1,N WKA(JS) = DZETA(X,X0,TT,XX,ZZ,AA,F,FINIT,ALPHA,N,JS, * WKA(N+1),H,JRAY,IQ) 120 CONTINUE DO 130 JS=1,N WK(JS) = F(X,X0,WKA,JS) 130 CONTINUE DO 230 I=1,N Q = H3*WK(I) IJK0 = N + I IJK1 = IB1 + I IJK2 = IB2 + I GO TO (140, 150, 160, 170, 180), SW 140 R = Q WK(IJK1) = Q GO TO 190 C 150 R = P5*(Q+WK(IJK1)) GO TO 190 C 160 R = THREE*Q WK(IJK2) = R R = .375*(R+WK(IJK1)) GO TO 190 C 170 R = WK(IJK1) + FOUR*Q WK(IJK1) = R R = OP5*(R-WK(IJK2)) GO TO 190 C 180 R = P5*(Q+WK(IJK1)) Q = ABS(R+R-OP5*(Q+WK(IJK2))) 190 X0(I) = WK(IJK0) + R IF (SW.NE.5) GO TO 230 C ************************************ C AUTOMATIC STEP CHANGE C ************************************ E = ABS(X0(I)) R = E5 IF (E.GE.1.E-3) R = E*E5 IF (HMIN.EQ.HMAX) GO TO 280 C ************************************ C TEST ADJUSTMENT OF THE STEP C ************************************ IF (Q.LT.R .OR. (.NOT.BX)) GO TO 220 BR = .TRUE. BT = .TRUE. BH = .FALSE. H = P5*H IF (ABS(H).GT.HMIN) GO TO 200 IF (.NOT.BR) GO TO 280 C ************************************ C THE STEP IS HALVED RESTORE X AND X0, C AND GO BACK FOR REPEATED INTEGRATION C WITH THIS NEW STEP C ************************************ H = SIGN(1.,H)*HMIN BX = .FALSE. IER = 1 200 DO 210 J=1,N IJK0 = N + J X0(J) = WK(IJK0) 210 CONTINUE X = XS GO TO 90 C 220 IF (Q.GE.0.03125*R) BE = .FALSE. 230 CONTINUE GO TO (240, 270, 250, 260, 270), SW 240 X = X + H3 GO TO 270 C 250 X = X + P5*H3 GO TO 270 C 260 X = X + P5*H 270 CONTINUE C ************************************ C TEST A POSSIBLE DOUBLING OF THE STEP C ************************************ IF (.NOT.(BE .AND. BH .AND. BR)) GO TO 280 H = H + H BX = .TRUE. 280 BH = .TRUE. GO TO 70 C 290 DO 300 I=1,N X0(I) = ZERO 300 CONTINUE RETURN END FUNCTION DZETA(TA, XA, T, X, XL, TAU, F, FINIT, ALPHA, N, DZE 10 * JCOMP, W, HP, JRAY, NN) C C F -INTERPOLATION ROUTINE. CALLED BY DMRODE AND C U - REQUIRES NO USER INTERFACE OTHER THAN USER C N - SUPPLIED EXTERNAL FUNCTION ALPHA (SEE DMRODE). C C - DZETA CALCULATES THE LAG(OR RETARDATION) VIA C T - ALPHA WITH TA,XA AS INPUT, THEN USES 2-POINT C I - HERMITE INTERPOLATION TO APPROXIMATE SOLUTION C O - AT LAG. IF STEP SIZE IS LARGER THAN LAG,HERMITE C N - INTERPOLANT CANNOT BE OBTAINED, SO DZETA CALLS C - A TRUE ONESTEP FDE SOLVER. SINCE DMRODE CALLS C - DZETA THE PARAMETER DESCRIPTION CAN BE FOUND C - IN DMRODE DIMENSION W(1), XA(1) DIMENSION T(1), X(N,1), XL(N,1) EXTERNAL F, FINIT, ALPHA DATA J /1/ IF ((J.GE.JRAY) .AND. (J.NE.1)) J = JRAY - 1 IF (J.LE.0) J = 1 C ************************************ C CALCULATE LAG C ************************************ B = ALPHA(XA,TA,JCOMP) IF (B.EQ.TA) DZETA = XA(JCOMP) IF (B.EQ.TA) RETURN IF (((HP.LT.0.) .AND. (B.GE.TAU)) .OR. ((HP.GE.0.) .AND. * (B.LE.TAU))) GO TO 100 C ************************************ C SEARCH T ARRAY IN ORDER TO FIND C INTERVAL CONTAINING ALPHA(XA,TA,I).. C IF NO SUCH INTERVAL IS FOUND CALL C START OR IF ARRAY EXCEEDS C DIMENSION IQ OF DMRODE, STOP C ************************************ DX = B - T(J) IF (HP.LT.0.) DX = -DX IF (DX) 10, 40, 30 10 IF (J.EQ.1) GO TO 40 J = J - 1 DX = B - T(J) IF (HP.LT.0.) DX = -DX IF (DX) 10, 40, 40 20 J = J + 1 30 IF (J.EQ.NN) GO TO 50 IF (J+1.GT.JRAY) GO TO 90 DDX = B - T(J+1) IF (HP.LT.0.) DDX = -DDX IF (DDX) 40, 20, 20 40 IF (J.EQ.NN) GO TO 50 GO TO 60 C 50 STOP 60 CONTINUE C ************************************ C 2-POINT HERMITE INTERPOLATION C ************************************ IF (T(J+1).GE.TA) GO TO 90 DX = T(J+1) - T(J) HH = B - T(J) DO 70 I=1,N II = I + N W(I) = X(I,J) W(II) = XL(I,J) 70 CONTINUE AC = F(T(J),W,W(N+1),JCOMP) DO 80 I=1,N II = I + N W(I) = X(I,J+1) W(II) = XL(I,J+1) 80 CONTINUE AF2P = F(T(J+1),W,W(N+1),JCOMP) DIVDF1 = (X(JCOMP,J+1)-X(JCOMP,J))/DX DIVDF3 = AC + AF2P - 2.*DIVDF1 C3 = (DIVDF1-AC-DIVDF3)/DX C4 = DIVDF3/DX**2 DZETA = X(JCOMP,J) + HH*(AC+HH*(C3+HH*C4)) GO TO 110 C 90 CONTINUE J = JRAY CALL START(TA, XA, T, X, XL, F, FINIT, ALPHA, N, JCOMP, W, J, * ANS, HP) DZETA = ANS RETURN C 100 DZETA = FINIT(B,JCOMP) 110 CONTINUE RETURN C C END SUBROUTINE START(TA, XA, T, X, XL, F, FINT, ALPHA, N, JCOMP, STA 10 * W, JEND, ANS, H1) C F -START IS A SUBROUTINE THAT CAN APPROXIMATE THE C U - SOLUTION AT A PAST POINT GIVEN ONLY THE SOLUTION C N - AT THE MOST RECENT MESH POINT AND THE FUNCTIONS C C - F AND ALPHA, AS DESCRIBED BY DMRODE. IT IS A C T - RUNGE-KUTTA TYPE ALGORITHM MODIFIED FOR RODES C I - AND IS GLOBALLY 3RD ORDER, BUT PREDICTS 4TH C O - ORDER APPROX. TO BE USED IN EVALUATION OF F C N - IN MERSON METHOD OF DMRODE. START, DZETA, AND C - DMRODE COMBINE TO GIVE 4TH ORDER SIMPLIFIED C - PREDICTOR-CORRECTOR TWO STEP ALGORITHM, WHERE C - START IS A TRUE ONE-STEP PREDICTOR USED ONLY C - WHEN STEP SIZE EXCEEDS THE LAG DIMENSION T(1), X(N,1), XL(N,1), XA(1), W(1) I1 = N I2 = I1 + N I3 = I2 + N I4 = I3 + N I5 = I4 + N I6 = I5 + N I7 = I6 + N I8 = I7 + N I9 = I8 + N JRAY = JEND C ************************************ C RECALCULATE LAG C ************************************ B = ALPHA(XA,TA,JCOMP) 10 TEST = B - T(JRAY) IF (H1.LT.0.) TEST = -TEST IF (TEST.GT.0.) GO TO 20 JRAY = JRAY - 1 GO TO 10 C ************************************ C CONSTRUCT EULER-TYPE PREDICTOR C CALCULATE APPROX. SOLN. AT TA+H C SEARCH FOR LAG CORRESPONDING TO C X(TA+H) JUST CALCULATED. C APPROXIMATE X AT LAG VIA EULER C ************************************ 20 CONTINUE H = H1 ISWITC = 1 IF (JRAY.NE.JEND) H = T(JRAY+1) - T(JRAY) DO 30 I=1,N L5 = I5 + I L6 = I6 + I W(L5) = X(I,JRAY) W(L6) = XL(I,JRAY) 30 CONTINUE DO 40 I=1,N L1 = I1 + I W(L1) = F(T(JRAY),W(I5+1),W(I6+1),I) 40 CONTINUE DO 50 I=1,N L5 = I5 + I L1 = I1 + I W(L5) = X(I,JRAY) + H*W(L1) C ************************************ C IN EACH COMPONENT CALCULATE THE C CORRESPONDING LAG,W(I), FOR THE C EULER PREDICTED VALUE W(L5), AND C EXECUTE SEARCH. C ************************************ 50 CONTINUE DO 90 I=1,N W(I) = ALPHA(W(I5+1),T(JRAY)+H,I) J = JRAY 60 TEST = W(I) - T(J) IF (H1.LT.0.) TEST = -TEST IF (TEST.GE.0.) GO TO 70 J = J - 1 IF (J.EQ.0) GO TO 70 GO TO 60 C 70 CONTINUE L9 = I9 + I W(L9) = J IF (J.EQ.0) GO TO 90 DO 80 II=1,N L3 = I3 + II L4 = I4 + II W(L3) = X(II,J) W(L4) = XL(II,J) 80 CONTINUE 90 CONTINUE DO 100 I=1,N L9 = I9 + I J = W(L9) H = H1 L6 = I6 + I IF (J.EQ.0) W(L6) = FINT(W(I),I) IF (J.EQ.0) GO TO 100 IF (J.NE.JEND) H = T(J+1) - T(J) C ************************************ C W(L6) IS THE SOLUTION AT ALPHA(W(I),T,I) C ************************************ W(L6) = X(I,J) + (W(I)-T(J))*F(T(J),W(I3+1),W(I4+1),I) 100 CONTINUE C ************************************ C REPEAT ENTIRE PROCEDURE ABOVE C USING EULER VALUE TO CALCULATE C HEUN-TYPE 2ND ORDER PREDICTOR C ************************************ DO 110 I=1,N L2 = I2 + I W(L2) = F(T(JRAY)+H,W(I5+1),W(I6+1),I) 110 CONTINUE DO 120 I=1,N L1 = I1 + I L2 = I2 + I L3 = I3 + I L7 = I7 + I C ************************************ C FOR EACH COMPONENT ESTIMATE THE C SOLUTION AT T+H AND T+H/2 AND STORE C ANSWERS IN W(L3) AND W(L7) RESP. C ************************************ W(L3) = X(I,JRAY) + H*W(L1) + .5*H*(W(L2)-W(L1)) W(L7) = X(I,JRAY) + .5*H*W(L1) + .125*H*(W(L2)-W(L1)) C ************************************ C IN STATEMENTS 120 THRU 280 THE COR- C RESPONDING LAGS AND THE SOLUTION AT C THESE LAGS ARE COMPUTED TO 2ND ORDER C BY HEUN METHOD AND STORED IN W(L4) AND C W(L8) RESPECTIVELY C ************************************ 120 CONTINUE DO 130 I=1,N W(I) = ALPHA(W(I3+1),T(JRAY)+H,I) 130 CONTINUE 140 DO 170 I=1,N J = JRAY 150 TEST = W(I) - T(J) IF (H1.LT.0.) TEST = -TEST IF (TEST.GE.0.) GO TO 160 J = J - 1 IF (J.EQ.0) GO TO 160 GO TO 150 C 160 CONTINUE L9 = I9 + I W(L9) = J 170 CONTINUE DO 290 I=1,N L9 = I9 + I J = W(L9) H = H1 IF (J.EQ.0) GO TO 250 IF (J.NE.JEND) H = T(J+1) - T(J) IF (J.EQ.JEND) GO TO 220 DO 180 II=1,N L5 = I5 + II L6 = I6 + II W(L5) = X(II,J) W(L6) = XL(II,J) 180 CONTINUE DO 190 II=1,N L1 = I1 + II W(L1) = F(T(J),W(I5+1),W(I6+1),II) 190 CONTINUE DO 200 II=1,N L5 = I5 + II L6 = I6 + II W(L5) = X(II,J+1) W(L6) = XL(II,J+1) 200 CONTINUE DO 210 II=1,N L2 = I2 + II W(L2) = F(T(J+1),W(I5+1),W(I6+1),II) 210 CONTINUE 220 GO TO (230, 240), ISWITC 230 L1 = I1 + I L2 = I2 + I L4 = I4 + I W(L4) = X(I,J) + (W(I)-T(J))*W(L1) + .5*(W(I)-T(J))**2* * (W(L2)-W(L1))/H GO TO 260 C 240 L1 = I1 + I L2 = I2 + I L8 = I8 + I W(L8) = X(I,J) + (W(I)-T(J))*W(L1) + .5*(W(I)-T(J))**2* * (W(L2)-W(L1))/H 250 CONTINUE L4 = I4 + I IF (ISWITC.EQ.1) W(L4) = FINT(W(I),I) L8 = I8 + I IF (ISWITC.EQ.2) W(L8) = FINT(W(I),I) GO TO 260 C 260 GO TO (270, 290), ISWITC 270 DO 280 II=1,N W(II) = ALPHA(W(I7+1),T(JRAY)+.5*H,II) 280 CONTINUE ISWITC = ISWITC + 1 GO TO 140 C C ************************************ C USE HEUN PREDICTOR FOR 4TH ORDER C APPROX. SOLUTION. C ************************************ 290 CONTINUE H = H1 IF (JRAY.NE.JEND) H = T(JRAY+1) - T(JRAY) D = F(T(JRAY)+H,W(I3+1),W(I4+1),JCOMP) C = F(T(JRAY)+.5*H,W(I7+1),W(I8+1),JCOMP) SS = B - T(JRAY) R = SS/H DO 300 I=1,N L5 = I5 + I L6 = I6 + I W(L5) = X(I,JRAY) W(L6) = XL(I,JRAY) 300 CONTINUE I = JCOMP L1 = I1 + I W(L1) = F(T(JRAY),W(I5+1),W(I6+1),JCOMP) C ************************************ C RETURN 4TH ORDER APROXIMATION FOR C USE IN F IN MERSON METHOD C ************************************ ANS = X(I,JRAY) + SS*(W(L1)+.5*R*(-3.*W(L1)+4.*C-D)+R*R/3.* * (2.*W(L1)-4.*C+2.*D)) RETURN C END