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