SUBROUTINE BSPVN(T,N,K,INDEX,X,ILEFT,VNIKX,WORK,IWORK) C***AUTHOR AMOS, D. E., (SNLA) C***PURPOSE Calculates the value of all (possibly) nonzero basis C functions at X. C***DESCRIPTION C C Written by Carl de Boor and modified by D. E. Amos C C Reference C SIAM J. Numerical Analysis, 14, No. 3, June, 1977, pp.441-472. C C Abstract C BSPVN is the BSPLVN routine of the reference. C C BSPVN calculates the value of all (possibly) nonzero basis C functions at X, where C T(K) .LE. X .LE. T(N+1) and J=IWORK is set inside the routine C on the first call when INDEX=1. ILEFT is such that T(ILEFT) C .LE. X .LT. T(ILEFT+1). A call to INTRV(T,N+K,X,ILO,ILEFT, C MFLAG) produces the proper ILEFT. C C Description of Arguments C Input C T - knot vector of length N+K C N - number of B-spline basis functions C = sum of knot multiplicities-K C K - highest possible order C INDEX - INDEX = 1 gives basis functions of order K C = 2 denotes previous entry with WORK, IWORK C values saved for subsequent calls to C BSPVN. C X - argument of basis functions, C T(K) .LE. X .LE. T(N+1) C C Output C VNIKX - vector of length K for spline values. C ILEFT - largest integer such that C T(ILEFT) .LE. X .LT. T(ILEFT+1) C WORK - a work vector of length 2*K C IWORK - a work parameter. Both WORK and IWORK contain C information necessary to continue for INDEX = 2. C When INDEX = 1 exclusively, these are scratch C variables and can be used for other purposes. C C***REFERENCES C. DE BOOR, *PACKAGE FOR CALCULATING WITH B-SPLINES*, C SIAM JOURNAL ON NUMERICAL ANALYSIS, VOLUME 14, NO. 3, C JUNE 1977, PP. 441-472. C***END PROLOGUE BSPVN C C INTEGER ILEFT, IMJP1, INDEX, IPJ, IWORK, JP1, JP1ML, K, L REAL T, VM, VMPREV, VNIKX, WORK, X DIMENSION T(N+K), VNIKX(K), WORK(2*K) DATA ILO/1/ C CONTENT OF J, DELTAM, DELTAP IS EXPECTED UNCHANGED BETWEEN CALLS. C WORK(I) = DELTAP(I), WORK(K+I) = DELTAM(I), I = 1,K C***FIRST EXECUTABLE STATEMENT BSPVN CALL INTRV (T, N+K, X, ILO, ILEFT, MFLAG) IF(K.LT.1) GO TO 900 IF(INDEX.LT.1 .OR. INDEX.GT.2) GO TO 900 IF(X.LT.T(ILEFT) .OR. X.GT.T(ILEFT+1)) GO TO 900 JHIGH= K GO TO (10, 20), INDEX 10 IWORK = 1 VNIKX(1) = 1.0E0 IF (IWORK.GE.JHIGH) GO TO 40 C 20 IPJ = ILEFT + IWORK WORK(IWORK) = T(IPJ) - X IMJP1 = ILEFT - IWORK + 1 WORK(K+IWORK) = X - T(IMJP1) VMPREV = 0.0E0 JP1 = IWORK + 1 DO 30 L=1,IWORK JP1ML = JP1 - L VM = VNIKX(L)/(WORK(L)+WORK(K+JP1ML)) VNIKX(L) = VM*WORK(L) + VMPREV VMPREV = VM*WORK(K+JP1ML) 30 CONTINUE VNIKX(JP1) = VMPREV IWORK = JP1 IF (IWORK.LT.JHIGH) GO TO 20 C 40 RETURN 900 WRITE (*,*) 'PREMATURE TERMINATION IN BSPVD' STOP END FUNCTION BVALU(T,A,N,K,IDERIV,X,INBV,WORK) C Evaluates the B-representation of a B-spline at X for the C function value or any of its derivatives. C Written by Carl de Boor and modified by D. E. Amos C C Reference C SIAM J. Numerical Analysis, 14, No. 3, June, 1977, pp.441-472. C C Abstract C BVALU is the BVALUE function of the reference. C C BVALU evaluates the B-representation (T,A,N,K) of a B-spline C at X for the function value on IDERIV = 0 or any of its C derivatives on IDERIV = 1,2,...,K-1. Right limiting values C (right derivatives) are returned except at the right end C point X=T(N+1) where left limiting values are computed. The C spline is defined on T(K) .LE. X .LE. T(N+1). BVALU returns C a fatal error message when X is outside of this interval. C C To compute left derivatives or left limiting values at a C knot T(I), replace N by I-1 and set X=T(I), I=K+1,N+1. C C BVALU calls INTRV C C Description of Arguments C Input C T - knot vector of length N+K C A - B-spline coefficient vector of length N C N - number of B-spline coefficients C N = sum of knot multiplicities-K C K - order of the B-spline, K .GE. 1 C IDERIV - order of the derivative, 0 .LE. IDERIV .LE. K-1 C IDERIV=0 returns the B-spline value C X - argument, T(K) .LE. X .LE. T(N+1) C INBV - an initialization parameter which must be set C to 1 the first time BVALU is called. C C Output C INBV - INBV contains information for efficient process- C ing after the initial call and INBV must not C be changed by the user. Distinct splines require C distinct INBV parameters. C WORK - work vector of length 3*K. C BVALU - value of the IDERIV-th derivative at X C Reference C. DE BOOR, *PACKAGE FOR CALCULATING WITH B-SPLINES*, C SIAM JOURNAL ON NUMERICAL ANALYSIS, VOLUME 14, NO. 3, C JUNE 1977, PP. 441-472. INTEGER I,IDERIV,IDERP1,IHI,IHMKMJ,ILO,IMK,IMKPJ, INBV, IPJ, 1 IP1, IP1MJ, J, JJ, J1, J2, K, KMIDER, KMJ, KM1, KPK, MFLAG, N REAL A, FKMJ, T, WORK, X C DIMENSION T(N+K), WORK(3*K) DIMENSION T(N+K), A(N), WORK(100) BVALU = 0.0E0 IF(K.LT.1) GO TO 102 IF(N.LT.K) GO TO 101 IF(IDERIV.LT.0 .OR. IDERIV.GE.K) GO TO 110 KMIDER = K - IDERIV C C *** FIND *I* IN (K,N) SUCH THAT T(I) .LE. X .LT. T(I+1) C (OR, .LE. T(I+1) IF T(I) .LT. T(I+1) = T(N+1)). KM1 = K - 1 CALL INTRV(T, N+1, X, INBV, I, MFLAG) IF (X.LT.T(K)) GO TO 120 IF (MFLAG.EQ.0) GO TO 20 IF (X.GT.T(I)) GO TO 130 10 IF (I.EQ.K) GO TO 140 I = I - 1 IF (X.EQ.T(I)) GO TO 10 C C *** DIFFERENCE THE COEFFICIENTS *IDERIV* TIMES C WORK(I) = AJ(I), WORK(K+I) = DP(I), WORK(K+K+I) = DM(I), I=1.K C 20 IMK = I - K DO 30 J=1,K IMKPJ = IMK + J WORK(J) = A(IMKPJ) 30 CONTINUE IF (IDERIV.EQ.0) GO TO 60 DO 50 J=1,IDERIV KMJ = K - J FKMJ = FLOAT(KMJ) DO 40 JJ=1,KMJ IHI = I + JJ IHMKMJ = IHI - KMJ WORK(JJ) = (WORK(JJ+1)-WORK(JJ))/(T(IHI)-T(IHMKMJ))*FKMJ 40 CONTINUE 50 CONTINUE C C *** COMPUTE VALUE AT *X* IN (T(I),(T(I+1)) OF IDERIV-TH DERIVATIVE, C GIVEN ITS RELEVANT B-SPLINE COEFF. IN AJ(1),...,AJ(K-IDERIV). 60 IF (IDERIV.EQ.KM1) GO TO 100 IP1 = I + 1 KPK = K + K J1 = K + 1 J2 = KPK + 1 DO 70 J=1,KMIDER IPJ = I + J WORK(J1) = T(IPJ) - X IP1MJ = IP1 - J WORK(J2) = X - T(IP1MJ) J1 = J1 + 1 J2 = J2 + 1 70 CONTINUE IDERP1 = IDERIV + 1 DO 90 J=IDERP1,KM1 KMJ = K - J ILO = KMJ DO 80 JJ=1,KMJ WORK(JJ) = (WORK(JJ+1)*WORK(KPK+ILO)+WORK(JJ) 1 *WORK(K+JJ))/(WORK(KPK+ILO)+WORK(K+JJ)) ILO = ILO - 1 80 CONTINUE 90 CONTINUE 100 BVALU = WORK(1) RETURN C C 101 CONTINUE WRITE (*,*) ' BVALU, N DOES NOT SATISFY N.GE.K' RETURN 102 CONTINUE WRITE (*,*) ' BVALU, K DOES NOT SATISFY K.GE.1' RETURN 110 CONTINUE WRITE (*,*) ' BVALU, IDERIV DOES NOT SATISFY 0.LE.IDERIV.LT.K' RETURN 120 CONTINUE WRITE (*,*) ' BVALU, X IS N0T GREATER THAN OR EQUAL TO T(K)' RETURN 130 CONTINUE WRITE (*,*) ' BVALU, X IS NOT LESS THAN OR EQUAL TO T(N+1)' RETURN 140 CONTINUE WRITE (*,*) 'BVALU, LEFT LIMITING VALUE CANN0T BE OBTAINED' RETURN END SUBROUTINE RUNTOCH(RUN,CHARRUN) C Convert run# to a character string C Input: Run # (integer) C Output: Character string (CHARACTER*2) CHARACTER*2 CHARRUN INTEGER CNOUGHT, RUN, TENS, ONES C CNOUGHT is the ASCII number that precedes that for the character 1 PARAMETER (CNOUGHT= 48) TENS= RUN/10 ONES= RUN - 10*TENS CHARRUN= CHAR(TENS+CNOUGHT)//CHAR(ONES+CNOUGHT) RETURN END FUNCTION FACT(N) C Factorial function INTEGER FACT FACT= 1 DO 10 I= 2,N FACT= FACT*I 10 CONTINUE RETURN END C SUBROUTINE LFNC (INIT,M,L,THETA,PB,W) C C DIMENSION OF W(5*L), PB(L+1) C ARGUMENTS C C PURPOSE ROUTINES LFNA AND LFNC BOTH CALCULATE SINGLE C PRECISION, NORMALIZED ASSOCIATED LEGENDRE C FUNCTIONS PBAR(N,M,THETA) FOR N=M,...,L-1 C SUBROUTINE LFNC DIFFERS FROM LFNA IN THAT THETA C IS USER SPECIFIED RATHER THAN EQUALLY SPACED C C USAGE CALL LFNC(INIT,M,L,THETA,PB,W) C C ARGUMENTS C C ON INPUT INIT C = 0 INITIALIZATION ONLY C = 1 COMPUTE PBAR(N,M,THETA) C C LFNC CALL WITH INIT = 0 INITIALIZES ARRAY W; C NO VALUES OF PBAR(N,M,THETA) ARE COMPUTED. C INIT=0 SHOULD BE USED ON THE FIRST CALL, OR C IF L, M, OR W VALUES DIFFER FROM THOSE IN THE C PREVIOUS CALL. C C M C NONNEGATIVE INTEGER, LESS THAN L, SPECIFYING C THE ORDER OF PBAR(N,M,THETA), N=M,...,L-1. C C L C THE MAXIMUM VALUE OF N+1, WHERE N IS DEFINED C IN THE 'PURPOSE' SECTION ABOVE. C C THETA C COLATITUDE, IN RADIANS, FOR PBAR(N,M,THETA). C C W C SINGLE PRECISION ARRAY REQUIRING SPECIFIC C VALUES ONLY WHEN THE VALUE OF THE INPUT C PARAMETERS L AND M ARE THE SAME IN THE C PREVIOUS CALL. IN THIS CASE THE VALUES OF C W MUST BE AS RETURNED IN THE PREVIOUS CALL. C C ON OUTPUT PB C SINGLE PRECISION ARRAY STORING PBAR(N,M,THETA) C IN PB(N+1) , N=M,...,L-1 C C W C SINGLE PRECISION ARRAY CONTAINING VALUES C WHICH MUST NOT BE DESTROYED UNTIL THE NEXT C CALL DIFFERS IN VALUE OF INPUT PARAMETERS C L OR M C C SPECIAL CONDITIONS ROUTINE LFNC IS MOST EFFICIENT IF THE VALUES C OF THE INPUT PARAMETERS L AND M REMAIN FIXED ON C CONSECUTIVE CALLS C C ALGORITHM ROUTINE LFNC SOLVES A TRIDIAGONAL SYSTEM OF C EQUATIONS TO OBTAIN VALUES OF PBAR(N,M,THETA) SUBROUTINE LFNC (INIT,M,L,THETA,PB,W) DIMENSION PB(1) ,W(1) IW1 = L+1 IW2 = IW1+L IW3 = IW2+L IW4 = IW3+L CALL LFNC1(INIT,M,L,THETA,PB,W,W(IW1),W(IW2),W(IW3),W(IW4)) RETURN END SUBROUTINE LFNC1(INIT,M,L,THETA,PB,PZ,P1,A,B,W) C Primary routine for calculating Legendre polynomials DIMENSION PB(1) ,PZ(1) ,P1(1) ,A(1) , 1 B(1) ,W(1) DO 10 N=1,L PB(N) = 0. 10 CONTINUE IF( INIT.NE.0) GO TO 70 LM1 = L-1 LM2 = L-2 LMM = L-M PI = 3.14159265358979 CALL ALFK(LM2,M,P1) CALL ALFK(LM1,M,PZ) K = 0 N = L-1 FNMM = N-M+1 FNPN = N+N+1 FNPM = N+M+1 50 N = N-1 IF(N.LT.M) RETURN 60 FNMM = FNMM-1. FNPN = FNPN-2. FNPM = FNPM-1. K = K+1 B(K) = SQRT(FNMM*FNPM/(FNPN*(FNPN+2.))) GO TO 50 70 CALL LFPT(LM1,M,THETA,PZ,PZI) PB(L) = PZI IF(L.EQ.1) RETURN CALL LFPT(LM2,M,THETA,P1,P1I) IF(LMM-2) 140,80,90 80 PB(L-1) = P1I GO TO 140 90 COST = COS(THETA) DO 100 K=1,LMM PB(K) = -COST A(K) = B(K) 100 CONTINUE IF (ABS(PZI) .LT. ABS(P1I)) GO TO 110 PB(1) = PZI R = -A(1)*PB(1) CALL TRIH (LMM-1,A,PB(2),A(2),R) GO TO 120 110 PB(1) = PZI PB(2) = P1I IF(LMM.EQ.1) GO TO 120 R = -A(2)*PB(2) CALL TRIH (LMM-2,A(2),PB(3),A(3),R) 120 NDO = (L+1)/2 DO 130 N=1,NDO N1 = L-N PHOLD = PB(N1+1) PB(N1+1) = PB(N) PB(N) = PHOLD 130 CONTINUE 140 RETURN END SUBROUTINE ALFK (N,M,CP) C For a given M and N, the unnormalized Legendre polynomial are produced. DIMENSION CP(1) CP(1) = 0. MA = IABS(M) IF(MA .GT. N) RETURN IF(N-1) 2,3,5 2 CP(1) = SQRT(2.) RETURN 3 IF(MA .NE. 0) GO TO 4 CP(1) = SQRT(1.5) RETURN 4 CP(1) = SQRT(.75) IF(M .EQ. -1) CP(1) = -CP(1) RETURN 5 IF(MOD(N+MA,2) .NE. 0) GO TO 10 NMMS2 = (N-MA)/2 FNUM = N+MA+1 FNMH = N-MA+1 PM1 = 1. GO TO 15 10 NMMS2 = (N-MA-1)/2 FNUM = N+MA+2 FNMH = N-MA+2 PM1 = -1. 15 T1 = 1. T2 = 1. IF(NMMS2 .LT. 1) GO TO 20 FDEN = 2. DO 18 I=1,NMMS2 T1 = FNUM*T1/FDEN FNUM = FNUM+2. FDEN = FDEN+2. 18 CONTINUE 20 IF(MA .EQ. 0) GO TO 26 DO 25 I=1,MA T2 = FNMH*T2/(FNMH+PM1) FNMH = FNMH+2. 25 CONTINUE 26 IF(MOD(MA/2,2) .NE. 0) T1 = -T1 CP2 = T1*SQRT((N+.5)*T2)/(2.**(N-1)) FNNP1 = N*(N+1) FNMSQ = FNNP1-2.*MA*MA L = (N+1)/2 IF(MOD(N,2) .EQ. 0 .AND. MOD(MA,2) .EQ. 0) L = L+1 CP(L) = CP2 IF(M .GE. 0) GO TO 29 IF(MOD(MA,2) .NE. 0) CP(L) = -CP(L) 29 IF(L .LE. 1) RETURN FK = N A1 = (FK-2.)*(FK-1.)-FNNP1 B1 = 2.*(FK*FK-FNMSQ) CP(L-1) = B1*CP(L)/A1 30 L = L-1 IF(L .LE. 1) RETURN FK = FK-2. A1 = (FK-2.)*(FK-1.)-FNNP1 B1 = -2.*(FK*FK-FNMSQ) C1 = (FK+1.)*(FK+2.)-FNNP1 CP(L-1) = -(B1*CP(L)+C1*CP(L+1))/A1 GO TO 30 END SUBROUTINE LFPT (N,M,THETA,CP,PB) C Uses recursion to produce the Legendre polynomial given their values C for N=0 and N=1 DIMENSION CP(1) PB = 0. MA = IABS(M) IF(MA .GT. N) RETURN IF (N) 10, 10, 30 10 IF (MA) 20, 20, 30 20 PB= SQRT(.5) GO TO 140 30 NP1 = N+1 NMOD = MOD(N,2) MMOD = MOD(MA,2) IF (NMOD) 40, 40, 90 40 IF (MMOD) 50, 50, 70 50 KDO = N/2+1 CDT = COS(THETA+THETA) SDT = SIN(THETA+THETA) CT = 1. ST = 0. SUM = .5*CP(1) DO 60 KP1=2,KDO CTH = CDT*CT-SDT*ST ST = SDT*CT+CDT*ST CT = CTH SUM = SUM+CP(KP1)*CT 60 CONTINUE PB= SUM GO TO 140 70 KDO = N/2 CDT = COS(THETA+THETA) SDT = SIN(THETA+THETA) CT = 1. ST = 0. SUM = 0. DO 80 K=1,KDO CTH = CDT*CT-SDT*ST ST = SDT*CT+CDT*ST CT = CTH SUM = SUM+CP(K)*ST 80 CONTINUE PB= SUM GO TO 140 90 KDO = (N+1)/2 IF (MMOD) 100,100,120 100 CDT = COS(THETA+THETA) SDT = SIN(THETA+THETA) CT = COS(THETA) ST = -SIN(THETA) SUM = 0. DO 110 K=1,KDO CTH = CDT*CT-SDT*ST ST = SDT*CT+CDT*ST CT = CTH SUM = SUM+CP(K)*CT 110 CONTINUE PB= SUM GO TO 140 120 CDT = COS(THETA+THETA) SDT = SIN(THETA+THETA) CT = COS(THETA) ST = -SIN(THETA) SUM = 0. DO 130 K=1,KDO CTH = CDT*CT-SDT*ST ST = SDT*CT+CDT*ST CT = CTH SUM = SUM+CP(K)*ST 130 CONTINUE PB= SUM 140 RETURN END SUBROUTINE TRIH (N,A,B,C,R) C A tri-diagonal matrix solver DIMENSION A(1) ,B(1) ,C(1) IF (N) 120,120, 10 10 IF (N-2) 20, 30, 40 20 B(1) = R/B(1) RETURN 30 QIH = A(2) BIH = B(2) GO TO 70 40 QIH = A(N) BIH = B(N) DO 60 IDO=3,N I = N-IDO+2 IF (ABS(BIH) .LT. ABS(C(I))) GO TO 50 RATIO = C(I)/BIH C(I) = 0. B(I+1) = QIH/BIH BIH = B(I)-RATIO*QIH QIH = A(I) GO TO 60 50 B(I+1) = B(I)/C(I) C(I) = A(I)/C(I) BIH1 = QIH-BIH*B(I+1) QIH = -BIH*C(I) BIH = BIH1 60 CONTINUE 70 IF (ABS(BIH) .LT. ABS(C(1))) GO TO 80 Q2 = QIH/BIH BIH = B(1)-C(1)/BIH*QIH B(1) = R/BIH B(2) = -Q2*B(1) GO TO 90 80 RATIO = BIH/C(1) BIH = QIH-RATIO*B(1) RIH = -RATIO*R B1 = RIH/BIH B(2) = (R-B(1)*B1)/C(1) B(1) = B1 90 IF (N-3) 120,100,100 100 DO 110 I=3,N B(I) = -B(I)*B(I-1)-C(I-1)*B(I-2) 110 CONTINUE 120 RETURN END c********************************************************************* SUBROUTINE INTRV(XT,LXT,X,ILO,ILEFT,MFLAG) C***PURPOSE Computes the largest integer ILEFT in 1.LE.ILEFT.LE.LXT C such that XT(ILEFT).LE.X where XT(*) is a subdivision C of the X interval. C***DESCRIPTION C C Written by Carl de Boor and modified by D. E. Amos C C Reference C SIAM J. Numerical Analysis, 14, No. 3, June 1977, pp. 441-472. C C Abstract C INTRV is the INTERV routine of the reference. C C INTRV computes the largest integer ILEFT in 1 .LE. ILEFT .LE. C LXT such that XT(ILEFT) .LE. X where XT(*) is a subdivision of C the X interval. Precisely, C C X .LT. XT(1) 1 -1 C if XT(I) .LE. X .LT. XT(I+1) then ILEFT=I , MFLAG=0 C XT(LXT) .LE. X LXT 1, C C That is, when multiplicities are present in the break point C to the left of X, the largest index is taken for ILEFT. C C Description of Arguments C Input C XT - XT is a knot or break point vector of length LXT C LXT - length of the XT vector C X - argument C ILO - an initialization parameter which must be set C to 1 the first time the spline array XT is C processed by INTRV. C C Output C ILO - ILO contains information for efficient process- C ing after the initial call, and ILO must not be C changed by the user. Distinct splines require C distinct ILO parameters. C ILEFT - largest integer satisfying XT(ILEFT) .LE. X C MFLAG - signals when X lies out of bounds C C Error Conditions C None C***REFERENCES C. DE BOOR, *PACKAGE FOR CALCULATING WITH B-SPLINES*, C SIAM JOURNAL ON NUMERICAL ANALYSIS, VOLUME 14, NO. 3, C JUNE 1977, PP. 441-472. C***ROUTINES CALLED (NONE) C***END PROLOGUE INTRV C C INTEGER IHI, ILEFT, ILO, ISTEP, LXT, MFLAG, MIDDLE REAL X, XT DIMENSION XT(LXT) C***FIRST EXECUTABLE STATEMENT INTRV IHI = ILO + 1 IF (IHI.LT.LXT) GO TO 10 IF (X.GE.XT(LXT)) GO TO 110 IF (LXT.LE.1) GO TO 90 ILO = LXT - 1 IHI = LXT C 10 IF (X.GE.XT(IHI)) GO TO 40 IF (X.GE.XT(ILO)) GO TO 100 C C *** NOW X .LT. XT(IHI) . FIND LOWER BOUND ISTEP = 1 20 IHI = ILO ILO = IHI - ISTEP IF (ILO.LE.1) GO TO 30 IF (X.GE.XT(ILO)) GO TO 70 ISTEP = ISTEP*2 GO TO 20 30 ILO = 1 IF (X.LT.XT(1)) GO TO 90 GO TO 70 C *** NOW X .GE. XT(ILO) . FIND UPPER BOUND 40 ISTEP = 1 50 ILO = IHI IHI = ILO + ISTEP IF (IHI.GE.LXT) GO TO 60 IF (X.LT.XT(IHI)) GO TO 70 ISTEP = ISTEP*2 GO TO 50 60 IF (X.GE.XT(LXT)) GO TO 110 IHI = LXT C C *** NOW XT(ILO) .LE. X .LT. XT(IHI) . NARROW THE INTERVAL 70 MIDDLE = (ILO+IHI)/2 IF (MIDDLE.EQ.ILO) GO TO 100 C NOTE. IT IS ASSUMED THAT MIDDLE = ILO IN CASE IHI = ILO+1 IF (X.LT.XT(MIDDLE)) GO TO 80 ILO = MIDDLE GO TO 70 80 IHI = MIDDLE GO TO 70 C *** SET OUTPUT AND RETURN 90 MFLAG = -1 ILEFT = 1 RETURN 100 MFLAG = 0 ILEFT = ILO RETURN 110 MFLAG = 1 ILEFT = LXT RETURN END FUNCTION GETROOT(X1,X2,IER, KNOT, COEFF, KORDER, NCOEFF,ALT,VAR) CD FUNCTION GETROOT CD CD Purpose: CD find a root known to be within a given bound, using combined Newton- CD Raphson and bisection CD X1 Real - none lower bound CD X2 Real - none upper bound CD XACC Real - none tolerence (+/- XACC) CD CD Output CD Parameters Type Units Dimensions Description CD IER Integer - none error code CD = 0, no error CD = 1, not bracketed CD = 2, too many iterations CD GETROOT Real - none root CD Method: CD See "Numerical Recipies", 1986 edition, section 9.4, CD "Netwon-Raphson Method Using Derivatives", page 254. c MAXIT...max iterations c DIMENSION COEFF(NCOEFF), WORK(100) REAL KNOT(KORDER+NCOEFF), LOGPFAC, LOGALT PARAMETER (MAXIT=100, XACC=1E-5, LOGPFAC= -25.836952) INTEGER VAR IER = 0 XL=X1 XH=X2 GETROOT=.5*(X1+X2) c c old stepsize c DXOLD=ABS(X2-X1) DX=DXOLD C CALL FUNCD(GETROOT,F,DF,A) INBV=1 IF (VAR .EQ. 6) THEN F= BVALU(KNOT,COEFF,NCOEFF,KORDER,0,GETROOT,INBV,WORK) - ALT DF=BVALU(KNOT,COEFF,NCOEFF,KORDER,1,GETROOT,INBV,WORK) ELSE F= BVALU(KNOT,COEFF,NCOEFF,KORDER,0,GETROOT,INBV,WORK) 1 - EXP(LOGPFAC - GETROOT - ALT) DF=BVALU(KNOT,COEFF,NCOEFF,KORDER,1,GETROOT,INBV,WORK) 1 + EXP(LOGPFAC - GETROOT - ALT) ENDIF c DO 11 J=1,MAXIT c c bisect if Newton is out of range or too slow c IF(((GETROOT-XH)*DF-F)*((GETROOT-XL)*DF-F).GE.0. * .OR. ABS(2.*F).GT.ABS(DXOLD*DF) ) THEN DXOLD=DX DX=0.5*(XH-XL) GETROOT=XL+DX c c no change in roots c IF(XL.EQ.GETROOT)RETURN ELSE c c use Newton step c DXOLD=DX DX=F/DF TEMP=GETROOT GETROOT=GETROOT-DX IF(TEMP.EQ.GETROOT)RETURN ENDIF c c convergence ? c IF(ABS(DX).LT.XACC) RETURN c c reevaluate function c INBV=1 IF (VAR .EQ. 6) THEN F= BVALU(KNOT,COEFF,NCOEFF,KORDER,0,GETROOT,INBV,WORK) - ALT DF=BVALU(KNOT,COEFF,NCOEFF,KORDER,1,GETROOT,INBV,WORK) ELSE F= BVALU(KNOT,COEFF,NCOEFF,KORDER,0,GETROOT,INBV,WORK) 1 - EXP(LOGPFAC - GETROOT - ALT) DF=BVALU(KNOT,COEFF,NCOEFF,KORDER,1,GETROOT,INBV,WORK) 1 + EXP(LOGPFAC - GETROOT - ALT) ENDIF C CALL FUNCD(GETROOT,F,DF,A) c c new bracket c IF(F.LT.0.) THEN XL=GETROOT FL=F ELSE XH=GETROOT FH=F ENDIF 11 CONTINUE c c took too long c IER = 2 RETURN END