V3FIT
tpspline.f
1 c tpspline -- dmc 20 Jan 1999
2 c
3 c set up coefficients for bicubic spline with following BC's:
4 c * 1st dimension: 3rd divided difference
5 c * derivatives periodic in second coordinate
6 c * derivatives periodic in third coordinate
7 c
8 c for more control over boundary conditions, use tpspline or tcspline.
9 c
10 c for evaluation of interpolant, see tcspeval.for
11 c
12  subroutine tpspline(x,inx,th,inth,ph,inph,fspl,inf4,inf5,
13  > wk,nwk,ilinx,ilinth,ilinph,ier)
14 c
15  implicit NONE
16  integer inx,inth,inph,inf4,inf5,nwk
17  real x(inx),th(inth),ph(inph)
18  real fspl(4,4,4,inf4,inf5,inph),wk(nwk)
19  integer ilinx,ilinth,ilinph,ier
20 c
21 c input:
22 c x(1...inx) -- abscissae, first dimension of data
23 c th(1...inth) -- abscissae, second (periodic) dimension of data
24 c ph(1...inph) -- abscissae, third (periodic) dimension of data
25 c fspl(1,1,1,1..inx,1..inth,1..inph) -- function values
26 c inf4 -- fspl dimensioning, inf4.ge.inx required.
27 c inf5 -- fspl dimensioning, inf5.ge.inth required.
28 c
29 c output:
30 c fspl(*,*,*,1..inx,1..inth,1..inph) -- bicubic spline coeffs (4x4)
31 c ...fspl(1,1,1,*,*,*) is not replaced.
32 c
33 c ilinx -- =1 on output if x(inx) pts are nearly evenly spaced (tol=1e-3)
34 c ilinth-- =1 on output if th(inth) evenly spaced (tol=1e-3)
35 c ilinph-- =1 on output if ph(inph) evenly spaced (tol=1e-3)
36 c
37 c ier -- completion code, 0 for normal
38 c
39 c workspace:
40 c wk -- must be at least 5*max(inx,inth,inph) large
41 c nwk -- size of workspace
42 c
43 c---------------------------------
44 c compute tricubic spline of 3d function, given values at the
45 c grid crossing points, f(1,1,1,i,j,k)=f(x(i),th(j),ph(k)).
46 c
47 c on evaluation: for point x btw x(i) and x(i+1), dx=x-x(i)
48 c and th btw th(j) and th(j+1), dt=th-th(j),
49 c and ph btw ph(k) and ph(k+1), dp=ph-ph(k),
50 c
51 c spline =
52 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)
53 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))
54 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))
55 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))
56 c +dp*(
57 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)
58 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))
59 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))
60 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)))
61 c +dp2*(
62 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)
63 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))
64 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))
65 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)))
66 c +dp3*(
67 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)
68 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))
69 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))
70 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)))
71 c
72 c where dx2=dx**2 and dx3=dx**3.
73 c where dt2=dt**2 and dt3=dt**3.
74 c where dp2=dp**2 and dp3=dp**3.
75 c
76 c---------------------------------
77  integer ierx,ierth,ierph
78  real zdumth(inth),zdumx(inx)
79 c---------------------------------
80 c
81  ier=0
82  if(nwk.lt.5*max(inx,inth,inph)) then
83  write(6,'('' ?tpspline: workspace too small.'')')
84  ier=1
85  endif
86  if(inx.lt.2) then
87  write(6,'('' ?tpspline: at least 2 x points required.'')')
88  ier=1
89  endif
90  if(inth.lt.2) then
91  write(6,'('' ?tpspline: need at least 2 theta points.'')')
92  ier=1
93  endif
94  if(inph.lt.2) then
95  write(6,'('' ?tpspline: need at least 2 phi points.'')')
96  ier=1
97  endif
98 c
99 c check ilinx & x vector
100 c
101  call splinck(x,inx,ilinx,1.0e-3,ierx)
102  if(ierx.ne.0) ier=2
103 c
104  if(ier.eq.2) then
105  write(6,'('' ?tpspline: x axis not strict ascending'')')
106  endif
107 c
108 c check ilinth & th vector
109 c
110  call splinck(th,inth,ilinth,1.0e-3,ierth)
111  if(ierth.ne.0) ier=3
112 c
113  if(ier.eq.3) then
114  write(6,'('' ?tpspline: theta axis not strict ascending'')')
115  endif
116 c
117 c check ilinth & th vector
118 c
119  call splinck(ph,inph,ilinph,1.0e-3,ierph)
120  if(ierph.ne.0) ier=4
121 c
122  if(ier.eq.4) then
123  write(6,'('' ?tpspline: phi axis not strict ascending'')')
124  endif
125 c
126  if(ier.ne.0) return
127 c
128 c------------------------------------
129 c
130  call tcspline(x,inx,th,inth,ph,inph,fspl,inf4,inf5,
131  > 7,zdumth,7,zdumth,inth,
132  > -1,zdumx,-1,zdumx,inx,
133  > -1,zdumx,-1,zdumx,inx,
134  > wk,nwk,ilinx,ilinth,ilinph,ier)
135 c
136  return
137  end