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