V3FIT
pc3ev.f
1  subroutine pc3ev(xget,yget,zget,x,nx,y,ny,z,nz,
2  > ilinx,iliny,ilinz,
3  > f,inf2,inf3,ict,fval,ier)
4 C
5 C evaluate a trilinear interpolant on a rectilinear grid
6 C derivatives are available, but, not continuous across grid planes.
7 C
8 C this subroutine calls two subroutines:
9 C herm3xyz -- find cell containing (xget,yget,zget)
10 C pc3fcn -- evaluate interpolant function and (optionally) derivatives
11 C
12 C input arguments:
13 C ================
14 C
15  real xget,yget,zget ! target of this interpolation
16  real x(nx) ! ordered x grid
17  real y(ny) ! ordered y grid
18  real z(nz) ! ordered z grid
19  integer ilinx ! ilinx=1 => assume x evenly spaced
20  integer iliny ! iliny=1 => assume y evenly spaced
21  integer ilinz ! ilinz=1 => assume z evenly spaced
22  integer nx, ny, nz
23 C
24  real f(inf2,inf3,nz) ! function data
25 C
26 C f 2nd dimension inf2 must be .ge. nx; 3rd dim inf3 .ge. ny
27 C contents of f:
28 C
29 C f(i,j,k) = f @ x(i),y(j),z(k)
30 C
31  integer ict(8) ! code specifying output desired
32 C
33 C ict(1)=1 -- return f (0, don't)
34 C ict(2)=1 -- return df/dx (0, don't)
35 C ict(3)=1 -- return df/dy (0, don't)
36 C ict(4)=1 -- return df/dz (0, don't)
37 C ict(5)=1 -- return d2f/dxdy (0, don't)
38 C ict(6)=1 -- return d2f/dxdz (0, don't)
39 C ict(7)=1 -- return d2f/dydz (0, don't)
40 C ict(8)=1 -- return d3f/dxdydz (0, don't)
41 C
42 C output arguments:
43 C =================
44 C
45  real fval(*) ! output data
46  integer ier ! error code =0 ==> no error
47 C
48 C fval(1) receives the first output (depends on ict(...) spec)
49 C fval(2) receives the second output (depends on ict(...) spec)
50 C fval(3) receives the third output (depends on ict(...) spec)
51 C fval(4) receives the 4th output (depends on ict(...) spec)
52 C fval(5-8) receive 5th thru 8th outputs (if required by ict(...) spec)
53 C
54 C examples:
55 C on input ict = [1,1,1,1,0,0,0,0]
56 C on output fval= [f,df/dx,df/dy,df/dz]
57 C
58 C on input ict = [1,0,0,0,0,0,0,0]
59 C on output fval= [f] ... elements 2-8 never referenced
60 C
61 C on input ict = [0,1,1,0,0,0,0,0]
62 C on output fval= [df/dx,df/dy] ... elements 3-8 never referenced
63 C
64 C on input ict = [0,0,0,0,1,0,0,0]
65 C on output fval= [d2f/dxdy] ... elements 2-8 never referenced.
66 C
67 C ier -- completion code: 0 means OK
68 C-------------------
69 C local:
70 C
71  integer i,j,k ! cell indices
72 C
73 C normalized displacement from (x(i),y(j)) corner of cell.
74 C xparam=0 @x(i) xparam=1 @x(i+1)
75 C yparam=0 @y(j) yparam=1 @y(j+1)
76 C zparam=0 @z(k) zparam=1 @z(k+1)
77 C
78  real xparam,yparam,zparam
79 C
80 C cell dimensions and
81 C inverse cell dimensions hxi = 1/(x(i+1)-x(i)), hyi = 1/(y(j+1)-y(j))
82 C
83  real hx,hy,hz
84  real hxi,hyi,hzi
85 C
86 C 0 .le. xparam .le. 1
87 C 0 .le. yparam .le. 1
88 C 0 .le. zparam .le. 1
89 C
90 C---------------------------------------------------------------------
91 C
92  call herm3xyz(xget,yget,zget,x,nx,y,ny,z,nz,ilinx,iliny,ilinz,
93  > i,j,k,xparam,yparam,zparam,hx,hxi,hy,hyi,hz,hzi,ier)
94  if(ier.ne.0) return
95 c
96  call pc3fcn(ict,1,1,
97  > fval,i,j,k,xparam,yparam,zparam,
98  > hx,hxi,hy,hyi,hz,hzi,
99  > f,inf2,inf3,nz)
100 C
101  return
102  end
103 C---------------------------------------------------------------------
104 C evaluate trilinear function interpolation -- 3d fcn
105 C --vectorized-- dmc 10 Feb 1999
106 C
107  subroutine pc3fcn(ict,ivec,ivecd,
108  > fval,ii,jj,kk,xparam,yparam,zparam,
109  > hx,hxi,hy,hyi,hz,hzi,
110  > fin,inf2,inf3,nz)
111 C
112  integer ict(8) ! requested output control
113  integer ivec ! vector length
114  integer ivecd ! vector dimension (1st dim of fval)
115 C
116  integer ii(ivec),jj(ivec),kk(ivec) ! target cells (i,j,k)
117  real xparam(ivec),yparam(ivec),zparam(ivec)
118  ! normalized displacements from (i,j,k) corners
119 C
120  real hx(ivec),hy(ivec),hz(ivec) ! grid spacing, and
121  real hxi(ivec),hyi(ivec),hzi(ivec) ! inverse grid spacing
122  ! 1/(x(i+1)-x(i)) & 1/(y(j+1)-y(j)) & 1/(z(k+1)-z(i))
123 C
124  real fin(inf2,inf3,nz) ! interpolant data (cf "pc3ev")
125 C
126  real fval(ivecd,*) ! output returned
127 C
128 C for detailed description of fin, ict and fval see subroutine pc3ev
129 C comments. Note ict is not vectorized; the same output
130 C is expected to be returned for all input vector data points.
131 C
132 C note that the index inputs ii,jj,kk and parameter inputs
133 C xparam,yparam,zparam,hx,hxi,hy,hyi,hz,hzi are vectorized, and the
134 C output array fval has a vector ** 1st dimension ** whose
135 C size must be given as a separate argument
136 C
137 C to use this routine in scalar mode, pass in ivec=ivecd=1
138 C
139 C---------------
140 C
141  real sum
142  integer v
143 C
144  do v=1,ivec
145  i=ii(v)
146  j=jj(v)
147  k=kk(v)
148 C
149 C ...in x direction
150 C
151  xp=xparam(v)
152  xpi=1.0-xp
153 C
154 C ...in y direction
155 C
156  yp=yparam(v)
157  ypi=1.0-yp
158 C
159 C ...in z direction
160 C
161  zp=zparam(v)
162  zpi=1.0-zp
163 C
164  iadr=0
165 C
166 C get desired values:
167 C
168  if(ict(1).eq.1) then
169 C
170 C function value:
171 C
172  iadr=iadr+1
173  sum=zpi*(
174  > xpi*(ypi*fin(i,j,k) +yp*fin(i,j+1,k))+
175  > xp*(ypi*fin(i+1,j,k)+yp*fin(i+1,j+1,k)))
176  > +zp*(
177  > xpi*(ypi*fin(i,j,k+1) +yp*fin(i,j+1,k+1))+
178  > xp*(ypi*fin(i+1,j,k+1)+yp*fin(i+1,j+1,k+1)))
179 C
180  fval(v,iadr)=sum
181  endif
182 C
183  if(ict(2).eq.1) then
184 C
185 C df/dx:
186 C
187  iadr=iadr+1
188  sum=zpi*(
189  > -(ypi*fin(i,j,k) +yp*fin(i,j+1,k))
190  > +(ypi*fin(i+1,j,k)+yp*fin(i+1,j+1,k)))
191  > +zp*(
192  > -(ypi*fin(i,j,k+1) +yp*fin(i,j+1,k+1))
193  > +(ypi*fin(i+1,j,k+1)+yp*fin(i+1,j+1,k+1)))
194 C
195  fval(v,iadr)=sum*hxi(v)
196  endif
197 C
198  if(ict(3).eq.1) then
199 C
200 C df/dy:
201 C
202  iadr=iadr+1
203  sum=zpi*(
204  > xpi*(-fin(i,j,k) +fin(i,j+1,k))+
205  > xp*(-fin(i+1,j,k)+fin(i+1,j+1,k)))
206  > +zp*(
207  > xpi*(-fin(i,j,k+1) +fin(i,j+1,k+1))+
208  > xp*(-fin(i+1,j,k+1)+fin(i+1,j+1,k+1)))
209 C
210  fval(v,iadr)=sum*hyi(v)
211  endif
212 C
213  if(ict(4).eq.1) then
214 C
215 C df/dz:
216 C
217  iadr=iadr+1
218  sum= -(
219  > xpi*(ypi*fin(i,j,k) +yp*fin(i,j+1,k))+
220  > xp*(ypi*fin(i+1,j,k)+yp*fin(i+1,j+1,k)))
221  > +(
222  > xpi*(ypi*fin(i,j,k+1) +yp*fin(i,j+1,k+1))+
223  > xp*(ypi*fin(i+1,j,k+1)+yp*fin(i+1,j+1,k+1)))
224 C
225  fval(v,iadr)=sum*hzi(v)
226  endif
227 C
228  if(ict(5).eq.1) then
229 C
230 C d2f/dxdy:
231 C
232  iadr=iadr+1
233  sum=zpi*(
234  > -(-fin(i,j,k) +fin(i,j+1,k))
235  > +(-fin(i+1,j,k)+fin(i+1,j+1,k)))
236  > +zp*(
237  > -(-fin(i,j,k+1) +fin(i,j+1,k+1))
238  > +(-fin(i+1,j,k+1)+fin(i+1,j+1,k+1)))
239 C
240  fval(v,iadr)=sum*hxi(v)*hyi(v)
241  endif
242 C
243  if(ict(6).eq.1) then
244 C
245 C d2f/dxdz:
246 C
247  iadr=iadr+1
248  sum= -(
249  > -(ypi*fin(i,j,k) +yp*fin(i,j+1,k))
250  > +(ypi*fin(i+1,j,k)+yp*fin(i+1,j+1,k)))
251  > +(
252  > -(ypi*fin(i,j,k+1) +yp*fin(i,j+1,k+1))
253  > +(ypi*fin(i+1,j,k+1)+yp*fin(i+1,j+1,k+1)))
254 C
255  fval(v,iadr)=sum*hxi(v)*hzi(v)
256  endif
257 C
258  if(ict(7).eq.1) then
259 C
260 C d2f/dydz:
261 C
262  iadr=iadr+1
263  sum= -(
264  > xpi*(-fin(i,j,k) +fin(i,j+1,k))+
265  > xp*(-fin(i+1,j,k)+fin(i+1,j+1,k)))
266  > +(
267  > xpi*(-fin(i,j,k+1) +fin(i,j+1,k+1))+
268  > xp*(-fin(i+1,j,k+1)+fin(i+1,j+1,k+1)))
269 C
270  fval(v,iadr)=sum*hyi(v)*hzi(v)
271  endif
272 C
273  if(ict(8).eq.1) then
274 C
275 C d3f/dxdydz:
276 C
277  iadr=iadr+1
278  sum= -(
279  > -(-fin(i,j,k) +fin(i,j+1,k))
280  > +(-fin(i+1,j,k)+fin(i+1,j+1,k)))
281  > +(
282  > -(-fin(i,j,k+1) +fin(i,j+1,k+1))
283  > +(-fin(i+1,j,k+1)+fin(i+1,j+1,k+1)))
284 C
285  fval(v,iadr)=sum*hxi(v)*hyi(v)*hzi(v)
286  endif
287 C
288  enddo ! vector loop
289 C
290  return
291  end
292 C---------------------------------------------------------------------
293 C evaluate trilinear function interpolation -- 3d fcn
294 C --vectorized-- dmc 10 Feb 1999
295 C --optimized for VARIATION along x axis ONLY--
296 C
297  subroutine pc3fcnx(ict,ivec,ivecd,
298  > fval,ii,jj,kk,xparam,yparam,zparam,
299  > hx,hxi,hy,hyi,hz,hzi,
300  > fin,inf2,inf3,nz)
301 C
302  integer ict(8) ! requested output control
303  integer ivec ! vector length
304  integer ivecd ! vector dimension (1st dim of fval)
305 C
306  integer ii(ivec),jj,kk ! target cells (i,j,k)
307  real xparam(ivec),yparam,zparam
308  ! normalized displacements from (i,j,k) corners
309 C
310  real hx(ivec),hy,hz ! grid spacing, and
311  real hxi(ivec),hyi,hzi ! inverse grid spacing
312  ! 1/(x(i+1)-x(i)) & 1/(y(j+1)-y(j)) & 1/(z(k+1)-z(i))
313 C
314  real fin(inf2,inf3,nz) ! interpolant data (cf "pc3ev")
315 C
316  real fval(ivecd,*) ! output returned
317 C
318 C for detailed description of fin, ict and fval see subroutine pc3ev
319 C comments. Note ict is not vectorized; the same output
320 C is expected to be returned for all input vector data points.
321 C
322 C note that the index inputs ii,jj,kk and parameter inputs
323 C xparam,yparam,zparam,hx,hxi,hy,hyi,hz,hzi are vectorized, and the
324 C output array fval has a vector ** 1st dimension ** whose
325 C size must be given as a separate argument
326 C
327 C to use this routine in scalar mode, pass in ivec=ivecd=1
328 C
329 C---------------
330 C
331  real sum
332  integer v
333 C
334  j=jj
335  k=kk
336 C
337 C ...in y direction
338 C
339  yp=yparam
340  ypi=1.0-yp
341 C
342 C ...in z direction
343 C
344  zp=zparam
345  zpi=1.0-zp
346 
347  do v=1,ivec
348  i=ii(v)
349 C
350 C ...in x direction
351 C
352  xp=xparam(v)
353  xpi=1.0-xp
354 C
355  iadr=0
356 C
357 C get desired values:
358 C
359  if(ict(1).eq.1) then
360 C
361 C function value:
362 C
363  iadr=iadr+1
364  sum=zpi*(
365  > xpi*(ypi*fin(i,j,k) +yp*fin(i,j+1,k))+
366  > xp*(ypi*fin(i+1,j,k)+yp*fin(i+1,j+1,k)))
367  > +zp*(
368  > xpi*(ypi*fin(i,j,k+1) +yp*fin(i,j+1,k+1))+
369  > xp*(ypi*fin(i+1,j,k+1)+yp*fin(i+1,j+1,k+1)))
370 C
371  fval(v,iadr)=sum
372  endif
373 C
374  if(ict(2).eq.1) then
375 C
376 C df/dx:
377 C
378  iadr=iadr+1
379  sum=zpi*(
380  > -(ypi*fin(i,j,k) +yp*fin(i,j+1,k))
381  > +(ypi*fin(i+1,j,k)+yp*fin(i+1,j+1,k)))
382  > +zp*(
383  > -(ypi*fin(i,j,k+1) +yp*fin(i,j+1,k+1))
384  > +(ypi*fin(i+1,j,k+1)+yp*fin(i+1,j+1,k+1)))
385 C
386  fval(v,iadr)=sum*hxi(v)
387  endif
388 C
389  if(ict(3).eq.1) then
390 C
391 C df/dy:
392 C
393  iadr=iadr+1
394  sum=zpi*(
395  > xpi*(-fin(i,j,k) +fin(i,j+1,k))+
396  > xp*(-fin(i+1,j,k)+fin(i+1,j+1,k)))
397  > +zp*(
398  > xpi*(-fin(i,j,k+1) +fin(i,j+1,k+1))+
399  > xp*(-fin(i+1,j,k+1)+fin(i+1,j+1,k+1)))
400 C
401  fval(v,iadr)=sum*hyi
402  endif
403 C
404  if(ict(4).eq.1) then
405 C
406 C df/dz:
407 C
408  iadr=iadr+1
409  sum= -(
410  > xpi*(ypi*fin(i,j,k) +yp*fin(i,j+1,k))+
411  > xp*(ypi*fin(i+1,j,k)+yp*fin(i+1,j+1,k)))
412  > +(
413  > xpi*(ypi*fin(i,j,k+1) +yp*fin(i,j+1,k+1))+
414  > xp*(ypi*fin(i+1,j,k+1)+yp*fin(i+1,j+1,k+1)))
415 C
416  fval(v,iadr)=sum*hzi
417  endif
418 C
419  if(ict(5).eq.1) then
420 C
421 C d2f/dxdy:
422 C
423  iadr=iadr+1
424  sum=zpi*(
425  > -(-fin(i,j,k) +fin(i,j+1,k))
426  > +(-fin(i+1,j,k)+fin(i+1,j+1,k)))
427  > +zp*(
428  > -(-fin(i,j,k+1) +fin(i,j+1,k+1))
429  > +(-fin(i+1,j,k+1)+fin(i+1,j+1,k+1)))
430 C
431  fval(v,iadr)=sum*hxi(v)*hyi
432  endif
433 C
434  if(ict(6).eq.1) then
435 C
436 C d2f/dxdz:
437 C
438  iadr=iadr+1
439  sum= -(
440  > -(ypi*fin(i,j,k) +yp*fin(i,j+1,k))
441  > +(ypi*fin(i+1,j,k)+yp*fin(i+1,j+1,k)))
442  > +(
443  > -(ypi*fin(i,j,k+1) +yp*fin(i,j+1,k+1))
444  > +(ypi*fin(i+1,j,k+1)+yp*fin(i+1,j+1,k+1)))
445 C
446  fval(v,iadr)=sum*hxi(v)*hzi
447  endif
448 C
449  if(ict(7).eq.1) then
450 C
451 C d2f/dydz:
452 C
453  iadr=iadr+1
454  sum= -(
455  > xpi*(-fin(i,j,k) +fin(i,j+1,k))+
456  > xp*(-fin(i+1,j,k)+fin(i+1,j+1,k)))
457  > +(
458  > xpi*(-fin(i,j,k+1) +fin(i,j+1,k+1))+
459  > xp*(-fin(i+1,j,k+1)+fin(i+1,j+1,k+1)))
460 C
461  fval(v,iadr)=sum*hyi*hzi
462  endif
463 C
464  if(ict(8).eq.1) then
465 C
466 C d3f/dxdydz:
467 C
468  iadr=iadr+1
469  sum= -(
470  > -(-fin(i,j,k) +fin(i,j+1,k))
471  > +(-fin(i+1,j,k)+fin(i+1,j+1,k)))
472  > +(
473  > -(-fin(i,j,k+1) +fin(i,j+1,k+1))
474  > +(-fin(i+1,j,k+1)+fin(i+1,j+1,k+1)))
475 C
476  fval(v,iadr)=sum*hxi(v)*hyi*hzi
477  endif
478 C
479  enddo ! vector loop
480 C
481  return
482  end
Generated on Thu Mar 5 2020 15:49:24 for V3FIT by   doxygen 1.8.17