V3FIT
r8herm1ev.f
1  subroutine r8herm1ev(xget,x,nx,ilinx,f,ict,fval,ier)
2 C
3 C evaluate a 1d cubic Hermite interpolant -- this is C1.
4 C
5 C this subroutine calls two subroutines:
6 C herm1x -- find cell containing (xget,yget)
7 C herm1fcn -- evaluate interpolant function and (optionally) derivatives
8 C
9 C input arguments:
10 C ================
11 C
12 !============
13 ! idecl: explicitize implicit INTEGER declarations:
14  IMPLICIT NONE
15  INTEGER, PARAMETER :: R8=selected_real_kind(12,100)
16  INTEGER nx
17 !============
18  real*8 xget ! target of this interpolation
19  real*8 x(nx) ! ordered x grid
20  integer ilinx ! ilinx=1 => assume x evenly spaced
21 C
22  real*8 f(0:1,nx) ! function data
23 C
24 C contents of f:
25 C
26 C f(0,i) = f @ x(i)
27 C f(1,i) = df/dx @ x(i)
28 C
29  integer ict(2) ! code specifying output desired
30 C
31 C ict(1)=1 -- return f (0, don't)
32 C ict(2)=1 -- return df/dx (0, don't)
33 C
34 C output arguments:
35 C =================
36 C
37  real*8 fval(*) ! output data
38  integer ier ! error code =0 ==> no error
39 C
40 C fval(1) receives the first output (depends on ict(...) spec)
41 C fval(2) receives the second output (depends on ict(...) spec)
42 C
43 C examples:
44 C on input ict = [1,1]
45 C on output fval= [f,df/dx]
46 C
47 C on input ict = [1,0]
48 C on output fval= [f] ... element 2 never referenced
49 C
50 C on input ict = [0,1]
51 C on output fval= [df/dx] ... element 2 never referenced
52 C
53 C ier -- completion code: 0 means OK
54 C-------------------
55 C local:
56 C
57  integer i ! cell index
58 C
59 C normalized displacement from (x(i)) corner of cell.
60 C xparam=0 @x(i) xparam=1 @x(i+1)
61 C
62  real*8 xparam
63 C
64 C cell dimensions and
65 C inverse cell dimensions hxi = 1/(x(i+1)-x(i))
66 C
67  real*8 hx
68  real*8 hxi
69 C
70 C 0 .le. xparam .le. 1
71 C
72 C---------------------------------------------------------------------
73 C
74  call r8herm1x(xget,x,nx,ilinx,i,xparam,hx,hxi,ier)
75  if(ier.ne.0) return
76 c
77  call r8herm1fcn(ict,1,1,fval,i,xparam,hx,hxi,f,nx)
78 C
79  return
80  end
81 C---------------------------------------------------------------------
82 c herm1xy -- look up x zone
83 c
84 c this is the "first part" of herm1ev, see comments, above.
85 c
86  subroutine r8herm1x(xget,x,nx,ilinx,i,xparam,hx,hxi,ier)
87 c
88 c input of herm1x
89 c ===============
90 c
91 !============
92 ! idecl: explicitize implicit INTEGER declarations:
93  IMPLICIT NONE
94  INTEGER, PARAMETER :: R8=selected_real_kind(12,100)
95  INTEGER nxm,ii
96 !============
97 ! idecl: explicitize implicit REAL declarations:
98  real*8 zxget,zxtol
99 !============
100  integer nx ! x array dimension
101 c
102  real*8 xget ! target point
103  real*8 x(nx) ! indep. coordinate, strict ascending
104 c
105  integer ilinx ! =1: x evenly spaced
106 c
107 c output of herm1x
108 c =================
109  integer i ! index to cell containing target pt
110 c on exit: 1.le.i.le.nx-1
111 c
112 c normalized position w/in (i) cell (always btw 0 and 1):
113 c
114  real*8 xparam ! (xget-x(i))/(x(i+1)-x(i))
115 c
116 c grid spacing
117 c
118  real*8 hx ! hx = x(i+1)-x(i)
119 c
120 c inverse grid spacing:
121 c
122  real*8 hxi ! 1/hx = 1/(x(i+1)-x(i))
123 c
124  integer ier ! return ier.ne.0 on error
125 c
126 c------------------------------------
127 c
128  ier=0
129 c
130 c range check
131 c
132  zxget=xget
133  if((xget.lt.x(1)).or.(xget.gt.x(nx))) then
134  zxtol=4.0e-7_r8*max(abs(x(1)),abs(x(nx)))
135  if((xget.lt.x(1)-zxtol).or.(xget.gt.x(nx)+zxtol)) then
136  ier=1
137  write(6,1001) xget,x(1),x(nx)
138  1001 format(' ?herm1ev: xget=',1pe11.4,' out of range ',
139  > 1pe11.4,' to ',1pe11.4)
140  else
141  if((xget.lt.x(1)-0.5_r8*zxtol).or.
142  > (xget.gt.x(nx)+0.5_r8*zxtol))
143  > write(6,1011) xget,x(1),x(nx)
144  1011 format(' %herm1ev: xget=',1pe15.8,' beyond range ',
145  > 1pe15.8,' to ',1pe15.8,' (fixup applied)')
146  if(xget.lt.x(1)) then
147  zxget=x(1)
148  else
149  zxget=x(nx)
150  endif
151  endif
152  endif
153  if(ier.ne.0) return
154 c
155 c now find interval in which target point lies..
156 c
157  nxm=nx-1
158 c
159  if(ilinx.eq.1) then
160  ii=1+nxm*(zxget-x(1))/(x(nx)-x(1))
161  i=min(nxm, ii)
162  if(zxget.lt.x(i)) then
163  i=i-1
164  else if(zxget.gt.x(i+1)) then
165  i=i+1
166  endif
167  else
168  if((1.le.i).and.(i.lt.nxm)) then
169  if((x(i).le.zxget).and.(zxget.le.x(i+1))) then
170  continue ! already have the zone
171  else
172  call r8zonfind(x,nx,zxget,i)
173  endif
174  else
175  i=nx/2
176  call r8zonfind(x,nx,zxget,i)
177  endif
178  endif
179 c
180  hx=(x(i+1)-x(i))
181 c
182  hxi=1.0_r8/hx
183 c
184  xparam=(zxget-x(i))*hxi
185 c
186  return
187  end
188 C---------------------------------------------------------------------
189 C evaluate C1 cubic Hermite function interpolation -- 1d fcn
190 C --vectorized-- dmc 10 Feb 1999
191 C
192  subroutine r8herm1fcn(ict,ivec,ivecd,
193  > fval,ii,xparam,hx,hxi,fin,nx)
194 C
195 C input:
196 C
197 !============
198 ! idecl: explicitize implicit INTEGER declarations:
199  IMPLICIT NONE
200  INTEGER, PARAMETER :: R8=selected_real_kind(12,100)
201  INTEGER nx,i,iadr
202 !============
203 ! idecl: explicitize implicit REAL declarations:
204  real*8 xp,xpi,xp2,xpi2,ax,axbar,bx,bxbar,axp,axbarp,bxp,bxbarp
205 !============
206  integer ict(2) ! requested output control
207  integer ivec ! vector length
208  integer ivecd ! vector dimension (1st dim of fval)
209 C
210  integer ii(ivec) ! target cells
211  real*8 xparam(ivec)
212  ! normalized displacements in cells
213 C
214  real*8 hx(ivec) ! grid spacing, and
215  real*8 hxi(ivec) ! inverse grid spacing 1/(x(i+1)-x(i))
216 C
217  real*8 fin(0:1,nx) ! Hermite dataset
218 C
219 C output:
220 C
221  real*8 fval(ivecd,*) ! interpolation results
222 C
223 C for detailed description of fin, ict and fval see subroutine
224 C herm1ev comments. Note ict is not vectorized -- the same output
225 C is expected to be returned for all input vector data points.
226 C
227 C note that the index inputs ii,jj and parameter inputs
228 C xparam,yparam,hx,hxi,hy,hyi are vectorized, and the
229 C
230 C output array fval has a vector ** 1st dimension ** whose
231 C size must be given as a separate argument; ivecd.ge.ivec
232 C expected!
233 C
234 C to use this routine in scalar mode, pass in ivec=ivecd=1
235 C
236 C---------------
237 C Hermite cubic basis functions
238 C -->for function value matching
239 C a(0)=0 a(1)=1 a'(0)=0 a'(1)=0
240 C abar(0)=1 abar(1)=0 abar'(0)=0 abar'(1)=0
241 C
242 C a(x)=-2*x**3 + 3*x**2 = x*x*(-2.0*x+3.0)
243 C abar(x)=1-a(x)
244 C a'(x)=-abar'(x) = 6.0*x*(1.0-x)
245 C
246 C -->for derivative matching
247 C b(0)=0 b(1)=0 b'(0)=0 b'(1)=1
248 C bbar(0)=0 bbar(1)=0 bbar'(0)=1 bbar'(1)=0
249 C
250 C b(x)=x**3-x**2 b'(x)=3*x**2-2*x
251 C bbar(x)=x**3-2*x**2+x bbar'(x)=3*x**2-4*x+1
252 C
253 C ...in x direction
254 C
255  real*8 sum
256  integer v
257 C
258  do v=1,ivec
259  i=ii(v)
260  xp=xparam(v)
261  xpi=1.0_r8-xp
262  xp2=xp*xp
263  xpi2=xpi*xpi
264 C
265  iadr=0
266 C
267 C get desired values:
268 C
269  if(ict(1).eq.1) then
270 C
271 C function value:
272 C
273  iadr=iadr+1
274 C
275 C Hermite basis functions
276 C
277  ax=xp2*(3.0_r8-2.0_r8*xp)
278  axbar=1.0_r8-ax
279  bx=-xp2*xpi
280  bxbar=xpi2*xp
281 C
282  sum=axbar*fin(0,i) + ax*fin(0,i+1)
283 C
284  sum=sum+hx(v)*(bxbar*fin(1,i) + bx*fin(1,i+1))
285 C
286  fval(v,iadr)=sum
287  endif
288 C
289  if(ict(2).eq.1) then
290 C
291 C df/dx:
292 C
293  iadr=iadr+1
294 C
295  axp=6.0_r8*xp*xpi
296  axbarp=-axp
297  bxp=xp*(3.0_r8*xp-2.0_r8)
298  bxbarp=xpi*(3.0_r8*xpi-2.0_r8)
299 C
300  sum=hxi(v)*(axbarp*fin(0,i) +axp*fin(0,i+1))
301 C
302  sum=sum+ bxbarp*fin(1,i) + bxp*fin(1,i+1)
303 C
304  fval(v,iadr)=sum
305  endif
306 C
307  enddo ! vector loop
308 C
309  return
310  end