V3FIT
r8evspline.f
1  subroutine r8evspline(xget,x,nx,ilinx,f,ict,fval,ier)
2 C
3 C use mkspline to set up spline coefficients...
4 C
5 C evaluate a 1d cubic Spline interpolant -- this is C2
6 C
7 C this subroutine calls two subroutines:
8 C herm1x -- find cell containing (xget)
9 C fvspline -- evaluate interpolant function and (optionally) derivatives
10 C
11 C input arguments:
12 C ================
13 C
14  IMPLICIT NONE
15  INTEGER, PARAMETER :: R8=selected_real_kind(12,100)
16  integer nx ! grid size
17  real*8 xget ! target of this interpolation
18  real*8 x(nx) ! ordered x grid
19  integer ilinx ! ilinx=1 => assume x evenly spaced
20 C
21  real*8 f(0:1,nx) ! function data
22 C
23 C f(0,i) = f @ x(i)
24 C f(1,i) = d2f/dx2 @ x(i)
25 C
26 C (these are spline coefficients selected for continuous 2-
27 C diffentiability, see mkspline.for)
28 C
29  integer ict(3) ! 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 ict(3)=1 -- return d2f/dx2 (0, don't)
34 c
35 c set ict(1)=3 to get d3f/dx3 (only)
36 C
37 C output arguments:
38 C =================
39 C
40  real*8 fval(*) ! output data
41  integer ier ! error code =0 ==> no error
42 C
43 C fval(1) receives the first output (depends on ict(...) spec)
44 C fval(2) receives the second output (depends on ict(...) spec)
45 C fval(3) receives the third output (depends on ict(...) spec)
46 C
47 C examples:
48 C on input ict = [1,1,1]
49 C on output fval= [f,df/dx,d2f/dx2]
50 C
51 C on input ict = [1,0,0]
52 C on output fval= [f] ... elements 2 -- 3 never referenced.
53 C
54 C on input ict = [0,0,1]
55 C on output fval= [d2f/dx2] ... elements 2 -- 3 never referenced.
56 C
57 C ier -- completion code: 0 means OK
58 C-------------------
59 C local:
60 C
61  integer i(1) ! cell indices
62 C
63 C normalized displacement from x(i) grid point.
64 C xparam=0 @x(i) xparam=1 @x(i+1)
65 C
66  real*8 xparam(1)
67 C
68 C cell dimensions and
69 C inverse cell dimensions hxi = 1/(x(i+1)-x(i))
70 C
71  real*8 hx(1)
72  real*8 hxi(1)
73 C
74 C 0 .le. xparam .le. 1
75 C
76 C ** the interface is very similar to herm2ev.for; can use herm2xy **
77 C---------------------------------------------------------------------
78 C
79  i(1)=0
80  call r8herm1x(xget,x,nx,ilinx,i(1),xparam(1),hx(1),hxi(1),ier)
81  if(ier.ne.0) return
82 c
83  call r8fvspline(ict,1,1,fval,i,xparam,hx,hxi,f,nx)
84 C
85  return
86  end
87 C---------------------------------------------------------------------
88 C evaluate C1 cubic Hermite function interpolation -- 2d fcn
89 C --vectorized-- dmc 10 Feb 1999
90 C
91  subroutine r8fvspline(ict,ivec,ivecd,
92  > fval,ii,xparam,hx,hxi,fin,nx)
93 C
94 C use mkspline to set up spline coefficients...
95 C
96 !============
97 ! idecl: explicitize implicit INTEGER declarations:
98  IMPLICIT NONE
99  INTEGER, PARAMETER :: R8=selected_real_kind(12,100)
100  INTEGER nx,iadr,i
101 !============
102 ! idecl: explicitize implicit REAL declarations:
103  real*8 xp,xpi,xp2,xpi2,cx,cxi,hx2,cxd,cxdi
104 !============
105  integer ict(3) ! requested output control
106  integer ivec ! vector length
107  integer ivecd ! vector dimension (1st dim of fval)
108 C
109  integer ii(ivec) ! target cells (i)
110  real*8 xparam(ivec)
111  ! normalized displacements from (i) corners
112 C
113  real*8 hx(ivec) ! grid spacing, and
114  real*8 hxi(ivec) ! inverse grid spacing 1/(x(i+1)-x(i))
115 C
116  real*8 fin(0:1,nx) ! interpolant data (cf "evspline")
117 C
118  real*8 fval(ivecd,*) ! output returned
119 C
120 C for detailed description of fin, ict and fval see subroutine
121 C evspline comments. Note ict is not vectorized; the same output
122 C is expected to be returned for all input vector data points.
123 C
124 C note that the index inputs ii and parameter inputs
125 C xparam,hx,hxi, are vectorized, and the
126 C output array fval has a vector ** 1st dimension ** whose
127 C size must be given as a separate argument
128 C
129 C to use this routine in scalar mode, pass in ivec=ivecd=1
130 C
131 C---------------
132 C Spline evaluation consists of a "mixing" of the interpolant
133 C data using the linear functionals xparam, xpi = 1-xparam,
134 C xparam**3-xparam, xpi**3-xpi ...
135 C and their derivatives as needed.
136 C
137  integer v
138 C
139  real*8 sum
140  real*8 sixth
141 C
142  data sixth/0.166666666666666667_r8/
143 C
144 C---------------
145 C ...in x direction
146 C
147  iadr=0
148 C
149  if(ict(1).le.2) then
150 C
151 C normal call
152 C
153  if(ict(1).eq.1) then
154 C
155 C get desired values:
156 C
157  iadr = iadr + 1
158  do v=1,ivec
159  i=ii(v)
160 C
161  xp=xparam(v)
162  xpi=1.0_r8-xp
163  xp2=xp*xp
164  xpi2=xpi*xpi
165 C
166  cx=xp*(xp2-1.0_r8)
167  cxi=xpi*(xpi2-1.0_r8)
168  hx2=hx(v)*hx(v)
169 C
170 C function value:
171 C
172  sum=xpi*fin(0,i) +xp*fin(0,i+1)
173  sum=sum+sixth*hx2*(cxi*fin(1,i) +cx*fin(1,i+1))
174 C
175  fval(v,iadr)=sum
176  enddo
177  endif
178 C
179  if(ict(2).eq.1) then
180 C
181 C df/dx:
182 C
183  iadr=iadr+1
184  do v=1,ivec
185  i=ii(v)
186 C
187  xp=xparam(v)
188  xpi=1.0_r8-xp
189  xp2=xp*xp
190  xpi2=xpi*xpi
191 C
192  cxd=3.0_r8*xp2-1.0_r8
193  cxdi=-3.0_r8*xpi2+1.0_r8
194 
195  sum=hxi(v)*(fin(0,i+1)-fin(0,i))
196  sum=sum+sixth*hx(v)*(cxdi*fin(1,i) +cxd*fin(1,i+1))
197 C
198  fval(v,iadr)=sum
199  enddo
200  endif
201 C
202  if(ict(3).eq.1) then
203 C
204 C d2f/dx2:
205 C
206  iadr=iadr+1
207  do v=1,ivec
208  i=ii(v)
209 C
210  xp=xparam(v)
211  xpi=1.0_r8-xp
212 C
213  sum=xpi*fin(1,i) +xp*fin(1,i+1)
214  fval(v,iadr)=sum
215  enddo
216  endif
217 C
218  else
219 C
220 C return fxxx = d3f/dx3
221 C
222  iadr=1
223  do v=1,ivec
224  i=ii(v)
225  fval(v,iadr)=hxi(v)*(fin(1,i+1)-fin(1,i))
226  enddo
227 C
228  endif
229 C
230  return
231  end
Generated on Thu Mar 5 2020 15:49:24 for V3FIT by   doxygen 1.8.17