V3FIT
bcspeval.f
1 c
2 c bcspeval -- eval bicubic spline function and/or derivatives
3 c
4  subroutine bcspeval(xget,yget,iselect,fval,
5  > x,nx,y,ny,ilinx,iliny,f,inf3,ier)
6 c
7  integer iselect(6)
8  integer ilinx,iliny,nx,ny,inf3,ier
9 c
10  real xget,yget
11  real fval(*)
12  real x(nx),y(ny),f(4,4,inf3,ny)
13 c
14 c modification -- dmc 11 Jan 1999 -- remove SAVE stmts; break routine
15 C into these parts:
16 C
17 C bcspevxy -- find grid cell of target pt.
18 C bcspevfn -- evaluate function using output of bcpsevxy
19 C
20 C in cases where multiple functions are defined on the same grid,
21 C time can be saved by using bcspevxy once and then bcspevfn
22 C multiple times.
23 c
24 c input:
25 c (xget,yget) location where interpolated value is desired
26 c x(1).le.xget.le.x(nx) expected
27 c y(1).le.yget.le.y(ny) expected
28 c
29 c iselect select desired output
30 c
31 c iselect(1)=1 -- want function value (f) itself
32 c iselect(2)=1 -- want df/dx
33 c iselect(3)=1 -- want df/dy
34 c iselect(4)=1 -- want d2f/dx2
35 c iselect(5)=1 -- want d2f/dy2
36 c iselect(6)=1 -- want d2f/dxdy
37 c
38 c example: iselect(1)=iselect(2)=iselect(3)=1
39 c f, df/dx, and df/dy all evaluated
40 c iselect(4)=iselect(5)=iselect(6)=0
41 c 2nd derivatives not evaluated.
42 c
43 c the number of non zero values iselect(1:6)
44 c determines the number of outputs...
45 c see fval (output) description.
46 c
47 c new dmc December 2005 -- access to higher derivatives (even if not
48 c continuous-- but can only go up to 3rd derivatives on any one coordinate.
49 c if iselect(1)=3 -- want 3rd derivatives
50 c iselect(2)=1 for d3f/dx3
51 c iselect(3)=1 for d3f/dx2dy
52 c iselect(4)=1 for d3f/dxdy2
53 c iselect(5)=1 for d3f/dy3
54 c number of non-zero values iselect(2:5) gives no. of outputs
55 c if iselect(1)=4 -- want 4th derivatives
56 c iselect(2)=1 for d4f/dx3dy
57 c iselect(3)=1 for d4f/dx2dy2
58 c iselect(4)=1 for d4f/dxdy3
59 c number of non-zero values iselect(2:4) gives no. of outputs
60 c if iselect(1)=5 -- want 5th derivatives
61 c iselect(2)=1 for d5f/dx3dy2
62 c iselect(3)=1 for d5f/dx2dy3
63 c number of non-zero values iselect(2:3) gives no. of outputs
64 c if iselect(1)=6 -- want 6th derivatives
65 c d6f/dx3dy3 -- one value is returned.
66 c
67 c x(1...nx) independent coordinate x, strict ascending
68 c y(1...ny) independent coordinate y, strict ascending
69 c
70 c ilinx -- =1: flag that x is linearly spaced (avoid search for speed)
71 c iliny -- =1: flag that y is linearly spaced (avoid search for speed)
72 c
73 c **CAUTION** actual even spacing of x, y is NOT CHECKED HERE!
74 c
75 c
76 c f the function values (at grid points) and spline coefs
77 c
78 c evaluation formula: for point x btw x(i) and x(i+1), dx=x-x(i)
79 c and y btw y(j) and y(j+1), dy=y-y(j),
80 c
81 c spline value =
82 c f(1,1,i,j) + dx*f(2,1,i,j) + dx**2*f(3,1,i,j) + dx**3*f(4,1,i,j)
83 c +dy*(f(1,2,i,j) + dx*f(2,2,i,j) + dx**2*f(3,2,i,j) + dx**3*f(4,2,i,j))
84 c +d2*(f(1,3,i,j) + dx*f(2,3,i,j) + dx**2*f(3,3,i,j) + dx**3*f(4,3,i,j))
85 c +d3*(f(1,4,i,j) + dx*f(2,4,i,j) + dx**2*f(3,4,i,j) + dx**3*f(4,4,i,j))
86 c
87 c where d2=dy**2 and d3=dy**3.
88 c
89 c output:
90 c up to 6 elements of fval, ordered as follows:
91 c fval(1)=function value or lowest order derivative requested
92 c fval(2)=next order derivative
93 c etc
94 c the ordering is a subset of the sequence given under the "iselect"
95 c description.
96 c
97 c ier = 0 -- successful completion; = 1 -- an error occurred.
98 c
99 c-------------------------------------------------------------------
100 c local
101 c
102  integer :: i(1)
103  integer :: j(1)
104 c
105  real dx(1),dy(1)
106 c
107 c--------------------------
108 c
109  i(1)=0
110  j(1)=0
111 c
112  call bcspevxy(xget,yget,x,nx,y,ny,ilinx,iliny,
113  > i(1),j(1),dx(1),dy(1),ier)
114  if(ier.ne.0) return
115 c
116  call bcspevfn(iselect,1,1,fval,i,j,dx,dy,f,inf3,ny)
117 c
118  return
119  end
120 c
121 c-------------------------------------------------------------------------
122 c bcspevxy -- look up x-y zone
123 c
124 c this is the "first part" of bcspeval, see comments, above.
125 c
126  subroutine bcspevxy(xget,yget,x,nx,y,ny,ilinx,iliny,
127  > i,j,dx,dy,ier)
128 c
129  integer nx,ny ! array dimensions
130 c
131  real xget,yget ! target point
132  real x(nx),y(ny) ! indep. coords.
133 c
134  integer ilinx ! =1: assume x evenly spaced
135  integer iliny ! =1: assume y evenly spaced
136 c
137 c output of bcspevxy
138 c
139  integer i,j ! index to cell containing target pt
140  real dx,dy ! displacement of target pt w/in cell
141  ! dx=x-x(i) dy=y-y(j)
142 C
143  integer ier ! return ier.ne.0 on error
144 c
145 c------------------------------------
146 c
147  ier=0
148 c
149 c range check
150 c
151  zxget=xget
152  zyget=yget
153 
154  if((xget.lt.x(1)).or.(xget.gt.x(nx))) then
155  zxtol=4.0e-7*max(abs(x(1)),abs(x(nx)))
156  if((xget.lt.x(1)-zxtol).or.(xget.gt.x(nx)+zxtol)) then
157  ier=1
158  write(6,1001) xget,x(1),x(nx)
159  1001 format(' ?bcspeval: xget=',1pe11.4,' out of range ',
160  > 1pe11.4,' to ',1pe11.4)
161  else
162  if((xget.lt.x(1)-0.5*zxtol).or.
163  > (xget.gt.x(nx)+0.5*zxtol))
164  > write(6,1011) xget,x(1),x(nx)
165  1011 format(' %bcspeval: xget=',1pe15.8,' beyond range ',
166  > 1pe15.8,' to ',1pe15.8,' (fixup applied)')
167  if(xget.lt.x(1)) then
168  zxget=x(1)
169  else
170  zxget=x(nx)
171  endif
172  endif
173  endif
174  if((yget.lt.y(1)).or.(yget.gt.y(ny))) then
175  zytol=4.0e-7*max(abs(y(1)),abs(y(ny)))
176  if((yget.lt.y(1)-zytol).or.(yget.gt.y(ny)+zytol)) then
177  ier=1
178  write(6,1002) yget,y(1),y(ny)
179  1002 format(' ?bcspeval: yget=',1pe11.4,' out of range ',
180  > 1pe11.4,' to ',1pe11.4)
181  else
182  if((yget.lt.y(1)-0.5*zytol).or.(yget.gt.y(ny)+0.5*zytol))
183  > write(6,1012) yget,y(1),y(ny)
184  1012 format(' %bcspeval: yget=',1pe15.8,' beyond range ',
185  > 1pe15.8,' to ',1pe15.8,' (fixup applied)')
186  if(yget.lt.y(1)) then
187  zyget=y(1)
188  else
189  zyget=y(ny)
190  endif
191  endif
192  endif
193  if(ier.ne.0) return
194 c
195 c now find interval in which target point lies..
196 c
197  nxm=nx-1
198  nym=ny-1
199 c
200  if(ilinx.eq.1) then
201  ii=1+nxm*(zxget-x(1))/(x(nx)-x(1))
202  i=min(nxm, ii)
203  if(zxget.lt.x(i)) then
204  i=i-1
205  else if(zxget.gt.x(i+1)) then
206  i=i+1
207  endif
208  else
209  if((1.le.i).and.(i.lt.nxm)) then
210  if((x(i).le.zxget).and.(zxget.le.x(i+1))) then
211  continue ! already have the zone
212  else
213  call zonfind(x,nx,zxget,i)
214  endif
215  else
216  i=nx/2
217  call zonfind(x,nx,zxget,i)
218  endif
219  endif
220 c
221  if(iliny.eq.1) then
222  jj=1+nym*(zyget-y(1))/(y(ny)-y(1))
223  j=min(nym, jj)
224  if(zyget.lt.y(j)) then
225  j=j-1
226  else if(zyget.gt.y(j+1)) then
227  j=j+1
228  endif
229  else
230  if((1.le.j).and.(j.lt.nym)) then
231  if((y(j).le.zyget).and.(zyget.le.y(j+1))) then
232  continue ! already have the zone
233  else
234  call zonfind(y,ny,zyget,j)
235  endif
236  else
237  j=ny/2
238  call zonfind(y,ny,zyget,j)
239  endif
240  endif
241 c
242  dx=zxget-x(i)
243  dy=zyget-y(j)
244 c
245  return
246  end
247 c------------------------------------------------------------------------
248 c bcspevfn -- OK now evaluate the bicubic spline
249 c
250  subroutine bcspevfn(ict,ivec,ivd,fval,iv,jv,dxv,dyv,f,inf3,ny)
251 c
252 c input:
253  integer ict(6) ! selector:
254 c ict(1)=1 for f (don't evaluate f if ict(1)=0)
255 c ict(2)=1 for df/dx ""
256 c ict(3)=1 for df/dy ""
257 c ict(4)=1 for d2f/dx2
258 c ict(5)=1 for d2f/dy2
259 c ict(6)=1 for d2f/dxdy
260 c
261 c note: if ict(1)=-1, evaluate f,d2f/dx2,d2f/dy2,d4f/dx2dy2
262 c
263 c the number of non zero values ict(1:6)
264 c determines the number of outputs...
265 c
266 c new dmc December 2005 -- access to higher derivatives (even if not
267 c continuous-- but can only go up to 3rd derivatives on any one coordinate.
268 c if ict(1)=3 -- want 3rd derivatives
269 c ict(2)=1 for d3f/dx3
270 c ict(3)=1 for d3f/dx2dy
271 c ict(4)=1 for d3f/dxdy2
272 c ict(5)=1 for d3f/dy3
273 c number of non-zero values ict(2:5) gives no. of outputs
274 c if ict(1)=4 -- want 4th derivatives
275 c ict(2)=1 for d4f/dx3dy
276 c ict(3)=1 for d4f/dx2dy2
277 c ict(4)=1 for d4f/dxdy3
278 c number of non-zero values ict(2:4) gives no. of outputs
279 c if ict(1)=5 -- want 5th derivatives
280 c ict(2)=1 for d5f/dx3dy2
281 c ict(3)=1 for d5f/dx2dy3
282 c number of non-zero values ict(2:3) gives no. of outputs
283 c if ict(1)=6 -- want 6th derivatives
284 c d6f/dx3dy3 -- one value is returned.
285 c
286  integer ivec,ivd ! vector dimensioning
287 c
288 c ivec-- number of vector pts (spline values to look up)
289 c ivd -- 1st dimension of fval, .ge.ivec
290 c
291 c output:
292  real fval(ivd,*) ! output array
293 c
294 c v = index to element in vector;
295 c fval(v,1) = first item requested by ict(...),
296 c fval(v,2) = 2nd item requested, ...etc...
297 c
298 c input:
299  integer iv(ivec),jv(ivec) ! grid cell indices -- vectors
300  real dxv(ivec),dyv(ivec) ! displacements w/in cell -- vectors
301 c
302  integer inf3 ! 3rd dimension of f -- .ge. nx
303  real f(4,4,inf3,ny) ! bicubic fcn spline coeffs array
304 c
305 c usage example:
306 c
307 c 1. for each element (xx(v),yy(v)) in a vector of (x,y) pairs,
308 c find the x and y zone indices and displacements with respect
309 c to the "lower left corner" of the zone; store these in vectors
310 c iv,jv and dxv,dyv.
311 c
312 c 2. set ict(1)=0, ict(2)=1, ict(3)=1, the rest zero -- get only
313 c the 1st derivatives.
314 c
315 c 3. ivec is the length of the vector; ivd is the 1st dimension
316 c of the array fval to receive the output
317 c
318 c real fval(ivd,6)
319 c real xv(ivd),yv(ivd)
320 c integer iv(ivd),jv(ivd)
321 c real dxv(ivd),dyv(ivd)
322 c integer ict(6)
323 c
324 c real fspline(4,4,nx,ny) ! spline coeffs
325 c data ict/0,1,1,0,0,0/ ! this call: want 1st derivatives
326 c ! only ... these will be output to
327 c ! fval(*,1) fval(*,2)
328 c ...
329 c do iv=1,ivec
330 c ... ! find indices and displacements
331 c enddo
332 c call bcspevfn(ict,ivec,ivd,fval,iv,jv,dxv,dyv,fspline,nx,ny)
333 c
334 c-------------------------------------------------------------------
335 c local:
336 c
337  integer v ! vector element index
338 c
339 c OK can now do evaluations
340 c
341  iaval=0 ! fval addressing
342 c
343  if(ict(1).le.2) then
344  if((ict(1).gt.0).or.(ict(1).eq.-1)) then
345 c evaluate f
346  iaval=iaval+1
347  do v=1,ivec
348  i=iv(v)
349  j=jv(v)
350  dx=dxv(v)
351  dy=dyv(v)
352  fval(v,iaval)=
353  > f(1,1,i,j)+dy*(f(1,2,i,j)+dy*(f(1,3,i,j)+dy*f(1,4,i,j)))
354  > +dx*(f(2,1,i,j)+dy*(f(2,2,i,j)+dy*(f(2,3,i,j)+dy*f(2,4,i,j)))
355  > +dx*(f(3,1,i,j)+dy*(f(3,2,i,j)+dy*(f(3,3,i,j)+dy*f(3,4,i,j)))
356  > +dx*(f(4,1,i,j)+dy*(f(4,2,i,j)+dy*(f(4,3,i,j)+dy*f(4,4,i,j))))))
357  enddo
358  endif
359 c
360  if((ict(2).gt.0).and.(ict(1).ne.-1)) then
361 c evaluate df/dx
362  iaval=iaval+1
363  do v=1,ivec
364  i=iv(v)
365  j=jv(v)
366  dx=dxv(v)
367  dy=dyv(v)
368  fval(v,iaval)=
369  > f(2,1,i,j)+dy*(f(2,2,i,j)+dy*(f(2,3,i,j)+dy*f(2,4,i,j)))
370  > +2.0*dx*(
371  > f(3,1,i,j)+dy*(f(3,2,i,j)+dy*(f(3,3,i,j)+dy*f(3,4,i,j)))
372  > +1.5*dx*(
373  > f(4,1,i,j)+dy*(f(4,2,i,j)+dy*(f(4,3,i,j)+dy*f(4,4,i,j)))
374  > ))
375  enddo
376  endif
377 c
378  if((ict(3).gt.0).and.(ict(1).ne.-1)) then
379 c evaluate df/dy
380  iaval=iaval+1
381  do v=1,ivec
382  i=iv(v)
383  j=jv(v)
384  dx=dxv(v)
385  dy=dyv(v)
386  fval(v,iaval)=
387  > f(1,2,i,j)+dy*(2.0*f(1,3,i,j)+dy*3.0*f(1,4,i,j))
388  > +dx*(f(2,2,i,j)+dy*(2.0*f(2,3,i,j)+dy*3.0*f(2,4,i,j))
389  > +dx*(f(3,2,i,j)+dy*(2.0*f(3,3,i,j)+dy*3.0*f(3,4,i,j))
390  > +dx*(f(4,2,i,j)+dy*(2.0*f(4,3,i,j)+dy*3.0*f(4,4,i,j))
391  > )))
392  enddo
393  endif
394 c
395  if((ict(4).gt.0).or.(ict(1).eq.-1)) then
396 c evaluate d2f/dx2
397  iaval=iaval+1
398  do v=1,ivec
399  i=iv(v)
400  j=jv(v)
401  dx=dxv(v)
402  dy=dyv(v)
403  fval(v,iaval)=
404  > 2.0*(
405  > f(3,1,i,j)+dy*(f(3,2,i,j)+dy*(f(3,3,i,j)+dy*f(3,4,i,j))))
406  > +6.0*dx*(
407  > f(4,1,i,j)+dy*(f(4,2,i,j)+dy*(f(4,3,i,j)+dy*f(4,4,i,j))))
408  enddo
409  endif
410 c
411  if((ict(5).gt.0).or.(ict(1).eq.-1)) then
412 c evaluate d2f/dy2
413  iaval=iaval+1
414  do v=1,ivec
415  i=iv(v)
416  j=jv(v)
417  dx=dxv(v)
418  dy=dyv(v)
419  fval(v,iaval)=
420  > 2.0*f(1,3,i,j)+6.0*dy*f(1,4,i,j)
421  > +dx*(2.0*f(2,3,i,j)+6.0*dy*f(2,4,i,j)
422  > +dx*(2.0*f(3,3,i,j)+6.0*dy*f(3,4,i,j)
423  > +dx*(2.0*f(4,3,i,j)+6.0*dy*f(4,4,i,j))))
424  enddo
425  endif
426 c
427  if((ict(6).gt.0).and.(ict(1).ne.-1)) then
428 c evaluate d2f/dxdy
429  iaval=iaval+1
430  do v=1,ivec
431  i=iv(v)
432  j=jv(v)
433  dx=dxv(v)
434  dy=dyv(v)
435  fval(v,iaval)=
436  > f(2,2,i,j)+dy*(2.0*f(2,3,i,j)+dy*3.0*f(2,4,i,j))
437  > +2.*dx*(f(3,2,i,j)+dy*(2.0*f(3,3,i,j)+dy*3.0*f(3,4,i,j))
438  >+1.5*dx*(f(4,2,i,j)+dy*(2.0*f(4,3,i,j)+dy*3.0*f(4,4,i,j))
439  > ))
440  enddo
441  endif
442 c
443  if(ict(1).eq.-1) then
444  iaval=iaval+1
445  do v=1,ivec
446  i=iv(v)
447  j=jv(v)
448  dx=dxv(v)
449  dy=dyv(v)
450  fval(v,iaval)=
451  > 4.0*f(3,3,i,j)+12.0*dy*f(3,4,i,j)
452  > +dx*(12.0*f(4,3,i,j)+36.0*dy*f(4,4,i,j))
453  enddo
454  endif
455 c
456 c-----------------------------------
457 c access to 3rd derivatives
458 c
459  else if(ict(1).eq.3) then
460  if(ict(2).eq.1) then
461 c evaluate d3f/dx3 (not continuous)
462  iaval=iaval+1
463  do v=1,ivec
464  i=iv(v)
465  j=jv(v)
466  dy=dyv(v)
467  fval(v,iaval)=
468  > +6.0*(
469  > f(4,1,i,j)+dy*(f(4,2,i,j)+dy*(f(4,3,i,j)+dy*f(4,4,i,j))))
470  enddo
471  endif
472 c
473  if(ict(3).eq.1) then
474 c evaluate d3f/dx2dy
475  iaval=iaval+1
476  do v=1,ivec
477  i=iv(v)
478  j=jv(v)
479  dx=dxv(v)
480  dy=dyv(v)
481  fval(v,iaval)=
482  > 2.0*(
483  > f(3,2,i,j)+dy*(2.0*f(3,3,i,j)+dy*3.0*f(3,4,i,j)))
484  > +6.0*dx*(
485  > f(4,2,i,j)+dy*(2.0*f(4,3,i,j)+dy*3.0*f(4,4,i,j)))
486  enddo
487  endif
488 c
489  if(ict(4).eq.1) then
490 c evaluate d3f/dxdy2
491  iaval=iaval+1
492  do v=1,ivec
493  i=iv(v)
494  j=jv(v)
495  dx=dxv(v)
496  dy=dyv(v)
497  fval(v,iaval)=
498  > (2.0*f(2,3,i,j)+6.0*dy*f(2,4,i,j)
499  > +2.0*dx*(2.0*f(3,3,i,j)+6.0*dy*f(3,4,i,j)
500  > +1.5*dx*(2.0*f(4,3,i,j)+6.0*dy*f(4,4,i,j))
501  > ))
502  enddo
503  endif
504 
505  if(ict(5).eq.1) then
506 c evaluate d3f/dy3 (not continuous)
507  iaval=iaval+1
508  do v=1,ivec
509  i=iv(v)
510  j=jv(v)
511  dx=dxv(v)
512  fval(v,iaval)=6.0*(f(1,4,i,j)+
513  > dx*(f(2,4,i,j)+dx*(f(3,4,i,j)+dx*f(4,4,i,j))))
514  enddo
515  endif
516 c
517 c-----------------------------------
518 c access to 4th derivatives
519 c
520  else if(ict(1).eq.4) then
521  if(ict(2).eq.1) then
522 c evaluate d4f/dx3dy (not continuous)
523  iaval=iaval+1
524  do v=1,ivec
525  i=iv(v)
526  j=jv(v)
527  dy=dyv(v)
528  fval(v,iaval)=
529  > +6.0*(
530  > f(4,2,i,j)+dy*2.0*(f(4,3,i,j)+dy*1.5*f(4,4,i,j)))
531  enddo
532  endif
533 c
534  if(ict(3).eq.1) then
535 c evaluate d4f/dx2dy2
536  iaval=iaval+1
537  do v=1,ivec
538  i=iv(v)
539  j=jv(v)
540  dx=dxv(v)
541  dy=dyv(v)
542  fval(v,iaval)=
543  > 4.0*f(3,3,i,j)+12.0*dy*f(3,4,i,j)
544  > +dx*(12.0*f(4,3,i,j)+36.0*dy*f(4,4,i,j))
545  enddo
546  endif
547 c
548  if(ict(4).eq.1) then
549 c evaluate d4f/dxdy3 (not continuous)
550  iaval=iaval+1
551  do v=1,ivec
552  i=iv(v)
553  j=jv(v)
554  dx=dxv(v)
555  fval(v,iaval)=
556  > 6.0*(f(2,4,i,j)
557  > +2.0*dx*(f(3,4,i,j)+1.5*dx*f(4,4,i,j)))
558  enddo
559  endif
560 c
561 c-----------------------------------
562 c access to 5th derivatives
563 c
564  else if(ict(1).eq.5) then
565  if(ict(2).eq.1) then
566 c evaluate d5f/dx3dy2 (not continuous)
567  iaval=iaval+1
568  do v=1,ivec
569  i=iv(v)
570  j=jv(v)
571  dy=dyv(v)
572  fval(v,iaval)=
573  > +12.0*(f(4,3,i,j)+dy*3.0*f(4,4,i,j))
574  enddo
575  endif
576 c
577  if(ict(3).eq.1) then
578 c evaluate d5f/dx3dy2 (not continuous)
579  iaval=iaval+1
580  do v=1,ivec
581  i=iv(v)
582  j=jv(v)
583  dx=dxv(v)
584  fval(v,iaval)=
585  > 12.0*(f(3,4,i,j)+dx*3.0*f(4,4,i,j))
586  enddo
587  endif
588 c
589 c-----------------------------------
590 c access to 6th derivatives
591 c
592  else if(ict(1).eq.6) then
593 c evaluate d6f/dx3dy3 (not continuous)
594  iaval=iaval+1
595  do v=1,ivec
596  i=iv(v)
597  j=jv(v)
598  fval(v,iaval)=
599  > 36.0*f(4,4,i,j)
600  enddo
601  endif
602 c
603  return
604  end
605 c----------------------
Generated on Thu Mar 5 2020 15:49:24 for V3FIT by   doxygen 1.8.17