V3FIT
evbicub.f
1  subroutine evbicub(xget,yget,x,nx,y,ny,ilinx,iliny,
2  > f,inf2,ict,fval,ier)
3 C
4 C evaluate a 2d cubic Spline interpolant on a rectilinear
5 C grid -- this is C2 in both directions.
6 C
7 C this subroutine calls two subroutines:
8 C herm2xy -- find cell containing (xget,yget)
9 C fvbicub -- evaluate interpolant function and (optionally) derivatives
10 C
11 C input arguments:
12 C ================
13 C
14  integer nx,ny ! grid sizes
15  real xget,yget ! target of this interpolation
16  real x(nx) ! ordered x grid
17  real y(ny) ! ordered y grid
18  integer ilinx ! ilinx=1 => assume x evenly spaced
19  integer iliny ! iliny=1 => assume y evenly spaced
20 C
21  real f(0:3,inf2,ny) ! function data
22 C
23 C f 2nd dimension inf2 must be .ge. nx
24 C contents of f:
25 C
26 C f(0,i,j) = f @ x(i),y(j)
27 C f(1,i,j) = d2f/dx2 @ x(i),y(j)
28 C f(2,i,j) = d2f/dy2 @ x(i),y(j)
29 C f(3,i,j) = d4f/dx2dy2 @ x(i),y(j)
30 C
31 C (these are spline coefficients selected for continuous 2-
32 C diffentiability, see mkbicub[w].for)
33 C
34  integer ict(6) ! code specifying output desired
35 C
36 C ict(1)=1 -- return f (0, don't)
37 C ict(2)=1 -- return df/dx (0, don't)
38 C ict(3)=1 -- return df/dy (0, don't)
39 C ict(4)=1 -- return d2f/dx2 (0, don't)
40 C ict(5)=1 -- return d2f/dy2 (0, don't)
41 C ict(6)=1 -- return d2f/dxdy (0, don't)
42 c the number of non zero values ict(1:6)
43 c determines the number of outputs...
44 c
45 c new dmc December 2005 -- access to higher derivatives (even if not
46 c continuous-- but can only go up to 3rd derivatives on any one coordinate.
47 c if ict(1)=3 -- want 3rd derivatives
48 c ict(2)=1 for d3f/dx3
49 c ict(3)=1 for d3f/dx2dy
50 c ict(4)=1 for d3f/dxdy2
51 c ict(5)=1 for d3f/dy3
52 c number of non-zero values ict(2:5) gives no. of outputs
53 c if ict(1)=4 -- want 4th derivatives
54 c ict(2)=1 for d4f/dx3dy
55 c ict(3)=1 for d4f/dx2dy2
56 c ict(4)=1 for d4f/dxdy3
57 c number of non-zero values ict(2:4) gives no. of outputs
58 c if ict(1)=5 -- want 5th derivatives
59 c ict(2)=1 for d5f/dx3dy2
60 c ict(3)=1 for d5f/dx2dy3
61 c number of non-zero values ict(2:3) gives no. of outputs
62 c if ict(1)=6 -- want 6th derivatives
63 c d6f/dx3dy3 -- one value is returned.
64 C
65 C output arguments:
66 C =================
67 C
68  real fval(*) ! output data
69  integer ier ! error code =0 ==> no error
70 C
71 C fval(1) receives the first output (depends on ict(...) spec)
72 C fval(2) receives the second output (depends on ict(...) spec)
73 C fval(3) receives the third output (depends on ict(...) spec)
74 C fval(4) receives the fourth output (depends on ict(...) spec)
75 C fval(5) receives the fourth output (depends on ict(...) spec)
76 C fval(6) receives the fourth output (depends on ict(...) spec)
77 C
78 C examples:
79 C on input ict = [1,1,1,0,0,1]
80 C on output fval= [f,df/dx,df/dy,d2f/dxdy], elements 5 & 6 not referenced.
81 C
82 C on input ict = [1,0,0,0,0,0]
83 C on output fval= [f] ... elements 2 -- 6 never referenced.
84 C
85 C on input ict = [0,0,0,1,1,0]
86 C on output fval= [d2f/dx2,d2f/dy2] ... elements 3 -- 6 never referenced.
87 C
88 C on input ict = [0,0,1,0,0,0]
89 C on output fval= [df/dy] ... elements 2 -- 6 never referenced.
90 C
91 C ier -- completion code: 0 means OK
92 C-------------------
93 C local:
94 C
95  integer i(1),j(1) ! cell indices
96 C
97 C normalized displacement from (x(i),y(j)) corner of cell.
98 C xparam=0 @x(i) xparam=1 @x(i+1)
99 C yparam=0 @y(j) yparam=1 @y(j+1)
100 C
101  real xparam(1),yparam(1)
102 C
103 C cell dimensions and
104 C inverse cell dimensions hxi = 1/(x(i+1)-x(i)), hyi = 1/(y(j+1)-y(j))
105 C
106  real hx(1),hy(1)
107  real hxi(1),hyi(1)
108 C
109 C 0 .le. xparam .le. 1
110 C 0 .le. yparam .le. 1
111 C
112 C ** the interface is very similar to herm2ev.for; can use herm2xy **
113 C---------------------------------------------------------------------
114 C
115  i(1)=0
116  j(1)=0
117  call herm2xy(xget,yget,x,nx,y,ny,ilinx,iliny,
118  > i(1),j(1),xparam(1),yparam(1),hx(1),hxi(1),hy(1),hyi(1),ier)
119  if(ier.ne.0) return
120 c
121  call fvbicub(ict,1,1,
122  > fval,i,j,xparam,yparam,hx,hxi,hy,hyi,f,inf2,ny)
123 C
124  return
125  end
126 C---------------------------------------------------------------------
127 C evaluate C1 cubic Hermite function interpolation -- 2d fcn
128 C --vectorized-- dmc 10 Feb 1999
129 C
130  subroutine fvbicub(ict,ivec,ivecd,
131  > fval,ii,jj,xparam,yparam,hx,hxi,hy,hyi,
132  > fin,inf2,ny)
133 C
134  integer ict(6) ! requested output control
135  integer ivec ! vector length
136  integer ivecd ! vector dimension (1st dim of fval)
137 C
138  integer ii(ivec),jj(ivec) ! target cells (i,j)
139  real xparam(ivec),yparam(ivec)
140  ! normalized displacements from (i,j) corners
141 C
142  real hx(ivec),hy(ivec) ! grid spacing, and
143  real hxi(ivec),hyi(ivec) ! inverse grid spacing 1/(x(i+1)-x(i))
144  ! & 1/(y(j+1)-y(j))
145 C
146  real fin(0:3,inf2,ny) ! interpolant data (cf "evbicub")
147 C
148  real fval(ivecd,*) ! output returned
149 C
150 C for detailed description of fin, ict and fval see subroutine
151 C evbicub comments. Note ict is not vectorized; the same output
152 C is expected to be returned for all input vector data points.
153 C
154 C note that the index inputs ii,jj and parameter inputs
155 C xparam,yparam,hx,hxi,hy,hyi are vectorized, and the
156 C output array fval has a vector ** 1st dimension ** whose
157 C size must be given as a separate argument
158 C
159 C to use this routine in scalar mode, pass in ivec=ivecd=1
160 C
161 C---------------
162 C Spline evaluation consists of a "mixing" of the interpolant
163 C data using the linear functionals xparam, xpi = 1-xparam,
164 C yparam, ypi = 1-yparam, and the cubic functionals
165 C xparam**3-xparam, xpi**3-xpi, yparam**3-yparam, ypi**3-ypi ...
166 C and their derivatives as needed.
167 C
168  integer v
169  real sum
170 C
171  real, parameter :: sixth = 0.166666666666666667
172 C
173 C---------------
174 C ...in x direction
175 C
176  z36th=sixth*sixth
177  iadr=0
178 C
179  if(ict(1).le.2) then
180 C
181 C get desired values:
182 C
183  if(ict(1).eq.1) then
184 C
185 C function value:
186 C
187  iadr=iadr+1
188  do v=1,ivec
189  i=ii(v)
190  j=jj(v)
191 C
192 C in x direction...
193 C
194  xp=xparam(v)
195  xpi=1.0-xp
196  xp2=xp*xp
197  xpi2=xpi*xpi
198 C
199  cx=xp*(xp2-1.0)
200  cxi=xpi*(xpi2-1.0)
201  hx2=hx(v)*hx(v)
202 C
203 C ...and in y direction
204 C
205  yp=yparam(v)
206  ypi=1.0-yp
207  yp2=yp*yp
208  ypi2=ypi*ypi
209 C
210  cy=yp*(yp2-1.0)
211  cyi=ypi*(ypi2-1.0)
212  hy2=hy(v)*hy(v)
213 C
214  sum=xpi*(ypi*fin(0,i,j) +yp*fin(0,i,j+1))+
215  > xp*(ypi*fin(0,i+1,j)+yp*fin(0,i+1,j+1))
216 C
217  sum=sum+sixth*hx2*(
218  > cxi*(ypi*fin(1,i,j) +yp*fin(1,i,j+1))+
219  > cx*(ypi*fin(1,i+1,j)+yp*fin(1,i+1,j+1)))
220 C
221  sum=sum+sixth*hy2*(
222  > xpi*(cyi*fin(2,i,j) +cy*fin(2,i,j+1))+
223  > xp*(cyi*fin(2,i+1,j)+cy*fin(2,i+1,j+1)))
224 C
225  sum=sum+z36th*hx2*hy2*(
226  > cxi*(cyi*fin(3,i,j) +cy*fin(3,i,j+1))+
227  > cx*(cyi*fin(3,i+1,j)+cy*fin(3,i+1,j+1)))
228 C
229  fval(v,iadr)=sum
230  enddo
231  endif
232 C
233  if(ict(2).eq.1) then
234 C
235 C df/dx:
236 C
237  iadr=iadr+1
238  do v=1,ivec
239  i=ii(v)
240  j=jj(v)
241 C
242 C in x direction...
243 C
244  xp=xparam(v)
245  xpi=1.0-xp
246  xp2=xp*xp
247  xpi2=xpi*xpi
248 
249  cxd=3.0*xp2-1.0
250  cxdi=-3.0*xpi2+1.0
251 C
252 C ...and in y direction
253 C
254  yp=yparam(v)
255  ypi=1.0-yp
256  yp2=yp*yp
257  ypi2=ypi*ypi
258 C
259  cy=yp*(yp2-1.0)
260  cyi=ypi*(ypi2-1.0)
261  hy2=hy(v)*hy(v)
262 C
263  sum=hxi(v)*(
264  > -(ypi*fin(0,i,j) +yp*fin(0,i,j+1))
265  > +(ypi*fin(0,i+1,j)+yp*fin(0,i+1,j+1)))
266 C
267  sum=sum+sixth*hx(v)*(
268  > cxdi*(ypi*fin(1,i,j) +yp*fin(1,i,j+1))+
269  > cxd*(ypi*fin(1,i+1,j)+yp*fin(1,i+1,j+1)))
270 C
271  sum=sum+sixth*hxi(v)*hy2*(
272  > -(cyi*fin(2,i,j) +cy*fin(2,i,j+1))
273  > +(cyi*fin(2,i+1,j)+cy*fin(2,i+1,j+1)))
274 C
275  sum=sum+z36th*hx(v)*hy2*(
276  > cxdi*(cyi*fin(3,i,j) +cy*fin(3,i,j+1))+
277  > cxd*(cyi*fin(3,i+1,j)+cy*fin(3,i+1,j+1)))
278 C
279  fval(v,iadr)=sum
280  enddo
281  endif
282 C
283  if(ict(3).eq.1) then
284 C
285 C df/dy:
286 C
287  iadr=iadr+1
288  do v=1,ivec
289  i=ii(v)
290  j=jj(v)
291 C
292 C in x direction...
293 C
294  xp=xparam(v)
295  xpi=1.0-xp
296  xp2=xp*xp
297  xpi2=xpi*xpi
298 C
299  cx=xp*(xp2-1.0)
300  cxi=xpi*(xpi2-1.0)
301  hx2=hx(v)*hx(v)
302 C
303 C ...and in y direction
304 C
305  yp=yparam(v)
306  ypi=1.0-yp
307  yp2=yp*yp
308  ypi2=ypi*ypi
309 
310  cyd=3.0*yp2-1.0
311  cydi=-3.0*ypi2+1.0
312 C
313  sum=hyi(v)*(
314  > xpi*(-fin(0,i,j) +fin(0,i,j+1))+
315  > xp*(-fin(0,i+1,j)+fin(0,i+1,j+1)))
316 C
317  sum=sum+sixth*hx2*hyi(v)*(
318  > cxi*(-fin(1,i,j) +fin(1,i,j+1))+
319  > cx*(-fin(1,i+1,j)+fin(1,i+1,j+1)))
320 C
321  sum=sum+sixth*hy(v)*(
322  > xpi*(cydi*fin(2,i,j) +cyd*fin(2,i,j+1))+
323  > xp*(cydi*fin(2,i+1,j)+cyd*fin(2,i+1,j+1)))
324 C
325  sum=sum+z36th*hx2*hy(v)*(
326  > cxi*(cydi*fin(3,i,j) +cyd*fin(3,i,j+1))+
327  > cx*(cydi*fin(3,i+1,j)+cyd*fin(3,i+1,j+1)))
328 C
329  fval(v,iadr)=sum
330  enddo
331  endif
332 C
333  if(ict(4).eq.1) then
334 C
335 C d2f/dx2:
336 C
337  iadr=iadr+1
338  do v=1,ivec
339  i=ii(v)
340  j=jj(v)
341 C
342 C in x direction...
343 C
344  xp=xparam(v)
345  xpi=1.0-xp
346 C
347 C ...and in y direction
348 C
349  yp=yparam(v)
350  ypi=1.0-yp
351  yp2=yp*yp
352  ypi2=ypi*ypi
353 C
354  cy=yp*(yp2-1.0)
355  cyi=ypi*(ypi2-1.0)
356  hy2=hy(v)*hy(v)
357 C
358  sum=(
359  > xpi*(ypi*fin(1,i,j) +yp*fin(1,i,j+1))+
360  > xp*(ypi*fin(1,i+1,j)+yp*fin(1,i+1,j+1)))
361 C
362  sum=sum+sixth*hy2*(
363  > xpi*(cyi*fin(3,i,j) +cy*fin(3,i,j+1))+
364  > xp*(cyi*fin(3,i+1,j)+cy*fin(3,i+1,j+1)))
365 C
366  fval(v,iadr)=sum
367  enddo
368  endif
369 C
370  if(ict(5).eq.1) then
371 C
372 C d2f/dy2:
373 C
374  iadr=iadr+1
375  do v=1,ivec
376  i=ii(v)
377  j=jj(v)
378 C
379 C in x direction...
380 C
381  xp=xparam(v)
382  xpi=1.0-xp
383  xp2=xp*xp
384  xpi2=xpi*xpi
385 C
386  cx=xp*(xp2-1.0)
387  cxi=xpi*(xpi2-1.0)
388  hx2=hx(v)*hx(v)
389 C
390 C ...and in y direction
391 C
392  yp=yparam(v)
393  ypi=1.0-yp
394 C
395  sum=(
396  > xpi*(ypi*fin(2,i,j) +yp*fin(2,i,j+1))+
397  > xp*(ypi*fin(2,i+1,j)+yp*fin(2,i+1,j+1)))
398 C
399  sum=sum+sixth*hx2*(
400  > cxi*(ypi*fin(3,i,j) +yp*fin(3,i,j+1))+
401  > cx*(ypi*fin(3,i+1,j)+yp*fin(3,i+1,j+1)))
402 C
403  fval(v,iadr)=sum
404  enddo
405  endif
406 C
407  if(ict(6).eq.1) then
408 C
409 C d2f/dxdy:
410 C
411  iadr=iadr+1
412  do v=1,ivec
413  i=ii(v)
414  j=jj(v)
415 C
416 C in x direction...
417 C
418  xp=xparam(v)
419  xpi=1.0-xp
420  xp2=xp*xp
421  xpi2=xpi*xpi
422 
423  cxd=3.0*xp2-1.0
424  cxdi=-3.0*xpi2+1.0
425 C
426 C ...and in y direction
427 C
428  yp=yparam(v)
429  ypi=1.0-yp
430  yp2=yp*yp
431  ypi2=ypi*ypi
432 
433  cyd=3.0*yp2-1.0
434  cydi=-3.0*ypi2+1.0
435 C
436  sum=hxi(v)*hyi(v)*(
437  > fin(0,i,j) -fin(0,i,j+1)
438  > -fin(0,i+1,j)+fin(0,i+1,j+1))
439 C
440  sum=sum+sixth*hx(v)*hyi(v)*(
441  > cxdi*(-fin(1,i,j) +fin(1,i,j+1))+
442  > cxd*(-fin(1,i+1,j)+fin(1,i+1,j+1)))
443 C
444  sum=sum+sixth*hxi(v)*hy(v)*(
445  > -(cydi*fin(2,i,j) +cyd*fin(2,i,j+1))
446  > +(cydi*fin(2,i+1,j)+cyd*fin(2,i+1,j+1)))
447 C
448  sum=sum+z36th*hx(v)*hy(v)*(
449  > cxdi*(cydi*fin(3,i,j) +cyd*fin(3,i,j+1))+
450  > cxd*(cydi*fin(3,i+1,j)+cyd*fin(3,i+1,j+1)))
451 C
452  fval(v,iadr)=sum
453  enddo
454  endif
455 C
456 C-------------------------------------------------
457 C
458  else if(ict(1).eq.3) then
459  if(ict(2).eq.1) then
460 c evaluate d3f/dx3 (not continuous)
461  iadr=iadr+1
462  do v=1,ivec
463  i=ii(v)
464  j=jj(v)
465  yp=yparam(v)
466  ypi=1.0-yp
467  yp2=yp*yp
468  ypi2=ypi*ypi
469  cy=yp*(yp2-1.0)
470  cyi=ypi*(ypi2-1.0)
471  hy2=hy(v)*hy(v)
472  sum=hxi(v)*(
473  > -(ypi*fin(1,i,j) +yp*fin(1,i,j+1))
474  > +(ypi*fin(1,i+1,j)+yp*fin(1,i+1,j+1)))
475 C
476  sum=sum+sixth*hy2*hxi(v)*(
477  > -(cyi*fin(3,i,j) +cy*fin(3,i,j+1))
478  > +(cyi*fin(3,i+1,j)+cy*fin(3,i+1,j+1)))
479 C
480  fval(v,iadr)=sum
481  enddo
482  endif
483 c
484  if(ict(3).eq.1) then
485 c evaluate d3f/dx2dy
486  iadr=iadr+1
487  do v=1,ivec
488  i=ii(v)
489  j=jj(v)
490  xp=xparam(v)
491  xpi=1.0-xp
492  yp=yparam(v)
493  ypi=1.0-yp
494  yp2=yp*yp
495  ypi2=ypi*ypi
496  cyd=3.0*yp2-1.0
497  cydi=-3.0*ypi2+1.0
498 C
499  sum=hyi(v)*(
500  > xpi*(-fin(1,i,j) +fin(1,i,j+1))+
501  > xp*(-fin(1,i+1,j) +fin(1,i+1,j+1)))
502 C
503  sum=sum+sixth*hy(v)*(
504  > xpi*(cydi*fin(3,i,j) +cyd*fin(3,i,j+1))+
505  > xp*(cydi*fin(3,i+1,j)+cyd*fin(3,i+1,j+1)))
506 C
507  fval(v,iadr)=sum
508  enddo
509  endif
510 c
511  if(ict(4).eq.1) then
512 c evaluate d3f/dxdy2
513  iadr=iadr+1
514  do v=1,ivec
515  i=ii(v)
516  j=jj(v)
517  xp=xparam(v)
518  xpi=1.0-xp
519  xp2=xp*xp
520  xpi2=xpi*xpi
521  cxd=3.0*xp2-1.0
522  cxdi=-3.0*xpi2+1.0
523  yp=yparam(v)
524  ypi=1.0-yp
525 C
526  sum=hxi(v)*(
527  > -(ypi*fin(2,i,j) +yp*fin(2,i,j+1))
528  > +(ypi*fin(2,i+1,j)+yp*fin(2,i+1,j+1)))
529 C
530  sum=sum+sixth*hx(v)*(
531  > cxdi*(ypi*fin(3,i,j) +yp*fin(3,i,j+1))+
532  > cxd*(ypi*fin(3,i+1,j)+yp*fin(3,i+1,j+1)))
533 C
534  fval(v,iadr)=sum
535  enddo
536  endif
537 
538  if(ict(5).eq.1) then
539 c evaluate d3f/dy3 (not continuous)
540  iadr=iadr+1
541  do v=1,ivec
542  i=ii(v)
543  j=jj(v)
544 C
545  xp=xparam(v)
546  xpi=1.0-xp
547  xp2=xp*xp
548  xpi2=xpi*xpi
549 C
550  cx=xp*(xp2-1.0)
551  cxi=xpi*(xpi2-1.0)
552  hx2=hx(v)*hx(v)
553 C
554  sum=hyi(v)*(
555  > xpi*(-fin(2,i,j) +fin(2,i,j+1))+
556  > xp*(-fin(2,i+1,j) +fin(2,i+1,j+1)))
557 C
558  sum=sum+sixth*hx2*hyi(v)*(
559  > cxi*(-fin(3,i,j) +fin(3,i,j+1))+
560  > cx*(-fin(3,i+1,j) +fin(3,i+1,j+1)))
561 C
562  fval(v,iadr)=sum
563  enddo
564  endif
565 c
566 c-----------------------------------
567 c access to 4th derivatives
568 c
569  else if(ict(1).eq.4) then
570  if(ict(2).eq.1) then
571 c evaluate d4f/dx3dy (not continuous)
572  iadr=iadr+1
573  do v=1,ivec
574  i=ii(v)
575  j=jj(v)
576  yp=yparam(v)
577  ypi=1.0-yp
578  yp2=yp*yp
579  ypi2=ypi*ypi
580  cyd=3.0*yp2-1.0
581  cydi=-3.0*ypi2+1.0
582 C
583  sum=hxi(v)*hyi(v)*(
584  > +( fin(1,i,j) -fin(1,i,j+1))
585  > +(-fin(1,i+1,j)+fin(1,i+1,j+1)))
586 C
587  sum=sum+sixth*hy(v)*hxi(v)*(
588  > -(cydi*fin(3,i,j) +cyd*fin(3,i,j+1))
589  > +(cydi*fin(3,i+1,j)+cyd*fin(3,i+1,j+1)))
590 C
591  fval(v,iadr)=sum
592  enddo
593  endif
594 c
595  if(ict(3).eq.1) then
596 c evaluate d4f/dx2dy2
597  iadr=iadr+1
598  do v=1,ivec
599  i=ii(v)
600  j=jj(v)
601 C
602  xp=xparam(v)
603  xpi=1.0-xp
604  yp=yparam(v)
605  ypi=1.0-yp
606 C
607  sum=xpi*(ypi*fin(3,i,j) +yp*fin(3,i,j+1))+
608  > xp*(ypi*fin(3,i+1,j)+yp*fin(3,i+1,j+1))
609 C
610  fval(v,iadr)=sum
611  enddo
612  endif
613 c
614  if(ict(4).eq.1) then
615 c evaluate d4f/dxdy3 (not continuous)
616  iadr=iadr+1
617  do v=1,ivec
618  i=ii(v)
619  j=jj(v)
620 C
621  xp=xparam(v)
622  xpi=1.0-xp
623  xp2=xp*xp
624  xpi2=xpi*xpi
625 C
626  cxd=3.0*xp2-1.0
627  cxdi=-3.0*xpi2+1.0
628 C
629  sum=hyi(v)*hxi(v)*(
630  > +( fin(2,i,j) -fin(2,i,j+1))
631  > +(-fin(2,i+1,j)+fin(2,i+1,j+1)))
632 C
633  sum=sum+sixth*hx(v)*hyi(v)*(
634  > cxdi*(-fin(3,i,j) +fin(3,i,j+1))+
635  > cxd*(-fin(3,i+1,j) +fin(3,i+1,j+1)))
636 C
637  fval(v,iadr)=sum
638  enddo
639  endif
640 c
641 c-----------------------------------
642 c access to 5th derivatives
643 c
644  else if(ict(1).eq.5) then
645  if(ict(2).eq.1) then
646 c evaluate d5f/dx3dy2 (not continuous)
647  iadr=iadr+1
648  do v=1,ivec
649  i=ii(v)
650  j=jj(v)
651 C
652  yp=yparam(v)
653  ypi=1.0-yp
654 C
655  sum=hxi(v)*(
656  > -(ypi*fin(3,i,j) +yp*fin(3,i,j+1))
657  > +(ypi*fin(3,i+1,j)+yp*fin(3,i+1,j+1)))
658 C
659  fval(v,iadr)=sum
660  enddo
661  endif
662 c
663  if(ict(3).eq.1) then
664 c evaluate d5f/dx2dy3 (not continuous)
665  iadr=iadr+1
666  do v=1,ivec
667  i=ii(v)
668  j=jj(v)
669 C
670  xp=xparam(v)
671  xpi=1.0-xp
672 C
673  sum=hyi(v)*(
674  > xpi*(-fin(3,i,j) +fin(3,i,j+1))+
675  > xp*(-fin(3,i+1,j)+fin(3,i+1,j+1)))
676 C
677  fval(v,iadr)=sum
678  enddo
679  endif
680 c
681 c-----------------------------------
682 c access to 6th derivatives
683 c
684  else if(ict(1).eq.6) then
685 c evaluate d6f/dx3dy3 (not continuous)
686  iadr=iadr+1
687  do v=1,ivec
688  i=ii(v)
689  j=jj(v)
690  sum=hxi(v)*hyi(v)*(
691  > +( fin(3,i,j) -fin(3,i,j+1))
692  > +(-fin(3,i+1,j)+fin(3,i+1,j+1)))
693  fval(v,iadr)=sum
694  enddo
695  endif
696 c
697  return
698  end
699 C---------------------------------------------------------------------
700 C evaluate C1 cubic Hermite function interpolation -- 2d fcn
701 C --vectorized-- dmc 10 Feb 1999
702 C --optimized for VARIATION along x axis ONLY--
703 C
704  subroutine fvbicubx(ict,ivec,ivecd,
705  > fval,ii,jj,xparam,yparam,hx,hxi,hy,hyi,
706  > fin,inf2,ny)
707 C
708  integer ict(6) ! requested output control
709  integer ivec ! vector length
710  integer ivecd ! vector dimension (1st dim of fval)
711 C
712  integer ii(ivec),jj ! target cells (i,j)
713  real xparam(ivec),yparam
714  ! normalized displacements from (i,j) corners
715 C
716  real hx(ivec),hy ! grid spacing, and
717  real hxi(ivec),hyi ! inverse grid spacing 1/(x(i+1)-x(i))
718  ! & 1/(y(j+1)-y(j))
719 C
720  real fin(0:3,inf2,ny) ! interpolant data (cf "evbicub")
721 C
722  real fval(ivecd,*) ! output returned
723 C
724 C for detailed description of fin, ict and fval see subroutine
725 C evbicub comments. Note ict is not vectorized; the same output
726 C is expected to be returned for all input vector data points.
727 C
728 C note that the index inputs ii,jj and parameter inputs
729 C xparam,yparam,hx,hxi,hy,hyi are vectorized, and the
730 C output array fval has a vector ** 1st dimension ** whose
731 C size must be given as a separate argument
732 C
733 C to use this routine in scalar mode, pass in ivec=ivecd=1
734 C
735 C---------------
736 C Spline evaluation consists of a "mixing" of the interpolant
737 C data using the linear functionals xparam, xpi = 1-xparam,
738 C yparam, ypi = 1-yparam, and the cubic functionals
739 C xparam**3-xparam, xpi**3-xpi, yparam**3-yparam, ypi**3-ypi ...
740 C and their derivatives as needed.
741 C
742  integer v
743  real sum
744 C
745  real, parameter :: sixth = 0.166666666666666667
746 c
747 C---------------
748 C ...in x direction
749 C
750  z36th=sixth*sixth
751  iadr=0
752 C
753  if(ict(1).le.2) then
754 C
755 C get desired values:
756 C
757  if(ict(1).eq.1) then
758 C
759 C function value:
760 C
761  j=jj
762 C
763 C ...and in y direction
764 C
765  yp=yparam
766  ypi=1.0-yp
767  yp2=yp*yp
768  ypi2=ypi*ypi
769 C
770  cy=yp*(yp2-1.0)
771  cyi=ypi*(ypi2-1.0)
772  hy2=hy*hy
773 C
774  iadr=iadr+1
775  do v=1,ivec
776  i=ii(v)
777 C
778 C in x direction...
779 C
780  xp=xparam(v)
781  xpi=1.0-xp
782  xp2=xp*xp
783  xpi2=xpi*xpi
784 C
785  cx=xp*(xp2-1.0)
786  cxi=xpi*(xpi2-1.0)
787  hx2=hx(v)*hx(v)
788 C
789  sum=xpi*(ypi*fin(0,i,j) +yp*fin(0,i,j+1))+
790  > xp*(ypi*fin(0,i+1,j)+yp*fin(0,i+1,j+1))
791 C
792  sum=sum+sixth*hx2*(
793  > cxi*(ypi*fin(1,i,j) +yp*fin(1,i,j+1))+
794  > cx*(ypi*fin(1,i+1,j)+yp*fin(1,i+1,j+1)))
795 C
796  sum=sum+sixth*hy2*(
797  > xpi*(cyi*fin(2,i,j) +cy*fin(2,i,j+1))+
798  > xp*(cyi*fin(2,i+1,j)+cy*fin(2,i+1,j+1)))
799 C
800  sum=sum+z36th*hx2*hy2*(
801  > cxi*(cyi*fin(3,i,j) +cy*fin(3,i,j+1))+
802  > cx*(cyi*fin(3,i+1,j)+cy*fin(3,i+1,j+1)))
803 C
804  fval(v,iadr)=sum
805  enddo
806  endif
807 C
808  if(ict(2).eq.1) then
809 C
810 C df/dx:
811 C
812  j=jj
813 C
814 C ...and in y direction
815 C
816  yp=yparam
817  ypi=1.0-yp
818  yp2=yp*yp
819  ypi2=ypi*ypi
820 C
821  cy=yp*(yp2-1.0)
822  cyi=ypi*(ypi2-1.0)
823  hy2=hy*hy
824 C
825  iadr=iadr+1
826  do v=1,ivec
827  i=ii(v)
828 C
829 C in x direction...
830 C
831  xp=xparam(v)
832  xpi=1.0-xp
833  xp2=xp*xp
834  xpi2=xpi*xpi
835 
836  cxd=3.0*xp2-1.0
837  cxdi=-3.0*xpi2+1.0
838 C
839  sum=hxi(v)*(
840  > -(ypi*fin(0,i,j) +yp*fin(0,i,j+1))
841  > +(ypi*fin(0,i+1,j)+yp*fin(0,i+1,j+1)))
842 C
843  sum=sum+sixth*hx(v)*(
844  > cxdi*(ypi*fin(1,i,j) +yp*fin(1,i,j+1))+
845  > cxd*(ypi*fin(1,i+1,j)+yp*fin(1,i+1,j+1)))
846 C
847  sum=sum+sixth*hxi(v)*hy2*(
848  > -(cyi*fin(2,i,j) +cy*fin(2,i,j+1))
849  > +(cyi*fin(2,i+1,j)+cy*fin(2,i+1,j+1)))
850 C
851  sum=sum+z36th*hx(v)*hy2*(
852  > cxdi*(cyi*fin(3,i,j) +cy*fin(3,i,j+1))+
853  > cxd*(cyi*fin(3,i+1,j)+cy*fin(3,i+1,j+1)))
854 C
855  fval(v,iadr)=sum
856  enddo
857  endif
858 C
859  if(ict(3).eq.1) then
860 C
861 C df/dy:
862 C
863  j=jj
864 C
865 C ...and in y direction
866 C
867  yp=yparam
868  ypi=1.0-yp
869  yp2=yp*yp
870  ypi2=ypi*ypi
871 
872  cyd=3.0*yp2-1.0
873  cydi=-3.0*ypi2+1.0
874 C
875  iadr=iadr+1
876  do v=1,ivec
877  i=ii(v)
878 C
879 C in x direction...
880 C
881  xp=xparam(v)
882  xpi=1.0-xp
883  xp2=xp*xp
884  xpi2=xpi*xpi
885 C
886  cx=xp*(xp2-1.0)
887  cxi=xpi*(xpi2-1.0)
888  hx2=hx(v)*hx(v)
889 C
890  sum=hyi*(
891  > xpi*(-fin(0,i,j) +fin(0,i,j+1))+
892  > xp*(-fin(0,i+1,j)+fin(0,i+1,j+1)))
893 C
894  sum=sum+sixth*hx2*hyi*(
895  > cxi*(-fin(1,i,j) +fin(1,i,j+1))+
896  > cx*(-fin(1,i+1,j)+fin(1,i+1,j+1)))
897 C
898  sum=sum+sixth*hy*(
899  > xpi*(cydi*fin(2,i,j) +cyd*fin(2,i,j+1))+
900  > xp*(cydi*fin(2,i+1,j)+cyd*fin(2,i+1,j+1)))
901 C
902  sum=sum+z36th*hx2*hy*(
903  > cxi*(cydi*fin(3,i,j) +cyd*fin(3,i,j+1))+
904  > cx*(cydi*fin(3,i+1,j)+cyd*fin(3,i+1,j+1)))
905 C
906  fval(v,iadr)=sum
907  enddo
908  endif
909 C
910  if(ict(4).eq.1) then
911 C
912 C d2f/dx2:
913 C
914  j=jj
915 C
916 C ...and in y direction
917 C
918  yp=yparam
919  ypi=1.0-yp
920  yp2=yp*yp
921  ypi2=ypi*ypi
922 C
923  cy=yp*(yp2-1.0)
924  cyi=ypi*(ypi2-1.0)
925  hy2=hy*hy
926 C
927  iadr=iadr+1
928  do v=1,ivec
929  i=ii(v)
930 C
931 C in x direction...
932 C
933  xp=xparam(v)
934  xpi=1.0-xp
935 C
936  sum=(
937  > xpi*(ypi*fin(1,i,j) +yp*fin(1,i,j+1))+
938  > xp*(ypi*fin(1,i+1,j)+yp*fin(1,i+1,j+1)))
939 C
940  sum=sum+sixth*hy2*(
941  > xpi*(cyi*fin(3,i,j) +cy*fin(3,i,j+1))+
942  > xp*(cyi*fin(3,i+1,j)+cy*fin(3,i+1,j+1)))
943 C
944  fval(v,iadr)=sum
945  enddo
946  endif
947 C
948  if(ict(5).eq.1) then
949 C
950 C d2f/dy2:
951 C
952  j=jj
953 C
954 C ...and in y direction
955 C
956  yp=yparam
957  ypi=1.0-yp
958 C
959  iadr=iadr+1
960  do v=1,ivec
961  i=ii(v)
962 C
963 C in x direction...
964 C
965  xp=xparam(v)
966  xpi=1.0-xp
967  xp2=xp*xp
968  xpi2=xpi*xpi
969 C
970  cx=xp*(xp2-1.0)
971  cxi=xpi*(xpi2-1.0)
972  hx2=hx(v)*hx(v)
973 C
974  sum=(
975  > xpi*(ypi*fin(2,i,j) +yp*fin(2,i,j+1))+
976  > xp*(ypi*fin(2,i+1,j)+yp*fin(2,i+1,j+1)))
977 C
978  sum=sum+sixth*hx2*(
979  > cxi*(ypi*fin(3,i,j) +yp*fin(3,i,j+1))+
980  > cx*(ypi*fin(3,i+1,j)+yp*fin(3,i+1,j+1)))
981 C
982  fval(v,iadr)=sum
983  enddo
984  endif
985 C
986  if(ict(6).eq.1) then
987 C
988 C d2f/dxdy:
989 C
990  j=jj
991 C
992 C ...and in y direction
993 C
994  yp=yparam
995  ypi=1.0-yp
996  yp2=yp*yp
997  ypi2=ypi*ypi
998 
999  cyd=3.0*yp2-1.0
1000  cydi=-3.0*ypi2+1.0
1001 C
1002  iadr=iadr+1
1003  do v=1,ivec
1004  i=ii(v)
1005 C
1006 C in x direction...
1007 C
1008  xp=xparam(v)
1009  xpi=1.0-xp
1010  xp2=xp*xp
1011  xpi2=xpi*xpi
1012 
1013  cxd=3.0*xp2-1.0
1014  cxdi=-3.0*xpi2+1.0
1015 C
1016  sum=hxi(v)*hyi*(
1017  > fin(0,i,j) -fin(0,i,j+1)
1018  > -fin(0,i+1,j)+fin(0,i+1,j+1))
1019 C
1020  sum=sum+sixth*hx(v)*hyi*(
1021  > cxdi*(-fin(1,i,j) +fin(1,i,j+1))+
1022  > cxd*(-fin(1,i+1,j)+fin(1,i+1,j+1)))
1023 C
1024  sum=sum+sixth*hxi(v)*hy*(
1025  > -(cydi*fin(2,i,j) +cyd*fin(2,i,j+1))
1026  > +(cydi*fin(2,i+1,j)+cyd*fin(2,i+1,j+1)))
1027 C
1028  sum=sum+z36th*hx(v)*hy*(
1029  > cxdi*(cydi*fin(3,i,j) +cyd*fin(3,i,j+1))+
1030  > cxd*(cydi*fin(3,i+1,j)+cyd*fin(3,i+1,j+1)))
1031 C
1032  fval(v,iadr)=sum
1033  enddo
1034  endif
1035 C
1036 C-------------------------------------------------
1037 C
1038  else if(ict(1).eq.3) then
1039  if(ict(2).eq.1) then
1040 c evaluate d3f/dx3 (not continuous)
1041  j=jj
1042  yp=yparam
1043  ypi=1.0-yp
1044  yp2=yp*yp
1045  ypi2=ypi*ypi
1046  cy=yp*(yp2-1.0)
1047  cyi=ypi*(ypi2-1.0)
1048  hy2=hy*hy
1049 C
1050  iadr=iadr+1
1051  do v=1,ivec
1052  i=ii(v)
1053  sum=hxi(v)*(
1054  > -(ypi*fin(1,i,j) +yp*fin(1,i,j+1))
1055  > +(ypi*fin(1,i+1,j)+yp*fin(1,i+1,j+1)))
1056 C
1057  sum=sum+sixth*hy2*hxi(v)*(
1058  > -(cyi*fin(3,i,j) +cy*fin(3,i,j+1))
1059  > +(cyi*fin(3,i+1,j)+cy*fin(3,i+1,j+1)))
1060 C
1061  fval(v,iadr)=sum
1062  enddo
1063  endif
1064 c
1065  if(ict(3).eq.1) then
1066 c evaluate d3f/dx2dy
1067  j=jj
1068 C
1069  yp=yparam
1070  ypi=1.0-yp
1071  yp2=yp*yp
1072  ypi2=ypi*ypi
1073  cyd=3.0*yp2-1.0
1074  cydi=-3.0*ypi2+1.0
1075 C
1076  iadr=iadr+1
1077  do v=1,ivec
1078  i=ii(v)
1079  xp=xparam(v)
1080  xpi=1.0-xp
1081 C
1082  sum=hyi*(
1083  > xpi*(-fin(1,i,j) +fin(1,i,j+1))+
1084  > xp*(-fin(1,i+1,j) +fin(1,i+1,j+1)))
1085 C
1086  sum=sum+sixth*hy*(
1087  > xpi*(cydi*fin(3,i,j) +cyd*fin(3,i,j+1))+
1088  > xp*(cydi*fin(3,i+1,j)+cyd*fin(3,i+1,j+1)))
1089 C
1090  fval(v,iadr)=sum
1091  enddo
1092  endif
1093 c
1094  if(ict(4).eq.1) then
1095 c evaluate d3f/dxdy2
1096  j=jj
1097  yp=yparam
1098  ypi=1.0-yp
1099 C
1100  iadr=iadr+1
1101  do v=1,ivec
1102  i=ii(v)
1103  xp=xparam(v)
1104  xpi=1.0-xp
1105  xp2=xp*xp
1106  xpi2=xpi*xpi
1107  cxd=3.0*xp2-1.0
1108  cxdi=-3.0*xpi2+1.0
1109 C
1110  sum=hxi(v)*(
1111  > -(ypi*fin(2,i,j) +yp*fin(2,i,j+1))
1112  > +(ypi*fin(2,i+1,j)+yp*fin(2,i+1,j+1)))
1113 C
1114  sum=sum+sixth*hx(v)*(
1115  > cxdi*(ypi*fin(3,i,j) +yp*fin(3,i,j+1))+
1116  > cxd*(ypi*fin(3,i+1,j)+yp*fin(3,i+1,j+1)))
1117 C
1118  fval(v,iadr)=sum
1119  enddo
1120  endif
1121 
1122  if(ict(5).eq.1) then
1123 c evaluate d3f/dy3 (not continuous)
1124  iadr=iadr+1
1125  j=jj
1126  do v=1,ivec
1127  i=ii(v)
1128 C
1129  xp=xparam(v)
1130  xpi=1.0-xp
1131  xp2=xp*xp
1132  xpi2=xpi*xpi
1133 C
1134  cx=xp*(xp2-1.0)
1135  cxi=xpi*(xpi2-1.0)
1136  hx2=hx(v)*hx(v)
1137 C
1138  sum=hyi*(
1139  > xpi*(-fin(2,i,j) +fin(2,i,j+1))+
1140  > xp*(-fin(2,i+1,j) +fin(2,i+1,j+1)))
1141 C
1142  sum=sum+sixth*hx2*hyi*(
1143  > cxi*(-fin(3,i,j) +fin(3,i,j+1))+
1144  > cx*(-fin(3,i+1,j) +fin(3,i+1,j+1)))
1145 C
1146  fval(v,iadr)=sum
1147  enddo
1148  endif
1149 c
1150 c-----------------------------------
1151 c access to 4th derivatives
1152 c
1153  else if(ict(1).eq.4) then
1154  if(ict(2).eq.1) then
1155 c evaluate d4f/dx3dy (not continuous)
1156  iadr=iadr+1
1157  j=jj
1158 C
1159  yp=yparam
1160  ypi=1.0-yp
1161  yp2=yp*yp
1162  ypi2=ypi*ypi
1163  cyd=3.0*yp2-1.0
1164  cydi=-3.0*ypi2+1.0
1165 C
1166  do v=1,ivec
1167  i=ii(v)
1168 C
1169  sum=hxi(v)*hyi*(
1170  > +( fin(1,i,j) -fin(1,i,j+1))
1171  > +(-fin(1,i+1,j)+fin(1,i+1,j+1)))
1172 C
1173  sum=sum+sixth*hy*hxi(v)*(
1174  > -(cydi*fin(3,i,j) +cyd*fin(3,i,j+1))
1175  > +(cydi*fin(3,i+1,j)+cyd*fin(3,i+1,j+1)))
1176 C
1177  fval(v,iadr)=sum
1178  enddo
1179  endif
1180 c
1181  if(ict(3).eq.1) then
1182 c evaluate d4f/dx2dy2
1183  iadr=iadr+1
1184  do v=1,ivec
1185  i=ii(v)
1186  j=jj
1187 C
1188  xp=xparam(v)
1189  xpi=1.0-xp
1190  yp=yparam
1191  ypi=1.0-yp
1192 C
1193  sum=xpi*(ypi*fin(3,i,j) +yp*fin(3,i,j+1))+
1194  > xp*(ypi*fin(3,i+1,j)+yp*fin(3,i+1,j+1))
1195 C
1196  fval(v,iadr)=sum
1197  enddo
1198  endif
1199 c
1200  if(ict(4).eq.1) then
1201 c evaluate d4f/dxdy3 (not continuous)
1202  j=jj
1203 C
1204  iadr=iadr+1
1205  do v=1,ivec
1206  i=ii(v)
1207 C
1208  xp=xparam(v)
1209  xpi=1.0-xp
1210  xp2=xp*xp
1211  xpi2=xpi*xpi
1212 C
1213  cxd=3.0*xp2-1.0
1214  cxdi=-3.0*xpi2+1.0
1215 C
1216  sum=hyi*hxi(v)*(
1217  > +( fin(2,i,j) -fin(2,i,j+1))
1218  > +(-fin(2,i+1,j)+fin(2,i+1,j+1)))
1219 C
1220  sum=sum+sixth*hx(v)*hyi*(
1221  > cxdi*(-fin(3,i,j) +fin(3,i,j+1))+
1222  > cxd*(-fin(3,i+1,j) +fin(3,i+1,j+1)))
1223 C
1224  fval(v,iadr)=sum
1225  enddo
1226  endif
1227 c
1228 c-----------------------------------
1229 c access to 5th derivatives
1230 c
1231  else if(ict(1).eq.5) then
1232  if(ict(2).eq.1) then
1233 c evaluate d5f/dx3dy2 (not continuous)
1234  j=jj
1235 C
1236  yp=yparam
1237  ypi=1.0-yp
1238 C
1239  iadr=iadr+1
1240  do v=1,ivec
1241  i=ii(v)
1242 C
1243  sum=hxi(v)*(
1244  > -(ypi*fin(3,i,j) +yp*fin(3,i,j+1))
1245  > +(ypi*fin(3,i+1,j)+yp*fin(3,i+1,j+1)))
1246 C
1247  fval(v,iadr)=sum
1248  enddo
1249  endif
1250 c
1251  if(ict(3).eq.1) then
1252 c evaluate d5f/dx2dy3 (not continuous)
1253  j=jj
1254  iadr=iadr+1
1255  do v=1,ivec
1256  i=ii(v)
1257 C
1258  xp=xparam(v)
1259  xpi=1.0-xp
1260 C
1261  sum=hyi*(
1262  > xpi*(-fin(3,i,j) +fin(3,i,j+1))+
1263  > xp*(-fin(3,i+1,j)+fin(3,i+1,j+1)))
1264 C
1265  fval(v,iadr)=sum
1266  enddo
1267  endif
1268 c
1269 c-----------------------------------
1270 c access to 6th derivatives
1271 c
1272  else if(ict(1).eq.6) then
1273 c evaluate d6f/dx3dy3 (not continuous)
1274  iadr=iadr+1
1275  j=jj
1276  do v=1,ivec
1277  i=ii(v)
1278  sum=hxi(v)*hyi*(
1279  > +( fin(3,i,j) -fin(3,i,j+1))
1280  > +(-fin(3,i+1,j)+fin(3,i+1,j+1)))
1281  fval(v,iadr)=sum
1282  enddo
1283  endif
1284 c
1285  return
1286  end
Generated on Thu Mar 5 2020 15:49:24 for V3FIT by   doxygen 1.8.17