C Written by Shunrong Zhang C Imported into Madrigal from modelRecovery/fullAnalytic/ C basis.f on 10/21/2005 - Was revision 1.2 C C $Id: basis.f 3728 2011-10-21 20:47:30Z brideout $ C SUBROUTINE POLBAS(N,IX,X,F) C C JMH - 9/88 C C POLBAS COMPUTES THE FIRST N MONOMIALS AT X, EG. 1.0,X,X**2 ETC. C THE RESULTS ARE RETURNED IN F. IX IS NOT USED. SEE FTL1L2 FOR AN C EXAMPLE OF A SUBROUTINE WHICH CALLS BASIS FUNCTION ROUTINES OF C THIS FORM. C C .. Scalar Arguments .. DOUBLE PRECISION X INTEGER IX,N C .. C .. Array Arguments .. DOUBLE PRECISION F(*) C .. C .. Local Scalars .. INTEGER I C .. F(1) = 1.0D0 DO 10 I = 2,N F(I) = X*F(I-1) 10 CONTINUE RETURN END C C SUBROUTINE HRMBAS(N,IX,X,F) C C JMH - 9/88 C C HRMBAS COMPUTES THE FIRST N TERMS OF A HARMONIC EXPANSION AT X, EG C 1.0, COS(X), SIN(X), COS(2*X), SIN(2*X), ... THE RESULTS ARE C RETURNED IN F. IX IS NOT USED. C C .. Scalar Arguments .. DOUBLE PRECISION X INTEGER IX,N C .. C .. Array Arguments .. DOUBLE PRECISION F(*) C .. C .. Local Scalars .. INTEGER I C .. C .. Intrinsic Functions .. INTRINSIC COS,SIN C .. F(1) = 1.0D0 DO 10 I = 2,N,2 F(I) = COS((I/2)*X) 10 CONTINUE DO 20 I = 3,N,2 F(I) = SIN((I/2)*X) 20 CONTINUE RETURN END C SUBROUTINE SPLBAS(N,IX,X,F) C C JMH - 9/88 C JMH - 7/02 - UPDATED C C SPLBAS COMPUTES THE N B-SPLINES OF ORDER K AND KNOT SEQUENCE K AT C X. THE SPLINE AND ITS DERIVATIVES ARE RETURNED IN F. SEE FTL1L2 C FOR AN EXAMPLE OF A SUBROUTINE WHICH CALLS BASIS FUNCTION C ROUTINES OF THIS FORM. COMMON/SPFTCM/ (IN basis.h) MUST BE SET UP C BEFORE SPLBAS IS CALLED. SEE SUBROUTINE KNOTS1. SPLBAS CALLS C SUBROUTINES INTERV AND BSPLVD WHICH ARE FROM CARL DEBOOR'S C B-SPLINE PACKAGE. C INCLUDE 'basis.h' INCLUDE 'basis1.h' C C C .. Scalar Arguments .. DOUBLE PRECISION X INTEGER IX,N C .. C .. Array Arguments .. DOUBLE PRECISION F(NDIM,*) C .. C .. Local Scalars .. INTEGER I,ILEFT,ILFTMK,J,K,MFLAG,MI,NDERIV C .. C .. Local Arrays .. DOUBLE PRECISION VAL(4,4) C .. C .. External Subroutines .. EXTERNAL BSPLVD,INTERV C .. K = 4 NDERIV = 4 DO 20 I = 1,N DO 10 J = 1,K F(I,J) = 0.0D0 10 CONTINUE 20 CONTINUE CALL INTERV(T,N+K,X,ILEFT,MFLAG) IF (ILEFT.GT.N) ILEFT = N ILFTMK = ILEFT - K CALL BSPLVD(T,K,X,ILEFT,VAL,NDERIV) DO 40 MI = 1,K I = ILFTMK + MI DO 30 J = 1,K F(I,J) = VAL(MI,J) 30 CONTINUE 40 CONTINUE RETURN C END C SUBROUTINE SPLBASP(N,IX,X,F) C C JMH - 2/02 C C SPLBASP COMPUTES THE N PERIODIC B-SPLINES OF ORDER K AND KNOT C SEQUENCE K AT X. THE RESULTS OR ONE OF THE FIRST THREE C DERIVATIVES ARE RETURNED IN F. IX SPECIFIES THE ORDER OF THE C DERIVATIVE, I.E IX=1 FOR THE FIRST DERIVATIVE. SEE FTL1L2 FOR AN C EXAMPLE OF A SUBROUTINE WHICH CALLS BASIS FUNCTION ROUTINES OF C THIS FORM. COMMON/SPFTCM/ MUST BE SET UP BEFORE SPLBAS IS C CALLED. SEE SUBROUTINE KNOTSP IN THE MATH LIBRARY FOR AN EXAMPLE C OF A ROUTINE WHICH DOES THIS. SPLBASP CALLS SUBROUTINES INTERV C AND BSPLVD WHICH ARE FROM CARL DEBOOR'S B-SPLINE PACKAGE AND ARE C ALSO IN THE MATH LIBRARY. C INCLUDE 'basis.h' INCLUDE 'basis1.h' C C .. Scalar Arguments .. DOUBLE PRECISION X INTEGER IX,N C .. C .. Array Arguments .. DOUBLE PRECISION F(NDIM,*) C .. C .. Local Scalars .. DOUBLE PRECISION XP INTEGER I,ILEFT,ILFTMK,J,K,MFLAG,MI,NDERIV,NN C .. C .. Local Arrays .. DOUBLE PRECISION VAL(4,4) C .. C .. External Subroutines .. EXTERNAL BSPLVD,INTERV C .. C .. Intrinsic Functions .. INTRINSIC DMOD C .. NN = N + 3 K = 4 NDERIV = 4 DO 20 I = 1,NN DO 10 J = 1,K F(I,J) = 0.0D0 10 CONTINUE 20 CONTINUE XP = DMOD(X,T(NN+1)) CALL INTERV(T,NN+K,X,ILEFT,MFLAG) IF (ILEFT.GT.NN) ILEFT = NN ILFTMK = ILEFT - K CALL BSPLVD(T,K,X,ILEFT,VAL,NDERIV) DO 40 MI = 1,K I = ILFTMK + MI IF (I.GT.N) I = I - N DO 30 J = 1,K F(I,J) = VAL(MI,J) 30 CONTINUE 40 CONTINUE RETURN C END C SUBROUTINE SPLBASP3(N,IX,X,F) C C SRZ - 7/02 C C SPLBASP3 COMPUTES THE 3-DIMENSIONAL TENSOR PRODUCT OF BASIS C FUNCTIONS SPLBASP AND SPLBAS. C INCLUDE 'basis.h' INCLUDE 'basis1.h' INCLUDE 'basis3.h' C C .. Scalar Arguments .. DOUBLE PRECISION X INTEGER IX,N C .. C .. Array Arguments .. DOUBLE PRECISION F(NDIM,*) C .. C .. Subroutine Arguments .. C .. C .. Local Scalars .. INTEGER I,IJ,J,L,LL C .. C .. Local Arrays .. DOUBLE PRECISION F1(NDIM,4),F2(NDIM,4),F3(NDIM,4) C .. C .. External Subroutines .. EXTERNAL SPLBAS,SPLBASP C .. KS = 4 DO 20 I = 1,N DO 10 J = 1,KS F(I,J) = 0.0D0 10 CONTINUE 20 CONTINUE DO 30 I = 1,N1 + 4 T(I) = T1(I) 30 CONTINUE CALL SPLBASP(N1,IX,X1(IX1(IX)),F1) DO 40 I = 1,N2 + 4 T(I) = T2(I) 40 CONTINUE CALL SPLBAS(N2,IX,X2(IX2(IX)),F2) DO I = 1,N3 + 4 T(I) = T3(I) END DO CALL SPLBASP(N3,IX,X3(IX3(IX)),F3) IJ = 0 DO LL = 1,N3 DO 70 I = 1,N1 DO 60 J = 1,N2 IJ = IJ + 1 DO 50 L = 1,KS F(IJ,L) = F1(I,L)*F2(J,L)*F3(LL,L) 50 CONTINUE 60 CONTINUE 70 CONTINUE END DO C C RETURN END c SUBROUTINE SPLBASP35(N,IX,X,F) C C SRZ - 7/02 C C SPLBASP3 COMPUTES THE 3-DIMENSIONAL TENSOR PRODUCT OF BASIS C FUNCTIONS SPLBASP AND SPLBAS. C INCLUDE 'basis.h' INCLUDE 'basis1.h' INCLUDE 'basis3.h' C C .. Scalar Arguments .. DOUBLE PRECISION X INTEGER IX,N C .. C .. Array Arguments .. DOUBLE PRECISION F(NDIM,*) C .. C .. Subroutine Arguments .. C .. C .. Local Scalars .. INTEGER I,IJ,J,L,LL C .. C .. Local Arrays .. DOUBLE PRECISION F1(NDIM,4),F2(NDIM,4),F3(NDIM,4) C .. C .. External Subroutines .. EXTERNAL POLBAS,SPLBAS,SPLBASP C .. KS = 4 DO 20 I = 1,N DO 10 J = 1,KS F(I,J) = 0.0D0 10 CONTINUE 20 CONTINUE DO 30 I = 1,N1 + 4 T(I) = T1(I) 30 CONTINUE CALL SPLBASP(N1,IX,X1(IX1(IX)),F1) DO 40 I = 1,N2 + 4 T(I) = T2(I) 40 CONTINUE CALL SPLBAS(N2,IX,X2(IX2(IX)),F2) DO I = 1,N3 + 4 T(I) = T3(I) END DO CALL POLBAS(N3,IX,X3(IX3(IX)),F3) IJ = 0 DO LL = 1,N3 DO 70 I = 1,N1 DO 60 J = 1,N2 IJ = IJ + 1 DO 50 L = 1,KS F(IJ,L) = F1(I,L)*F2(J,L)*F3(LL,L) 50 CONTINUE 60 CONTINUE 70 CONTINUE END DO C C RETURN END SUBROUTINE SPLBASH3(N,IX,X,F) C C SRZ - 7/02 C C SPLBASP3 COMPUTES THE 3-DIMENSIONAL TENSOR PRODUCT OF BASIS C FUNCTIONS SPLBASP AND HRMBAS. C INCLUDE 'basis.h' INCLUDE 'basis1.h' INCLUDE 'basis3.h' C C .. Scalar Arguments .. DOUBLE PRECISION X INTEGER IX,N C .. C .. Array Arguments .. DOUBLE PRECISION F(NDIM,*) C .. C .. Subroutine Arguments .. C .. C .. Local Scalars .. DOUBLE PRECISION PI INTEGER I,IJ,J,L,LL C .. C .. Local Arrays .. DOUBLE PRECISION F1(NDIM,4),F2(NDIM,4),F3(NDIM,4) C .. C .. External Subroutines .. EXTERNAL HRMBAS,SPLBAS C .. KS = 4 PI = 3.1415926535897d0 DO 20 I = 1,N DO 10 J = 1,KS F(I,J) = 0.0D0 10 CONTINUE 20 CONTINUE DO 30 I = 1,N1 + 4 T(I) = T1(I) 30 CONTINUE T(1) = 2.d0*PI/T1(1) CALL HRMBAS(N1,IX,X1(IX1(IX))*T(1),F1) DO 40 I = 1,N2 + 4 T(I) = T2(I) 40 CONTINUE CALL SPLBAS(N2,IX,X2(IX2(IX)),F2) DO I = 1,N3 + 4 T(I) = T3(I) END DO T(1) = 2.d0*PI/T3(1) CALL HRMBAS(N3,IX,X3(IX3(IX))*T(1),F3) IJ = 0 DO LL = 1,N3 DO 70 I = 1,N1 DO 60 J = 1,N2 IJ = IJ + 1 DO 50 L = 1,KS F(IJ,L) = F1(I,L)*F2(J,L)*F3(LL,L) 50 CONTINUE 60 CONTINUE 70 CONTINUE END DO C C RETURN END SUBROUTINE SPLBASP4(N,IX,X,F) C C SRZ - 7/02 C C SPLBASP4 COMPUTES THE 4-DIMENSIONAL TENSOR PRODUCT OF BASIS C FUNCTIONS HRMBAS, SPLBAS, HRMBAS AND POLBAS C C SRZ -4/05 updated to consider partial derivatives of the second C variable (SPLBAS) C INCLUDE 'basis.h' INCLUDE 'basis1.h' INCLUDE 'basis4.h' C C .. Scalar Arguments .. DOUBLE PRECISION X INTEGER IX,N C .. C .. Array Arguments .. DOUBLE PRECISION F(NDIM,*) C .. C .. Subroutine Arguments .. C .. C .. Local Scalars .. INTEGER I,IJ,J,L,L3,L4 C .. C .. Local Arrays .. DOUBLE PRECISION F1(NDIM,4),F2(NDIM,4),F3(NDIM,4),F4(NDIM,4) C .. C .. External Subroutines .. EXTERNAL HRMBAS,POLBAS,SPLBAS C .. KS = 4 DO 20 I = 1,N DO 10 J = 1,KS F(I,J) = 0.0D0 10 CONTINUE 20 CONTINUE DO 30 I = 1,N1 + 4 T(I) = T1(I) 30 CONTINUE C PERIODIC B-SPLINE C CALL SPLBASP(N1,IX,X1(IX1(IX)),F1) C HARMONIC BASIS CALL HRMBAS(N1,IX,X1(IX1(IX)),F1) DO 40 I = 1,N2 + 4 T(I) = T2(I) 40 CONTINUE C B-SPLINE CALL SPLBAS(N2,IX,X2(IX2(IX)),F2) DO 50 I = 1,N3 + 4 T(I) = T3(I) 50 CONTINUE C HARMONIC BASIS CALL HRMBAS(N3,IX,X3(IX3(IX)),F3) C PERIODIC B-SPLINE C CALL SPLBASP(N3,IX,X3(IX3(IX)),F3) DO 60 I = 1,N4 + 4 T(I) = T4(I) 60 CONTINUE C B-SPLINE IF (T4(1).NE.12345D0) THEN CALL SPLBAS(N4,IX,X4(IX4(IX)),F4) ELSE C POLYNOMIAL CALL POLBAS(N4,IX,X4(IX4(IX)),F4) ENDIF IJ = 0 DO 110 L4 = 1,N4 DO 100 L3 = 1,N3 DO 90 I = 1,N1 DO 80 J = 1,N2 IJ = IJ + 1 DO 70 L = 1,KS F(IJ,L) = F1(I,L)*F2(J,L)*F3(L3,L)*F4(L4,L) C PARTIAL DERIVATIVES FOR VARIABLE 2 IF (L.GT.1) F(IJ,L) = F1(I,1)*F2(J,L)*F3(L3,1)* * F4(L4,1) 70 CONTINUE 80 CONTINUE 90 CONTINUE 100 CONTINUE 110 CONTINUE C C RETURN END C C SUBROUTINE KNOTS2(K,N,X1,X2,T) C C JMH - 9/88 C JMH - 11/01 MODIFIED FROM KNOTS C C F:KNOTS2 - COMPUTES AN ARRAY OF KNOTS FOR A SPLINE FUNCTION. C THE FIRST AND LAST INTER-KNOT SPACINGS ARE HALF AS C LONG AS THE EQUAL INTERIOR KNOT SPACINGS. C C INPUT PARAMETERS: C K - ORDER (=DEGREE+1) OF B-SPLINES (K=4 FOR CUBIC SPLINES) C X1 - LOCATION OF FIRST SPLINE KNOT C X2 - LOCATION OF LAST SPLINE KNOT C C OUTPUT PARAMETERS C T - ARRAY OF KNOTS C N - NUMBER OF B-SPLINES OF ORDER K FOR THE GIVEN KNOT C SEQUENCE (N = NBKPT+K-2). C C OTHER VARIABLES C NBKPT - NUMBER OF BREAKPOINTS C C EXAMPLE (K=4, NBKPT=5, N=7): C C I1=1 I2=8 C | | | | | C | | <-DT-> | | | C | | | | | C 1 5 6 7 8 C 2 9 C 3 10 C 4 11 C C .. Scalar Arguments .. DOUBLE PRECISION X1,X2 INTEGER K,N C .. C .. Array Arguments .. DOUBLE PRECISION T(*) C .. C .. Local Scalars .. DOUBLE PRECISION DT,EPS INTEGER I,I1,I2,NBKPT C .. C .. Data statements .. DATA EPS/1.D-6/ C .. C NBKPT = N - K + 2 DT = (X2-X1)/(NBKPT-2) I1 = K I2 = NBKPT + K - 1 T(I1) = X1 - EPS T(I1+1) = X1 + 0.5D0*DT DO 10 I = I1 + 2,I2 - 1 T(I) = T(I-1) + DT 10 CONTINUE T(I2) = X2 + EPS DO 20 I = I1 - 1,1,-1 T(I) = T(I1) 20 CONTINUE DO 30 I = I2 + 1,I2 + K - 1,1 T(I) = T(I2) 30 CONTINUE N = NBKPT + K - 2 RETURN END C SUBROUTINE KNOTSP(K,N,X1,X2,T) C C JMH - 6/02 C C F:KNOTSP - COMPUTES AN ARRAY OF KNOTS FOR A PERIODIC SPLINE C FUNCTION. C C INPUT PARAMETERS: C K - ORDER (=DEGREE+1) OF B-SPLINES (K=4 FOR CUBIC SPLINES) C X1 - LOCATION OF FIRST SPLINE KNOT C X2 - LOCATION OF LAST SPLINE KNOT C C OUTPUT PARAMETERS C T - ARRAY OF KNOTS C N - NUMBER OF B-SPLINES OF ORDER K FOR THE GIVEN KNOT C SEQUENCE (N = NBKPT+K-2). C C OTHER VARIABLES C NBKPT - NUMBER OF BREAKPOINTS C C EXAMPLE (K=4, NBKPT=5, N=7): C C I1=1 I2=8 C | | | | | C | | <-DT-> | | | C | | | | | C 1 5 6 7 8 C 2 9 C 3 10 C 4 11 C C .. Scalar Arguments .. DOUBLE PRECISION X1,X2 INTEGER K,N C .. C .. Array Arguments .. DOUBLE PRECISION T(*) C .. C .. Local Scalars .. DOUBLE PRECISION DT,PI,TWOPI INTEGER I,I1,I2,J,NBKPT C .. C .. Intrinsic Functions .. INTRINSIC DBLE C .. C .. Data statements .. DATA PI/3.141592653D0/,TWOPI/6.283185307D0/ C .. C PI = PI NBKPT = N - K + 5 DT = TWOPI/DBLE(NBKPT-1) I1 = K I2 = NBKPT + K - 1 T(I1) = 0.D0 DO 10 I = I1 + 1,I2 T(I) = T(I-1) + DT 10 CONTINUE T(I2) = TWOPI DO 20 I = I1 - 1,1,-1 J = I1 - I T(I) = T(I1) - (T(I2)-T(I2-J)) 20 CONTINUE DO 30 I = I2 + 1,I2 + K - 1,1 J = I - I2 T(I) = T(I2) + (T(I1+J)-T(I1)) 30 CONTINUE C DO 40 I = 1,NBKPT + 2*(K-1) C WRITE (6,FMT='('' I,T(I) ='',I4,F10.3)') I,T(I) C 40 CONTINUE RETURN END C SUBROUTINE INTERV(XT,LXT,X,ILEFT,MFLAG) C OMPUTES LARGEST ILEFT IN (1,LXT) SUCH THAT XT(ILEFT) .LE. X C .. Scalar Arguments .. DOUBLE PRECISION X INTEGER ILEFT,LXT,MFLAG C .. C .. Array Arguments .. DOUBLE PRECISION XT(LXT) C .. C .. Local Scalars .. INTEGER IHI,ILO,ISTEP,MIDDLE C .. C .. Data statements .. DATA ILO/1/ C .. IHI = ILO + 1 IF (IHI.LT.LXT) GO TO 10 IF (X.GE.XT(LXT)) GO TO 120 IF (LXT.LE.1) GO TO 100 ILO = LXT - 1 IHI = LXT C 10 IF (X.GE.XT(IHI)) GO TO 50 IF (X.GE.XT(ILO)) GO TO 110 C C **** NOW X .LT. XT(IHI) . FIND LOWER BOUND 20 ISTEP = 1 30 IHI = ILO ILO = IHI - ISTEP IF (ILO.LE.1) GO TO 40 IF (X.GE.XT(ILO)) GO TO 80 ISTEP = ISTEP*2 GO TO 30 40 ILO = 1 IF (X.LT.XT(1)) GO TO 100 GO TO 80 C **** NOW X .GE. XT(ILO) . FIND UPPER BOUND 50 ISTEP = 1 60 ILO = IHI IHI = ILO + ISTEP IF (IHI.GE.LXT) GO TO 70 IF (X.LT.XT(IHI)) GO TO 80 ISTEP = ISTEP*2 GO TO 60 70 IF (X.GE.XT(LXT)) GO TO 120 IHI = LXT C C **** NOW XT(ILO) .LE. X .LT. XT(IHI) . NARROW THE INTERVAL 80 MIDDLE = (ILO+IHI)/2 IF (MIDDLE.EQ.ILO) GO TO 110 C NOTE. IT IS ASSUMED THAT MIDDLE = ILO IN CASE IHI = ILO+1 IF (X.LT.XT(MIDDLE)) GO TO 90 ILO = MIDDLE GO TO 80 90 IHI = MIDDLE GO TO 80 C **** SET OUTPUT AND RETURN 100 MFLAG = -1 ILEFT = 1 RETURN 110 MFLAG = 0 ILEFT = ILO RETURN 120 MFLAG = 1 ILEFT = LXT RETURN END C SUBROUTINE BSPLVD(T,K,X,ILEFT,VNIKX,NDERIV) C CALCULATES VALUE AND DERIV.S OF ALL B-SPLINES WHICH DO NOT VANISH C AT X. FILL VNIKX(J,IDERIV), J=IDERIV, ... ,K WITH NONZERO VALUES C OF B-SPLINES OF ORDER K+1-IDERIV , IDERIV=NDERIV, ... ,1, BY C REPEATED CALLS TO BSPLVN C DIMENSION T(ILEFT+K) C .. Scalar Arguments .. DOUBLE PRECISION X INTEGER ILEFT,K,NDERIV C .. C .. Array Arguments .. DOUBLE PRECISION T(*),VNIKX(K,NDERIV) C .. C .. Local Scalars .. DOUBLE PRECISION FACTOR,FKMD,V INTEGER I,IDERIV,IPKMD,J,JLOW,JP1MID,KMD,KP1,L,LDUMMY,M,MHIGH C .. C .. Local Arrays .. DOUBLE PRECISION A(5,5) C .. C .. External Subroutines .. EXTERNAL BSPLVN C .. C .. Intrinsic Functions .. INTRINSIC DBLE,MAX0,MIN0 C .. IDERIV = MAX0(MIN0(NDERIV,K),1) KP1 = K + 1 CALL BSPLVN(T,KP1-IDERIV,1,X,ILEFT,VNIKX) IF (IDERIV.EQ.1) GO TO 110 MHIGH = IDERIV DO 20 M = 2,MHIGH JP1MID = 1 DO 10 J = IDERIV,K VNIKX(J,IDERIV) = VNIKX(JP1MID,1) JP1MID = JP1MID + 1 10 CONTINUE IDERIV = IDERIV - 1 CALL BSPLVN(T,KP1-IDERIV,2,X,ILEFT,VNIKX) 20 CONTINUE C JLOW = 1 DO 40 I = 1,K DO 30 J = JLOW,K A(I,J) = 0.D0 30 CONTINUE JLOW = I A(I,I) = 1.D0 40 CONTINUE KMD = K DO 100 M = 2,MHIGH KMD = KMD - 1 FKMD = DBLE(KMD) I = ILEFT J = K DO 60 LDUMMY = 1,KMD IPKMD = I + KMD FACTOR = FKMD/(T(IPKMD)-T(I)) DO 50 L = 1,J A(L,J) = (A(L,J)-A(L,J-1))*FACTOR 50 CONTINUE I = I - 1 J = J - 1 60 CONTINUE C 70 DO 90 I = 1,K V = 0.D0 JLOW = MAX0(I,M) DO 80 J = JLOW,K V = A(I,J)*VNIKX(J,M) + V 80 CONTINUE VNIKX(I,M) = V 90 CONTINUE 100 CONTINUE 110 RETURN END C SUBROUTINE BSPLVN(T,JHIGH,INDEX,X,ILEFT,VNIKX) C ALCULATES THE VALUE OF ALL POSSIBLY NONZERO B-SPLINES AT 'X' OF C ORDER MAX(JHIGH,(J+1)(INDEX-1)) ON 'T'. C DIMENSION T(ILEFT+JHIGH) C .. Scalar Arguments .. DOUBLE PRECISION X INTEGER ILEFT,INDEX,JHIGH C .. C .. Array Arguments .. DOUBLE PRECISION T(*),VNIKX(JHIGH) C .. C .. Local Scalars .. DOUBLE PRECISION VM,VMPREV INTEGER IMJP1,IPJ,J,JP1,JP1ML,L C .. C .. Local Arrays .. DOUBLE PRECISION DELTAM(5),DELTAP(5) C .. C .. Data statements .. DATA J/1/ SAVE DELTAM,DELTAP C .. C ONTENT OF J, DELTAM, DELTAP IS EXPECTED UNCHANGED BETWEEN CALLS. GO TO (10,20) INDEX 10 J = 1 VNIKX(1) = 1.D0 IF (J.GE.JHIGH) GO TO 40 C 20 IPJ = ILEFT + J DELTAP(J) = T(IPJ) - X IMJP1 = ILEFT - J + 1 DELTAM(J) = X - T(IMJP1) VMPREV = 0.D0 JP1 = J + 1 DO 30 L = 1,J JP1ML = JP1 - L VM = VNIKX(L)/(DELTAP(L)+DELTAM(JP1ML)) VNIKX(L) = VM*DELTAP(L) + VMPREV VMPREV = VM*DELTAM(JP1ML) 30 CONTINUE VNIKX(JP1) = VMPREV J = JP1 IF (J.LT.JHIGH) GO TO 20 C 40 RETURN END SUBROUTINE SPLBAS3(N,IX,X,F) C C SRZ - 4/05 C C SPLBAS3 COMPUTES THE 3-DIMENSIONAL TENSOR PRODUCT OF BASIS C FUNCTIONS HRMBAS, SPLBAS AND HRMBAS C INCLUDE 'basis.h' INCLUDE 'basis1.h' INCLUDE 'basis3.h' C C C .. Scalar Arguments .. DOUBLE PRECISION X INTEGER IX,N C .. C .. Array Arguments .. DOUBLE PRECISION F(NDIM,*) C .. C .. Subroutine Arguments .. C .. C .. Local Scalars .. INTEGER I,IJ,J,L,L3 C .. C .. Local Arrays .. DOUBLE PRECISION F1(NDIM,4),F2(NDIM,4),F3(NDIM,4) C .. C .. External Subroutines .. EXTERNAL HRMBAS,SPLBAS C .. KS = 4 DO 20 I = 1,N DO 10 J = 1,KS F(I,J) = 0.0D0 10 CONTINUE 20 CONTINUE DO 30 I = 1,N1 + 4 T(I) = T1(I) 30 CONTINUE C PERIODIC B-SPLINE C CALL SPLBASP(N1,IX,X1(IX1(IX)),F1) C HARMONIC BASIS CALL HRMBAS(N1,IX,X1(IX1(IX)),F1) DO 40 I = 1,N2 + 4 T(I) = T2(I) 40 CONTINUE C B-SPLINE CALL SPLBAS(N2,IX,X2(IX2(IX)),F2) DO 50 I = 1,N3 + 4 T(I) = T3(I) 50 CONTINUE C HARMONIC BASIS CALL HRMBAS(N3,IX,X3(IX3(IX)),F3) C PERIODIC B-SPLINE C CALL SPLBASP(N3,IX,X3(IX3(IX)),F3) IJ = 0 DO 100 L3 = 1,N3 DO 90 I = 1,N1 DO 80 J = 1,N2 IJ = IJ + 1 DO 70 L = 1,KS F(IJ,L) = F1(I,L)*F2(J,L)*F3(L3,L) C PARTIAL DERIVATIVES FOR VARIABLE 2 IF (L.GT.1) F(IJ,L) = F1(I,1)*F2(J,L)*F3(L3,1) 70 CONTINUE 80 CONTINUE 90 CONTINUE 100 CONTINUE C C RETURN END C C