V3FIT
r8tpsplinb.f
1 c tpsplinb -- dmc 20 Jan 1999
2 c
3 c set up coefficients for bicubic spline with following BC's:
4 c * 1st (x) dimension BCs under user control
5 c * derivatives periodic in second & third coordinates
6 c
7  subroutine r8tpsplinb(x,inx,th,inth,ph,inph,fspl,inf4,inf5,
8  > ibcxmin,bcxmin,ibcxmax,bcxmax,inb1,
9  > wk,nwk,ilinx,ilinth,ilinph,ier)
10 c
11  IMPLICIT NONE
12  INTEGER, PARAMETER :: R8=selected_real_kind(12,100)
13  integer inx,inth,inph,inf4,inf5,nwk,inb1
14  real*8 x(inx),th(inth),ph(inph)
15  real*8 fspl(4,4,4,inf4,inf5,inph),wk(nwk)
16  integer ibcxmin,ibcxmax
17  real*8 bcxmin(inb1,inph),bcxmax(inb1,inph)
18  integer ilinx,ilinth,ilinph,ier
19 c
20 c input:
21 c x(1...inx) -- abscissae, first dimension of data
22 c th(1...inth) -- abscissae, second (periodic) dimension of data
23 c ph(1...inph) -- abscissae, third (periodic) dimension of data
24 c fspl(1,1,1,1..inx,1..inth,1..inph) -- function values
25 c inf4 -- fspl dimensioning, inf4.ge.inx required.
26 c inf5 -- fspl dimensioning, inf5.ge.inth required.
27 c
28 c boundary conditions input:
29 c ibcxmin -- indicator for boundary condition at x(1):
30 c bcxmin(...) -- boundary condition data
31 c =0 -- use "not a knot", bcxmin(...) ignored
32 c =1 -- match slope, specified at x(1),th(ith),ph(iph) by bcxmin(ith,iph)
33 c =2 -- match 2nd derivative, specified at x(1),th(ith),ph(iph)
34 c by bcxmin(ith,iph
35 c =3 -- boundary condition is slope=0 (df/dx=0) at x(1), all th(j)
36 c =4 -- boundary condition is d2f/dx2=0 at x(1), all th(j)
37 c =5 -- match 1st derivative to 1st divided difference
38 c =6 -- match 2nd derivative to 2nd divided difference
39 c =7 -- match 3rd derivative to 3rd divided difference
40 c (for more detailed definition of BCs 5-7, see the
41 c comments of subroutine mkspline)
42 c NOTE bcxmin(...) referenced ONLY if ibcxmin=1 or ibcxmin=2
43 c
44 c ibcxmax -- indicator for boundary condition at x(nx):
45 c bcxmax(...) -- boundary condition data
46 c (interpolation as with ibcxmin, bcxmin)
47 c
48 c output:
49 c fspl(*,*,*,1..inx,1..inth,1..inph) -- bicubic spline coeffs (4x4)
50 c ...fspl(1,1,1,*,*,*) is not replaced.
51 c
52 c ilinx -- =1 on output if x(inx) pts are nearly evenly spaced (tol=1e-3)
53 c ilinth-- =1 on output if th(inth) evenly spaced (tol=1e-3)
54 c ilinph-- =1 on output if ph(inph) evenly spaced (tol=1e-3)
55 c
56 c ier -- completion code, 0 for normal
57 c
58 c workspace:
59 c wk -- must be at least 5*max(inx,inth,inph) large
60 c nwk -- size of workspace
61 c
62 c---------------------------------
63 c compute tricubic spline of 3d function, given values at the
64 c grid crossing points, f(1,1,1,i,j,k)=f(x(i),th(j),ph(k)).
65 c
66 c on evaluation: for point x btw x(i) and x(i+1), dx=x-x(i)
67 c and th btw th(j) and th(j+1), dt=th-th(j),
68 c and ph btw ph(k) and ph(k+1), dp=ph-ph(k),
69 c
70 c spline =
71 c f(1,1,1,i,j,k)+dx*f(2,1,1,i,j,k)+dx2*f(3,1,1,i,j,k)+dx3*f(4,1,1,i,j,k)
72 c +dt*(f(1,2,1,i,j,k)+dx*f(2,2,1,i,j,k)+dx2*f(3,2,1,i,j,k)+dx3*f(4,2,1,i,j,k))
73 c +dt2*(f(1,3,1,i,j,k)+dx*f(2,3,1,i,j,k)+dx2*f(3,3,1,i,j,k)+dx3*f(4,3,1,i,j,k))
74 c +dt3*(f(1,4,1,i,j,k)+dx*f(2,4,1,i,j,k)+dx2*f(3,4,1,i,j,k)+dx3*f(4,4,1,i,j,k))
75 c +dp*(
76 c f(1,1,2,i,j,k)+dx*f(2,1,2,i,j,k)+dx2*f(3,1,2,i,j,k)+dx3*f(4,1,2,i,j,k)
77 c +dt*(f(1,2,2,i,j,k)+dx*f(2,2,2,i,j,k)+dx2*f(3,2,2,i,j,k)+dx3*f(4,2,2,i,j,k))
78 c +dt2*(f(1,3,2,i,j,k)+dx*f(2,3,2,i,j,k)+dx2*f(3,3,2,i,j,k)+dx3*f(4,3,2,i,j,k))
79 c +dt3*(f(1,4,2,i,j,k)+dx*f(2,4,2,i,j,k)+dx2*f(3,4,2,i,j,k)+dx3*f(4,4,2,i,j,k)))
80 c +dp2*(
81 c f(1,1,3,i,j,k)+dx*f(2,1,3,i,j,k)+dx2*f(3,1,3,i,j,k)+dx3*f(4,1,3,i,j,k)
82 c +dt*(f(1,2,3,i,j,k)+dx*f(2,2,3,i,j,k)+dx2*f(3,2,3,i,j,k)+dx3*f(4,2,3,i,j,k))
83 c +dt2*(f(1,3,3,i,j,k)+dx*f(2,3,3,i,j,k)+dx2*f(3,3,3,i,j,k)+dx3*f(4,3,3,i,j,k))
84 c +dt3*(f(1,4,3,i,j,k)+dx*f(2,4,3,i,j,k)+dx2*f(3,4,3,i,j,k)+dx3*f(4,4,3,i,j,k)))
85 c +dp3*(
86 c f(1,1,4,i,j,k)+dx*f(2,1,4,i,j,k)+dx2*f(3,1,4,i,j,k)+dx3*f(4,1,4,i,j,k)
87 c +dt*(f(1,2,4,i,j,k)+dx*f(2,2,4,i,j,k)+dx2*f(3,2,4,i,j,k)+dx3*f(4,2,4,i,j,k))
88 c +dt2*(f(1,3,4,i,j,k)+dx*f(2,3,4,i,j,k)+dx2*f(3,3,4,i,j,k)+dx3*f(4,3,4,i,j,k))
89 c +dt3*(f(1,4,4,i,j,k)+dx*f(2,4,4,i,j,k)+dx2*f(3,4,4,i,j,k)+dx3*f(4,4,4,i,j,k)))
90 c
91 c where dx2=dx**2 and dx3=dx**3.
92 c where dt2=dt**2 and dt3=dt**3.
93 c where dp2=dp**2 and dp3=dp**3.
94 c
95 c---------------------------------
96  integer ierx,ierth,ierph
97  REAL*8 zdumx(inx)
98 c---------------------------------
99 c
100  ier=0
101  if(nwk.lt.5*max(inx,inth,inph)) then
102  write(6,'('' ?tpsplinb: workspace too small.'')')
103  ier=1
104  endif
105  if(inx.lt.2) then
106  write(6,'('' ?tpsplinb: at least 2 x points required.'')')
107  ier=1
108  endif
109  if(inth.lt.2) then
110  write(6,'('' ?tpsplinb: need at least 2 theta points.'')')
111  ier=1
112  endif
113  if(inph.lt.2) then
114  write(6,'('' ?tpsplinb: need at least 2 phi points.'')')
115  ier=1
116  endif
117 c
118  if((ibcxmin.eq.1).or.(ibcxmax.eq.1).or.(ibcxmin.eq.2).or.
119  > (ibcxmax.eq.2)) then
120  if(inb1.lt.inth) then
121  ier=1
122  write(6,
123  >'('.lt.' ?tpsplinb: 1st dim of bcxmin/max arrays inth'')')
124  endif
125  endif
126 c
127  call ibc_ck(ibcxmin,'tpsplinb','xmin',0,7,ier)
128  call ibc_ck(ibcxmax,'tpsplinb','xmax',0,7,ier)
129 c
130 c check ilinx & x vector
131 c
132  call r8splinck(x,inx,ilinx,1.0e-3_r8,ierx)
133  if(ierx.ne.0) ier=2
134 c
135  if(ier.eq.2) then
136  write(6,'('' ?tpsplinb: x axis not strict ascending'')')
137  endif
138 c
139 c check ilinth & th vector
140 c
141  call r8splinck(th,inth,ilinth,1.0e-3_r8,ierth)
142  if(ierth.ne.0) ier=3
143 c
144  if(ier.eq.3) then
145  write(6,'('' ?tpsplinb: theta axis not strict ascending'')')
146  endif
147 c
148 c check ilinth & th vector
149 c
150  call r8splinck(ph,inph,ilinph,1.0e-3_r8,ierph)
151  if(ierph.ne.0) ier=4
152 c
153  if(ier.eq.4) then
154  write(6,'('' ?tpsplinb: phi axis not strict ascending'')')
155  endif
156 c
157  if(ier.ne.0) return
158 c
159 c------------------------------------
160 c
161  call r8tcspline(x,inx,th,inth,ph,inph,fspl,inf4,inf5,
162  > ibcxmin,bcxmin,ibcxmax,bcxmax,inb1,
163  > -1,zdumx,-1,zdumx,inx,
164  > -1,zdumx,-1,zdumx,inx,
165  > wk,nwk,ilinx,ilinth,ilinph,ier)
166 c
167  return
168  end
Generated on Thu Mar 5 2020 15:49:24 for V3FIT by   doxygen 1.8.17