V3FIT
r8bpspline.f
1 c bpspline -- dmc 30 May 1996
2 c
3 c set up coefficients for bicubic spline with following BC's:
4 c * LHS and RHS BCs -- d3f/dx3 from 3rd divided difference, 1st coordinate
5 c * derivatives periodic in second coordinate
6 c
7 c see similar routine, bpsplinb, for more control over x boundary
8 c conditions...
9 c
10  subroutine r8bpspline(x,inx,th,inth,fspl,inf3,
11  > wk,nwk,ilinx,ilinth,ier)
12 c
13  IMPLICIT NONE
14  INTEGER, PARAMETER :: R8=selected_real_kind(12,100)
15  integer inx,inth,inf3,nwk
16  real*8 x(inx),th(inth),fspl(4,4,inf3,inth),wk(nwk)
17  integer ilinx,ilinth
18  integer ier
19 c
20 c input:
21 c x(1...inx) -- abscissae, first dimension of data
22 c th(1...inth) -- abscissae, second (periodic) dimension of data
23 c fspl(1,1,1..inx,1..inth) -- function values
24 c inf3 -- fspl dimensioning, inf3.ge.inx required.
25 c
26 c output:
27 c fspl(*,*,1..inx,1..inth) -- bicubic spline coeffs (4x4)
28 c ...fspl(1,1,*,*) is not replaced.
29 c
30 c ilinx -- =1 on output if x(inx) pts are nearly evenly spaced (tol=1e-3)
31 c ilinth-- =1 on output if th(inth) evenly spaced (tol=1e-3)
32 c
33 c ier -- completion code, 0 for normal
34 c
35 c workspace:
36 c wk -- must be at least 5*max(inx,inth) large
37 c nwk -- size of workspace
38 c
39 c---------------------------------
40 c compute bicubic spline of 2d function, given values at the grid
41 c grid crossing points, f(1,1,i,j)=f(x(i),th(j)).
42 c
43 c on evaluation: for point x btw x(i) and x(i+1), dx=x-x(i)
44 c and th btw th(j) and th(j+1), dt=th-th(j),
45 c
46 c spline =
47 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)
48 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))
49 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))
50 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))
51 c
52 c where d2=dt**2 and d3=dt**3.
53 c
54 c---------------------------------
55  integer ierx,ierth
56  real*8 zdum(1)
57 c---------------------------------
58 c
59  ier=0
60  if(nwk.lt.5*max(inx,inth)) then
61  write(6,'('' ?bpspline: workspace too small.'')')
62  ier=1
63  endif
64  if(inx.lt.2) then
65  write(6,'('' ?bpspline: at least 2 x points required.'')')
66  ier=1
67  endif
68  if(inth.lt.2) then
69  write(6,'('' ?bpspline: need at least 2 theta points.'')')
70  ier=1
71  endif
72 c
73 c check ilinx & x vector
74 c
75  call r8splinck(x,inx,ilinx,1.0e-3_r8,ierx)
76  if(ierx.ne.0) ier=2
77 c
78  if(ier.eq.2) then
79  write(6,'('' ?bpspline: x axis not strict ascending'')')
80  endif
81 c
82 c check ilinth & th vector
83 c
84  call r8splinck(th,inth,ilinth,1.0e-3_r8,ierth)
85  if(ierth.ne.0) ier=3
86 c
87  if(ier.eq.3) then
88  write(6,'('' ?bpspline: th axis not strict ascending'')')
89  endif
90 c
91  if(ier.ne.0) return
92 c
93 c------------------------------------
94 c
95  zdum=0.0_r8
96 c
97 c not-a-knot BCs for x, periodic for theta
98 c
99  call r8bcspline(x,inx,th,inth,fspl,inf3,
100  > 7,zdum,7,zdum,
101  > -1,zdum,-1,zdum,
102  > wk,nwk,ilinx,ilinth,ier)
103 c
104  return
105  end