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