V3FIT
r8bpsplinb.f
1 c bpsplinb -- dmc 30 May 1996
2 c
3 c set up coefficients for bicubic spline with following BC's:
4 c * LHS and RHS BC under user control (see comments)
5 c * derivatives periodic in second coordinate (use pspline.for)
6 c
7  subroutine r8bpsplinb(x,inx,th,inth,fspl,inf3,
8  > ibcxmin,bcxmin,ibcxmax,bcxmax,
9  > wk,nwk,ilinx,ilinth,ier)
10 c
11  IMPLICIT NONE
12  INTEGER, PARAMETER :: R8=selected_real_kind(12,100)
13  integer inx,inth,inf3,nwk
14  real*8 x(inx),th(inth),fspl(4,4,inf3,inth),wk(nwk)
15  integer ibcxmin,ibcxmax
16  real*8 bcxmin(inth),bcxmax(inth)
17 c
18 c input:
19 c x(1...inx) -- abscissae, first dimension of data
20 c th(1...inth) -- abscissae, second (periodic) dimension of data
21 c fspl(1,1,1..inx,1..inth) -- function values
22 c inf3 -- fspl dimensioning, inf3.ge.inx required.
23 c
24 c boundary conditions input:
25 c ibcxmin -- indicator for boundary condition at x(1):
26 c bcxmin(...) -- boundary condition data
27 c =0 -- use "not a knot", bcxmin(...) ignored
28 c =1 -- match slope, specified at x(1),th(ith) by bcxmin(ith)
29 c =2 -- match 2nd derivative, specified at x(1),th(ith) by bcxmin(ith)
30 c =3 -- boundary condition is slope=0 (df/dx=0) at x(1), all th(j)
31 c =4 -- boundary condition is d2f/dx2=0 at x(1), all th(j)
32 c =5 -- match 1st derivative to 1st divided difference
33 c =6 -- match 2nd derivative to 2nd divided difference
34 c =7 -- match 3rd derivative to 3rd divided difference
35 c (for more detailed definition of BCs 5-7, see the
36 c comments of subroutine mkspline)
37 c NOTE bcxmin(...) referenced ONLY if ibcxmin=1 or ibcxmin=2
38 c
39 c ibcxmax -- indicator for boundary condition at x(nx):
40 c bcxmax(...) -- boundary condition data
41 c (interpolation as with ibcxmin, bcxmin)
42 c
43 c output:
44 c fspl(*,*,1..inx,1..inth) -- bicubic spline coeffs (4x4)
45 c ...fspl(1,1,*,*) is not replaced.
46 c
47  integer ilinx,ilinth,ier
48 c
49 c ilinx -- =1 on output if x(inx) pts are nearly evenly spaced (tol=1e-3)
50 c ilinth-- =1 on output if th(inth) evenly spaced (tol=1e-3)
51 c
52 c ier -- completion code, 0 for normal
53 c
54 c workspace:
55 c wk -- must be at least 5*max(inx,inth) large
56 c nwk -- size of workspace
57 c
58 c---------------------------------
59 c compute bicubic spline of 2d function, given values at the grid
60 c grid crossing points, f(1,1,i,j)=f(x(i),th(j)).
61 c
62 c on evaluation: for point x btw x(i) and x(i+1), dx=x-x(i)
63 c and th btw th(j) and th(j+1), dt=th-th(j),
64 c
65 c spline =
66 c f(1,1,i,j) + dx*f(2,1,i,j) + dx**2*f(3,1,i,j) + dx**3*f(4,1,i,j)
67 c +dt*(f(1,2,i,j) + dx*f(2,2,i,j) + dx**2*f(3,2,i,j) + dx**3*f(4,2,i,j))
68 c +d2*(f(1,3,i,j) + dx*f(2,3,i,j) + dx**2*f(3,3,i,j) + dx**3*f(4,3,i,j))
69 c +d3*(f(1,4,i,j) + dx*f(2,4,i,j) + dx**2*f(3,4,i,j) + dx**3*f(4,4,i,j))
70 c
71 c where d2=dt**2 and d3=dt**3.
72 c
73 c---------------------------------
74 c
75  integer ierx,ierth
76  REAL*8 zdum(1)
77 c
78 c---------------------------------
79 c
80  ier=0
81  if(nwk.lt.5*max(inx,inth)) then
82  write(6,'('' ?bpsplinb: workspace too small.'')')
83  ier=1
84  endif
85  if(inx.lt.2) then
86  write(6,'('' ?bpsplinb: at least 2 x points required.'')')
87  ier=1
88  endif
89  if(inth.lt.2) then
90  write(6,'('' ?bpsplinb: need at least 2 theta points.'')')
91  ier=1
92  endif
93 c
94  call ibc_ck(ibcxmin,'bcspline','xmin',0,7,ier)
95  call ibc_ck(ibcxmax,'bcspline','xmax',0,7,ier)
96 c
97 c check ilinx & x vector
98 c
99  call r8splinck(x,inx,ilinx,1.0e-3_r8,ierx)
100  if(ierx.ne.0) ier=2
101 c
102  if(ier.eq.2) then
103  write(6,'('' ?bpsplinb: x axis not strict ascending'')')
104  endif
105 c
106 c check ilinth & th vector
107 c
108  call r8splinck(th,inth,ilinth,1.0e-3_r8,ierth)
109  if(ierth.ne.0) ier=3
110 c
111  if(ier.eq.3) then
112  write(6,'('' ?bpsplinb: th axis not strict ascending'')')
113  endif
114 c
115  if(ier.ne.0) return
116 c
117 c------------------------------------
118  zdum=0.0_r8
119 c
120  call r8bcspline(x,inx,th,inth,fspl,inf3,
121  > ibcxmin,bcxmin,ibcxmax,bcxmax,
122  > -1,zdum,-1,zdum,
123  > wk,nwk,ilinx,ilinth,ier)
124 c
125  end