1      REAL FUNCTION B4VAL(XVAL,YVAL,ZVAL,VVAL,IDX,IDY,IDZ,IDV,
2     *  TX,TY,TZ,TV,NX,NY,NZ,NV,KX,KY,KZ,KV,BCOEF,WORK)
3C***BEGIN PROLOGUE  B4VAL
4C***DATE WRITTEN   25 MAY 1982
5C***REVISION DATE  25 MAY 1982
6C***CATEGORY NO.  E1A
7C***KEYWORDS  INTERPOLATION, THREE-DIMENSIONS, GRIDDED DATA, SPLINES,
8C             PIECEWISE POLYNOMIALS
9C***AUTHOR  BOISVERT, RONALD, NBS
10C             SCIENTIFIC COMPUTING DIVISION
11C             NATIONAL BUREAU OF STANDARDS
12C             WASHINGTON, DC 20234
13C***PURPOSE  B4VAL EVALUATES THE PIECEWISE POLYNOMIAL INTERPOLATING
14C            FUNCTION CONSTRUCTED BY THE ROUTINE B4INK OR ONE OF ITS
15C            PARTIAL DERIVATIVES.
16C***DESCRIPTION
17C
18C   B4VAL evaluates the tensor product piecewise polynomial interpolant
19C   constructed by the routine B4INK or one of its derivatives  at  the
20C   point (XVAL,YVAL,ZVAL). To evaluate  the  interpolant  itself,  set
21C   IDX=IDY=IDZ=0, to evaluate the first partial with respect to x, set
22C   IDX=1,IDY=IDZ=0, and so on.
23C
24C   B4VAL returns 0.0E0 if (XVAL,YVAL,ZVAL) is out of range. That is,
25C            XVAL.LT.TX(1) .OR. XVAL.GT.TX(NX+KX) .OR.
26C            YVAL.LT.TY(1) .OR. YVAL.GT.TY(NY+NY) .OR.
27C            ZVAL.LT.TZ(1) .OR. ZVAL.GT.TZ(NZ+KZ)
28C   If the knots TX, TY, and TZ were chosen  by  B4INK,  then  this  is
29C   equivalent to
30C            XVAL.LT.X(1) .OR. XVAL.GT.X(NX)+EPSX .OR.
31C            YVAL.LT.Y(1) .OR. YVAL.GT.Y(NY)+EPSY .OR.
32C            ZVAL.LT.Z(1) .OR. ZVAL.GT.Z(NZ)+EPSZ
33C   where EPSX = 0.1*(X(NX)-X(NX-1)), EPSY =  0.1*(Y(NY)-Y(NY-1)),  and
34C   EPSZ = 0.1*(Z(NZ)-Z(NZ-1)).
35C
36C   The input quantities TX, TY, TZ, NX, NY, NZ, KX, KY, KZ, and  BCOEF
37C   should remain unchanged since the last call of B4INK.
38C
39C
40C   I N P U T
41C   ---------
42C
43C   XVAL    Real scalar
44C           X coordinate of evaluation point.
45C
46C   YVAL    Real scalar
47C           Y coordinate of evaluation point.
48C
49C   ZVAL    Real scalar
50C           Z coordinate of evaluation point.
51C
52C   IDX     Integer scalar
53C           X derivative of piecewise polynomial to evaluate.
54C
55C   IDY     Integer scalar
56C           Y derivative of piecewise polynomial to evaluate.
57C
58C   IDZ     Integer scalar
59C           Z derivative of piecewise polynomial to evaluate.
60C
61C   TX      Real 1D array (size NX+KX)
62C           Sequence of knots defining the piecewise polynomial in
63C           the x direction.  (Same as in last call to B4INK.)
64C
65C   TY      Real 1D array (size NY+KY)
66C           Sequence of knots defining the piecewise polynomial in
67C           the y direction.  (Same as in last call to B4INK.)
68C
69C   TZ      Real 1D array (size NZ+KZ)
70C           Sequence of knots defining the piecewise polynomial in
71C           the z direction.  (Same as in last call to B4INK.)
72C
73C   NX      Integer scalar
74C           The number of interpolation points in x.
75C           (Same as in last call to B4INK.)
76C
77C   NY      Integer scalar
78C           The number of interpolation points in y.
79C           (Same as in last call to B4INK.)
80C
81C   NZ      Integer scalar
82C           The number of interpolation points in z.
83C           (Same as in last call to B4INK.)
84C
85C   KX      Integer scalar
86C           Order of polynomial pieces in x.
87C           (Same as in last call to B4INK.)
88C
89C   KY      Integer scalar
90C           Order of polynomial pieces in y.
91C           (Same as in last call to B4INK.)
92C
93C   KZ      Integer scalar
94C           Order of polynomial pieces in z.
95C           (Same as in last call to B4INK.)
96C
97C   BCOEF   Real 2D array (size NX by NY by NZ)
98C           The B-spline coefficients computed by B4INK.
99C
100C   WORK    Real 1D array (size KY*KZ+3*max(KX,KY,KZ)+KZ)
101C           A working storage array.
102C
103C***REFERENCES  CARL DE BOOR, A PRACTICAL GUIDE TO SPLINES,
104C                 SPRINGER-VERLAG, NEW YORK, 1978.
105C***ROUTINES CALLED  INTRV,BVALU
106C***END PROLOGUE  B4VAL
107C
108C  ------------
109C  DECLARATIONS
110C  ------------
111C
112C  PARAMETERS
113C
114      INTEGER
115     *        IDX, IDY, IDZ, IDV, NX, NY, NZ, NV, KX, KY, KZ, KV
116      REAL
117     *     XVAL, YVAL, ZVAL, VVAL, TX(1), TY(1), TZ(1), TV(1),
118     *     BCOEF(NX,NY,NZ,NV), WORK(1)
119C
120C  LOCAL VARIABLES
121C
122      INTEGER
123     *        ILOY, ILOZ, ILOV, INBVX, INBV1,
124     *        LEFTY, LEFTZ, LEFTV, MFLAG, KCOLY, KCOLZ, KCOLV,
125     *        IV, IZ, IZM1, IW, I, J, K, L
126      REAL
127     *     BVALU
128C
129      DATA ILOY /1/,  ILOZ /1/,  ILOV /1/, INBVX /1/
130C     SAVE ILOY    ,  ILOZ    ,  ILOV    , INBVX
131C
132C
133C***FIRST EXECUTABLE STATEMENT
134      B4VAL = 0.0E0
135      CALL INTRV(TY,NY+KY,YVAL,ILOY,LEFTY,MFLAG)
136      IF (MFLAG .NE. 0)  GO TO 100
137      CALL INTRV(TZ,NZ+KZ,ZVAL,ILOZ,LEFTZ,MFLAG)
138      IF (MFLAG .NE. 0)  GO TO 100
139      CALL INTRV(TV,NV+KV,VVAL,ILOV,LEFTV,MFLAG)
140      IF (MFLAG .NE. 0)  GO TO 100
141         IV = 1 + KY*KZ*KV
142         IW = IV + KZ*KV
143         KCOLV = LEFTV - KV
144         I = 1
145         DO 40 L=1,KV
146            KCOLV = KCOLV + 1
147            KCOLZ = LEFTZ - KZ
148            DO 40 K=1,KZ
149               KCOLZ = KCOLZ + 1
150               KCOLY = LEFTY - KY
151               DO 40 J=1,KY
152                  KCOLY = KCOLY + 1
153                  WORK(I) = BVALU(TX,BCOEF(1,KCOLY,KCOLZ,KCOLV),NX,KX,
154     *                            IDX,XVAL,INBVX,WORK(IW))
155                  WRITE (*,'(A, I10, E20.10)') 'I=', I, WORK(I)
156                  I = I + 1
157 40      CONTINUE
158         INBV1 = 1
159         KCOLY = LEFTY - KY + 1
160         I = 1
161         J = IV
162         DO 50 L=1,KV
163            DO 50 K=1,KZ
164               WORK(J) = BVALU(TY(KCOLY),WORK(I),KY,KY,IDY,YVAL,
165     *                         INBV1,WORK(IW))
166               WRITE (*,'(A, I10, E20.10)') 'J=', J, WORK(J)
167               J = J + 1
168               I = I + KY
169   50    CONTINUE
170         INBV1 = 1
171         KCOLZ = LEFTZ - KZ + 1
172         J = IV
173         K = 1
174         DO 60 L=1,KV
175            WORK(K) = BVALU(TZ(KCOLZ),WORK(J),KZ,KZ,IDZ,ZVAL,
176     *                      INBV1,WORK(IW))
177            WRITE (*,'(A, I10, E20.10)') 'K=', K, WORK(K)
178            K = K + 1
179            J = J + KZ
180  60     CONTINUE
181         INBV1 = 1
182         KCOLV = LEFTV - KV + 1
183         B4VAL = BVALU(TV(KCOLV),WORK(1),KV,KV,IDV,VVAL,INBV1,WORK(IW))
184         WRITE (*,'(A, I10, E20.10)') 'R=', 0, B4VAL
185  100 CONTINUE
186      RETURN
187      END
188