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