##// END OF EJS Templates
DATABASES » madrigal
RSS
Add restart to nginx container
Add restart to nginx container
jespinoza - - Load All Authors

File last commit:

r0:b84e1135c2c4
r22:1c8c96a4f254 master
Show More
basis.f
878 lines | 20.9 KiB | text/x-fortran | FortranFixedLexer
/ source / madf / geolib / basis.f
History | Source | Raw |Copy content |Copy permalink
jespinoza
Initial
 r0 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