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