Princeton Splines -- D. McCune Feb. 1999, updated Feb. 2004.
(software contributions also by W. Houlberg and A. Pletzer).
This document describes PSPLINE -- a library of Spline and Hermite
cubic interpolation routines for 1d, 2d, and 3d datasets on
rectilinear grids.
Abstract --
PSPLINE -- a collection of Spline and Hermite interpolation tools for
1D, 2D, and 3D datasets on rectilinear grids. Spline routines give
full control over boundary conditions -- the user may specify "periodic",
"not a knot", 1st derivative match, or 2nd derivative match boundary
conditions on either end of each grid dimension. Hermite routines take
as input the function value and derivatives at each grid point, giving
back a representation of the function between grid points. The splines
are continuously twice differentiable in all directions across all
grid cell boundaries and over the entire grid domain. Hermite functions
are continuously once differentiable in all directions over the entire
grid domain. For a representation of dimensionality N, an N-dimensional
Hermite or Spline function requires 2**N*(nx1*nx2*...*NxN) memory
words; there is also a somewhat faster but much more memory intensive
Spline representation requiring 4**N*(nx1*nx2*...*nxN) memory words.
Speed of evaluation of Hermite and Spline functions are similar--
both are constructed of cubic polynomial segments or "patches", which
are often much faster to evaluate than representations using Fourier
series or otherwise involving transcendental functions.
Both spline and hermite methods offer highly efficient options for
differentiable representations of higher dimensional objects-- see the
case study.
Many routines in pspline used to insist on a minimum gridsize of nx=4.
It is now allowed to use nx=2 (the smallest conceivable size, containing
only a single cubic or linear element). The test codes (e.g. pspltsub
in pspltest) have been extended with a new subroutine "smalltest" that
verifies handling of boundary condition controls in spline routines with
small grids (some special case coding was required for implementation).
A defect in a 1d spline fitting routine caused some spline fits using the
higher order boundary conditions
d2f/dx2 .ne. 0
d2f/dx2 from 2nd divided difference (edge parabolic fit)
d3f/dx3 from 3rd divided difference (edge cubic fit)
To be C0 instead of C2. This has been repaired.
The compact bicubic and tricubic splines are now constructed directly,
rather than by conversion from the non-compact representations. In the
case of the tricubics, this saves a factor of nearly 8 on memory usage
and a factor of about 2-5 on cpu time for construction of the spline
coefficients, depending on problem size, workstation cache size, etc.
I would like to acknowledge Kenneth Esler, a graduate student in
condensed matter physics at University of Illinois -- Urbana Champaign,
who downloaded pspline from the NTCC modules library.
In the course of his use of tricubic splines in his condensed matter
physics application, he encountered performance problems due to the
memory-intensive "non-compact" construction method used in the
previous version of the code. He sent a C++ code sample suggesting
a more compact construction method. I found his suggestion to be
usable; a fortran-90 implementation for compact construction of
bicubic and tricubic splines was developed and is now included in
the library. The fortran-90 interface (ezspline) uses this new method.
In the fortran-77 interface, use:
subroutine mkbicub -- construct compact bicubic C2 spline
subroutine mktricub -- construct compact tricubic C2 spline
These routines were completed on Feb. 26 2004 -- by Douglas McCune.
The advice of PPPL physicist Dr. Leonid Zakharov is much appreciated,
for pointing out the existence of the 2**N "compact" spline
representation. Dr. Zakharov and others also pointed out the
desirability of including cross derivative terms in the higher
dimensional Hermite formulations.
I also want to express thanks to ORNL physicist Dr. Wayne Houlberg
for reviewing the PSPLINE software and documentation, and for
writing "v_spline" which re-implemented the 1d spline tridiagonal
solver in a cleaner, more efficient form, with additional boundary
condition options.
Finally, I want to express appreciation to PPPL Computational
Physicist Dr. Alex Pletzer for his contribution of the EZspline
software interface to pspline.
PPPL Computational Plasma Physics Group, Feb. 5 1999 -- D. McCune
Updated: June 10, 1999
A 3d Spline Representation of B field greatly speeds stellarator fast
ion orbit calculations --
CPPG has developed two sets of highly general tricubic spline
interpolation tools for 3d objects. Periodic boundary condition
options make the 3d splines especially well suited to representation
of objects such as B(psi,theta,zeta) which are functions of poloidal
and toroidal angle parameters which arise in magnetic flux coordinate
systems. The spline tools have now been installed in R. White's
ORBIT code [1,2].
Two types of tricubic splines are supported:
* "explicit" spline-- the fastest, but very memory intensive:
64*Npsi*Ntheta*Nzeta memory words for the spline coefficients
* "compact" spline-- not quite as fast but eight times less
memory consumption: 8*Npsi*Ntheta*Nzeta memory words.
Evaluation of the 3d spline, which is continuously twice differentiable
in all coordinates, consists (for the explicit spline) in computing a
sum of 21 cubic polynomials-- regardless of the complexity of the object
represented (however, more complex objects require finer meshes and much
more memory). The fixed cost of evaluation distinguishes the spline
method from the Fourier harmonics summation representation typically
used in the past, which requires additional harmonics to represent
more complicated objects or fields, making each evaluation
proportionately more expensive.
In a trial case-- a stellarator equilibrium B field representation
requiring 304 Fourier harmonics to describe the poloidal and toroidal
variation-- replacement of the Fourier representation with a 3d spline
representation resulted in evaluation speed 80x faster on a DEC alpha.
By controlling the fineness of the mesh, the spline can reproduce the
Fourier representation to any desired accuracy. On the other hand,
the memory cost of the spline is significantly higher than the moments
representation, which needed just 4 x Nmoments x Npsi memory words.
The 3d spline has been vectorized and installed and tested in Roscoe
White's 3d Stellarator fast ion orbit code ORBIT3D on Killeen (a vector
Cray supercomputer at NERSC). For a reference problem (1000 particles
and 100 transits), the following table shows the speedups obtained,
relative to the older harmonics based code ORBITMN:
# of harmonics Speedup factor (ORBIT3D vs. ORBITMN)
======================================================
30 4 times faster
100 10 " "
300 50 " "
A typical 3d stellator equilibrium would require on the order of 100
harmonics (but 300 harmonic cases do arise). It is expected that
ORBIT3D, as a replacement to ORBITMN, will be heavily used to study
fast ion confinement in stellarator B field equilibria.
The 3d spline may have application in other codes as well -- e.g.
particle simulations. Because the data structure has locality in
three dimensions, it may be well suited to a 3d domain decomposition
for parallelization.
[June 10, 1999 addendum] ORBIT3D now allows switching between use
of the "explicit" and the "compact" spline method. Although lookups
using the explicit spline were 1.3 times faster than compact splines
on a standalone test application built on a DEC Alpha workstation,
inside ORBIT3D, using hand inlined vectorized interpolant evaluation
code for both methods, the net difference was only about 10%:
ORBIT3D on Cray-J90, 5000 ptcls, 100 transits
kspline=1 1227 sec (explicit spline)
kspline=2 1374 sec (compact spline) [3].
Probably for most 3d applications the compact spline's factor of 8
savings in memory utilization would outweigh the modest loss in
computational performance.
[1] R. B. White and M. S. Chance, Phys Fluids 27, 2455 (1984)
[2] R. B. White, Phys Fluids B 2, 845 (1990)
[3] R. B. White, priv. communication, June 10, 1999.
Cubic Splines are piecewise cubic functions that pass through a
set of given data points. Cubic splines are continuous and
continuously twice differentiable.
1D Cubic Splines have long been popular interpolation tools. Use of
splines always involves two steps: (1) setup of spline coefficients,
and (2) evaluation of spline at desired target points.
The setup of 1D spline coefficients involves a highly stable numerical
procedure-- solution of a diagonally dominant tri-diagonal matrix [1,2].
The following inputs are required:
the grid: x(i), i=1 to Nx x(i+1).gt.x(i) for all i.
the data: f(i), i=1 to Nx
a boundary condition at x(1)
a boundary condition at x(Nx).
This data uniquely specifies a cubic spline function, which is
output from the spline setup routine as a set of spline coefficients.
There are two choices for representation of the spline coefficients--
an "explicit" representation requiring 4*Nx memory words, and a Hermite-
like "compact" representation which stores just the data and 2nd
derivative at the grid points, requiring 2*Nx memory words.
The "explicit" representation is somewhat faster to use and is probably
the best choice for 1d splines. But, for higher dimensional splines,
these two representations require, respectively (where N is the
dimensionality or rank of the object):
explicit -- 4**N*(Nx1*Nx2*...*NxN)
compact -- 2**N*(Nx1*Nx2*...*NxN)
Thus for N=3, the explicit representation requires 8 times the memory
of the compact representation; this memory cost savings is likely to
outweigh the modestly higher computational cost of using the compact
representation scheme.
In 1d, the setup of a spline using the explicit representation yields
an array C(4,Nx) such that the spline function is evaluated by the
formulae
piecewise cubic spline function value:
s(x)=C(1,i)+(x-x(i))*C(2,i)+(x-x(i))**2*C(3,i)+(x-x(i))**3*C(4,i)
for x(i).le.x.le.x(i+1)
piecewise quadratic spline function 1st derivative:
s'(x)=C(2,i)+2*(x-x(i))*C(3,i)+3*(x-x(i))**2*C(4,i)
for x(i).le.x.le.x(i+1)
piecewise linear spline function 2nd derivative:
s''(x)=2*C(3,i)+6*(x-x(i))*C(4,i)
for x(i).le.x.le.x(i+1)
and satisfying the boundary conditions and the requirements:
s(x(i))=f(i), i = 1 to N
s'(x) is continuous accross all interior x(i)
s''(x) is continuous accross all interior x(i)
The compact representation requires an array D(2,Nx) where D(1,j)
gives the function value at grid point x(j), and D(2,j) gives the
interpolant's 2nd derivative at x(j). Then, evaluation takes the
following form.
If x(i).le.x.le.x(i+1), let
h = (x(i+1)-x(i))
p = (x-x(i))/h
q = 1-p
s(x)= q*D(1,i) + p*D(1,i+1) +
(h**2/6)*{[q**3-q]*D(2,i) + [p**3-p]*D(2,i+1)}
s'(x) = (D(1,i+1) - D(1,i))/h +
(h/6)*{[3*p**2-1]*D(2,i+1) - [3*q**2-1]*D(2,i)}
s''(x) = q*D(2,i) + p*D(2,i+1)
where (by the choice of D) the given data points are interpolated and
the specified boundary conditions are satisfied.
Typical boundary conditions given are s'(x(1)) or s''(x(Nx)). Many
spline routines also allow the boundary condition to be inferred from
the data in the vicinity of the boundary, e.g. match a numerical
third derivative from the first four data points as in "spline.for" [1]
or a "not a knot" boundary condition forcing the third derivative to be
continuous at x(2) and x(N-1) as in the default usage of "cubspl.for" [2].
(Note that the third derivative of a cubic spline is a step function
and generally not continuous -- imposition of continuity near the
boundary constitutes the "not a knot" boundary condition).
Another useful option is the "periodic" boundary condition
s'(x(1))=s'(x(Nx))
s''(x(1))=s''(x(Nx)).
(note s(x(1))=s(x(Nx)), i.e. f(1)=f(Nx), is not required).
This requires a simple modification to the spline coefficients
algorithm (typically left as an exercise in the texts) as the
coefficients solution matrix is no longer strictly tridiagonal
(it has non-zero off diagonal corner elements). The PSPLINE
library originally took its name from "pspline.for" which contains
an algorithm for producing a periodic spline.
v_spline.for gathers together all the 1d spline coefficient calculations
into a single routine.
references --
[1] Forsythe, Malcolm, and Moler, "Computer Methods for Mathematical
Computations", Prentice-Hall, 1977. (see chapter on interpolation).
[2] C. De Boor, "A Practical Guide to Splines".
[3] W. Press, S. Teukolsky, W. Vetterling, B. Flannery, "Numerical
Recipies in Fortran 77 (2nd edition)", pp 107-110 -- discusses the
compact representation.
Also useful:
[4] Engeln-Muellges, Uhlig, Numerical Algorithms with Fortran, Springer,
1996, p.251
It turns out that 1D splines can be extended to define interpolation
of data given on a 2d rectilinear grid:
f(i,j), i=1 to Nx, j=1 to Ny,
vs. x(i), i=1 to Nx, x(i).lt.x(i+1), and
vs. y(j), j=1 to Ny, y(j).lt.y(j+1).
As in the 1d case, boundary condition data is required: typically,
df/dx or d2f/dx2 vs. y(j) at x(1) and x(Nx), and, df/dy or d2f/dy2
vs. x(i) at y(1) and Y(Ny). As in the 1d case, sometimes the data
itself can be used to infer a reasonable boundary condition, obviating
the requirement for its explicit specification. Or, perhaps, a periodic
boundary condition can be used. In the case of a bicubic spline
representation of a function f(r,theta) on a polar coordinate system,
the periodic boundary condition would be natural, and the preferred
choice.
For a given choice of Nx*Ny data points f(i,j) and 2*Nx + 2*Ny
boundary conditions, it turns out there is exactly one bicubic
spline function which is continuously twice differentiable over
the entire domain x(1) .le. x .le. x(Nx) & y(1) .le. y .le. y(Ny).
Existence / Uniqueness proofs and a construction method will be
sketched below.
In the PSPLINE library, a choice between two simple representations
is available for bicubic splines. In the "explicit" representation,
for each data point (x(i),y(j)), a 4x4 array of 16 coefficients Ckl
are retained.
The bicubic spline function s(x,y) is then evaluated as follows:
for x(i).le.x.le.x(i+1), let dx=(x-x(i))
for y(j).le.y.le.y(j+1), let dy=(y-y(j))
let C11...C44 denote the 16 coeffcients associated with (i,j).
then:
s(x,y)= C11 + dx*C21 + dx**2*C31 + dx**3*C41
+dy*(C12 + dx*C22 + dx**2*C32 + dx**3*C42)
+dy**2*(C13 + dx*C23 + dx**2*C33 + dx**3*C43)
+dy**3*(C14 + dx*C24 + dx**2*C34 + dx**3*C44)
This can be thought of as the evaluation of 5 cubic pieces, each
of which can be cast in the form A+dx*(B+dx*(C+dx*D)), i.e. three
A*x+B operations: so then a total of 15 A*x+B operations for
a single evaluation of s(x,y)-- thus, a very computationally
efficient representation.
The formulae for ds/dx, ds/dy, d2s/dx2, d2s/dy2, d2s/dxdy are
all derived from the formula for s(x,y) in the obvious way.
The explicit representation requires 16*Nx*Ny memory words. The
compact representation requires only 4*Nx*Ny memory words, and
is only slightly slower computationally. In the compact
representation one uses
D(1,i,j) = the function value at (x(i),y(j))
D(2,i,j) = the 2nd derivative d2s/dx2 at (x(i),y(j))
D(3,i,j) = the 2nd derivative d2s/dy2 at (x(i),y(j))
D(4,i,j) = the 4th derivative d4s/dx2dy2 at (x(i),y(j))
where the spline derivative coefficients are carefully chosen
to achieve the required level of continuity, differentiability,
and satisfaction of boundary conditions.
Evaluation of the spline using the compact representation is
as follows:
for x(i).le.x.le.x(i+1) and y(j).le.y.le.y(j+1), let
hx = x(i+1)-x(i) hy = y(j+1)-y(j)
px = (x-x(i))/hx py = (y-y(j))/hy
qx = 1-px qy = 1-py
rx = (px**3-px) ry = (py**3-py)
sx = (qx**3-qx) sy = (qy**3-qy)
s(x,y) = qx*(qy*D(1,i,j)+py*D(1,i,j+1)) +
px*(qy*D(1,i+1,j)+py*D(1,i+1,j+1)) +
(hx**2/6)*{sx*(qy*D(2,i,j)+py*D(2,i,j+1)) +
rx*(qy*D(2,i+1,j)+py*D(2,i+1,j+1))} +
(hy**2/6)*{qx*(sy*D(3,i,j)+ry*D(3,i,j+1)) +
px*(sy*D(3,i+1,j)+ry*D(3,i+1,j+1))} +
(hx**2*hy**2/36)*
{sx*(sy*D(4,i,j)+ry*D(4,i,j+1)) +
rx*(sy*D(4,i+1,j)+ry*D(4,i+1,j+1))}
Of course the key is to evaluate the set of coefficients {Ckl}
or {D} such that the continuously twice differentiable property is
achieved. The fact that this can be done, and how it is done,
is considered in the subtopics. The subtopics focus on obtaining
{Ckl}; the derivation of {D} from {Ckl} is straight-forward.
A proof is sketched here for the following result:
given: an x grid {x(1),x(2),...,x(Nx)}, x(i).lt.x(i+1),
a y grid {y(1),y(2),...,y(Ny)}, y(j).lt.y(j+1)
a data array f(i,j) specifying the desired
interpolation values,
for each y grid point j:
a boundary condition bcx1(j) at x(1) and
bcx2(j) at x(Nx),
and for each x grid point i:
a boundary contition bcy1(i) at y(1) and
bcy2(i) y(Ny),
then: there exists one and only one bicubic spline
function s(x,y) satisifying
s(x(i),y(j)) = f(i,j)
AND all the boundary conditions
AND s(x,y) continuous in [x(1),x(Nx)]x[y(1),y(Nx)]
AND ds/dx continuous in [x(1),x(Nx)]x[y(1),y(Nx)]
AND ds/dy continuous in [x(1),x(Nx)]x[y(1),y(Nx)]
AND d2s/dx2 continuous in [x(1),x(Nx)]x[y(1),y(Nx)]
AND d2s/dy2 continuous in [x(1),x(Nx)]x[y(1),y(Nx)]
AND d2s/dxdy continuous in [x(1),x(Nx)]x[y(1),y(Nx)].
Theoretical results on splines can be had by persuing methods of
classical analysis -- such as, construction of set of independent
"basis" splines from which any spline can be built. The "basis"
splines discussed here are not efficient for computation, but are
useful as a theoretical tool.
It has been noted that boundary conditions for splines can be
expressed in a variety of ways. For simplicity of this discussion
(and without loss of generality) it will be assumed that the
boundary conditions at x(1) and x(Nx) specify ds/dx at each y(j),
and that the boundary conditions at y(1) and y(Ny) specify ds/dy
at each x(i).
Then, the basis set can be formed as follows. First define the
the 1d "x basis" splines (for i=1 to Nx) s[i](x) on
{x(1),x(2),...,x(Nx)} satisfying
s[i](x(j)) = delta(i,j) (i.e. 1 if i=j, 0 otherwise)
and
ds[i]/dx = 0 at x(1)
and
ds[i]/dx = 0 at x(Nx).
Because of the properties of 1d splines, s[i] exists and is
unique. Then define also s[0](x) satisfying
s[0](x(j)) = 0 for all j
and
ds[0]/dx = 1 at x(1)
and
ds[0]/dx = 0 at x(Nx).
Similarly, define s[Nx+1](x) satisfying
s[Nx+1](x(j)) = 0 for all j
and
ds[Nx+1]/dx = 0 at x(1)
and
ds[Nx+1]/dx = 1 at x(Nx+1).
Again the properties of 1d splines show that s[0](x) and s[Nx+1](x)
exist and are unique.
Similarly, define a set of Ny+2 1d "y basis" splines s[j](y) on
{y(1),y(2),...,y(Ny)} satisfying the precisely analogous conditions
on the values of s[j](y(k)) and ds[j]/dy at y(1) and y(Ny). Again,
all of the 1d "y basis" splines exist and are unique.
From the 1d basis splines construct 2d basis splines s[i,j](x,y) by
the simple formula
s[i,j](x,y) = s[i](x) * s[j](y).
It can readily be seen that s[i,j](x(k),y(l)) = delta(i,k)*delta(j,l)
(for the interior), and ds/dx[0,j](x(1)) = delta(j,l) for the x(1)
boundary, and so on for the x(Nx), y(1) and y(Ny) boundaries.
Because s[i](x) and s[j](y) are continuously twice differentiable
in x and y respectively, it is clear that s[i,j] is continuously
twice differentiable on [x(1),x(Nx)]x[y(1),y(Nx)].
Therefore, the sum
s(x,y) = [sum i=1 to Nx, j=1 to Ny] f(i,j)*s[i,j](x,y)
+ [sum j=1 to Ny] bcx1(j)*s[0,j](x,y)
+ [sum j=1 to Ny] bcx2(j)*s[Nx+1,j](x,y)
+ [sum i=1 to Nx] bcy1(i)*s[i,0](x,y)
+ [sum i=1 to Nx] bcy2(i)*s[i,Ny+1](x,y)
is also continuously twice differentiable on
[x(1),x(Nx)]x[y(1),y(Nx)]. But by the definition of the s[i,j] basis
splines, s(x,y) also satisfies all of the required conditions
s(x(i),y(j)) = f(i,j) (i=1 to Nx, j=1 to Ny)
ds/dx(x(1),y(j)) = bcx1(j)
ds/dx(x(Nx),y(j)) = bcx2(j)
ds/dy(x(i),y(1)) = bcy2(i)
ds/dy(x(i),y(Ny)) = bcy2(i)
This proves existence.
Uniqueness is shown in the usual way, using a proof by contradiction.
Assume there exists another bicubic spline s2(x,y) that satisfies all
the conditions. Then the difference
r(x,y) = s(x,y) - s2(x,y)
must satisfy r(x(i),y(j)) = 0 for all i,j, and also dr/dx = 0 at
x(1) and x(Nx), and dr/dy = 0 at y(1) and y(Ny).
The general bicubic spline r(x,y) can be expressed in the form
r(x,y)= C11 + dx*C21 + dx**2*C31 + dx**3*C41
+dy*(C12 + dx*C22 + dx**2*C32 + dx**3*C42)
+dy**2*(C13 + dx*C23 + dx**2*C33 + dx**3*C43)
+dy**3*(C14 + dx*C24 + dx**2*C34 + dx**3*C44)
where each Ckl is associated with a subregion
[x(k),x(k+1)]x[x(l),x(l+1)] of the gridded space.
Consideration of the restriction of r(x,y) to x(i) or y(j)
quickly leads to the conclusion that C11=C21=C31=C41=C12=C13=C14=0
(because a 1d spline with zero values and zero slope boundary
conditions is itself identically zero).
This leaves in each grid square [x(k),x(k+1)]x[x(l),x(l+1)] a set
of 9 coefficients C22,C32,C42,C23,C33,C43,C24,C34,C44 not yet
demonstrated to be zero. But, consideration of continuity and
differentiability requirements across grid square boundaries
quickly leads to an overdetermined system of linear equations
involving these coefficients, the only solution of which is
that all these coefficients, too, are zero.
Thus r(x,y)=0 everywhere, and s(x,y)=s2(x,y) -- contradicting
the assumption of non-uniqueness.
Thus uniqueness is proven.
Given that there exists a unique continuous twice differentiable
bicubic spline interpolating given data and satisfying given
boundary conditions, how is it to be constructed efficiently?
The bicubic spline existence proof showed a method of construction,
but, it is not practical. Each spline basis function is an object
with size of proportional to (Nx+Ny), and there are (Nx+2)*(Ny+2)
such 2d basis splines to be considered. Any algorithm using basis
splines is at best going to have computational performance scaling
of Nx*Ny*(Nx+Ny). An construction algorithm that scales as Nx*Ny
should be sought first. It can be found.
Before describing this in detail, a digression on boundary condition
"sets" is necessary.
The set of boundary conditions on the spline at x(1) and x(Nx) shall
be denoted as the X-set, and the set at y(1) and y(Ny) shall be denoted
the Y-set.
Boundary condition sets can be divided into two categories:
"homogeneous" and "inhomogeneous". A homogeneous boundary condition
set A has the property that if s1(x,y) satisfies A, then (for any real
number c) c*S1 also satisfies A, and, if both s1(x,y) and s2(x,y)
satisfy A then so also does s1+s2.
The following are examples of homogeneous boundary condition X-sets:
* periodic boundary condition:
ds/dx(x(1),y)=ds/dx(x(Nx),y), d2s/dx2(x(1),y)=d2s/dx2(x(Nx),y)
* any set constructed by taking one condition from column A
and one from column B:
A B
========================== =============================
d3s/dx3(x(2),y) continuous d3s/dx3(x(Nx-1),y) continuous[*]
ds/dx(x(1),y)=0 ds/dx(x(Nx),y)=0
d2s/dx2(x(1),y)=0 d2s/dx2(x(Nx),y)=0
d[k]s/dx[n] = k'th divided d[k]s/dx[k] = k'th divided
difference at x(1) [**] difference at x(n) [**]
[*] these are the "not a knot" boundary conditions, which are
homogeneous.
[**] e.g. for the "first" divided difference at x(1):
ds/dx(x(1),y(j)) = [f(2,j)-f(1,j)]/(x(2)-x(1))
2nd and 3rd divided difference expressions for d2s/dx2 and
d3s/dx3 are also supported (software contributed by Wayne
Houlberg).
* inhomogeneous boundary conditions all involve explicitly setting
a derivative or 2nd derivative to a fixed non-zero value at a
boundary point.
Homogeneous boundary conditions are "nice" because they facilitate
the construction of bicubic splines. To see this, consider a
particular bicubic spline construction problem (non-compact
representation). We are seeking a continuous twice differentiable
spline of the form
s(x,y)= C11 + dx*C21 + dx**2*C31 + dx**3*C41
+dy*(C12 + dx*C22 + dx**2*C32 + dx**3*C42)
+dy**2*(C13 + dx*C23 + dx**2*C33 + dx**3*C43)
+dy**3*(C14 + dx*C24 + dx**2*C34 + dx**3*C44)
where the 16 coefficients {Ckl} are defined separately in each
grid square [x(i),x(i+1)]x[y(j),y(j+1)].
and we are given:
* an x grid {x(1),x(2),...,x(Nx)}, x(i).lt.x(i+1),
* a y grid {y(1),y(2),...,y(Ny)}, y(j).lt.y(j+1)
* a data array f(i,j) specifying the desired
interpolation values,
* The X-set of boundary conditions:
for each y grid point j:
a boundary condition bcx1(j) at x(1) and
bcx2(j) at x(Nx),
* The Y-set of boundary conditions:
for each x grid point i:
a boundary contition bcy1(i) at y(1) and
bcy2(i) y(Ny).
which we can denote as problem S.
If we assume that the Y-set is homogeneous, we can construct the
unique twice differentiable bicubic spline by the following
algorithm:
Algorithm A:
(a) at each y(j), compute the 1d spline through the data points
f(x(1),y(j)), f(x(2),y(j)), ... ,f(x(Nx),y(j)) and satisfying
the (not necessarily linear) boundary conditions at (x(1),y(j))
and (x(Nx),y(j)). In each grid square this yields the
coefficients {C11, C21, C31, C41}. Note that for grid cell
(i,j), C11 is just f(i,j).
(b) using the HOMOGENEOUS Y-set boundary condition at each x(i),
find the 1d splines s1(y), s2(y), s3(y) and s4(y) satisfying:
s1(x(i),y(j))=C11(i,j)
s2(x(i),y(j))=C21(i,j)
s3(x(i),y(j))=C31(i,j)
s4(x(i),y(j))=C41(i,j)
which yields four spline coefficients for each 1d spline.
These 4 x 4 coefficients are precisely the {Ckl} sought.
Because the time for computation of an N point 1d spline, (which
is just the solution of a tridiagonal system) is proportional to c*N,
the above algorithm scales 5*c*Nx*Ny -- i.e. proportional to Nx*Ny,
as desired.
Why does this method work? It can be seen by thinking (again) in
terms of basis splines. Consider a simplified problem P[j0] for which
(1) the X boundary conditions and function values are zero for
all y(j) except when j=j0, and
(2) The Y-set of boundary conditions is linear.
Let s[j0](x) be the 1d spline at y(j0) that satisfies the given data
and boundary conditions along along y(j0). Then, there is a 1d basis
spline b(y) which satisfies the lineary Y-set boundary conditions, and,
b(y(j)) = 1 if j=j0
= 0 otherwise.
The solution to P[j0] is just
s(x,y) = s[j0](x)*b(y),
and this solution IS CONSTRUCTED by algorithm A, because all of the
1d y spline problems presented yield solutions proportional to b(y),
and because of the Y-set boundary condition linearity allows these
solutions to be combined with preservation of the boundary condition.
The Y-set linearity also allows the full problem S to be considered
as a sum of P-type problems. It is easily shown that
sum[j0=1 to Ny] A(P[j0]) = A(sum[j0=1 to Ny] P(j0)) = A(S).
=================
Thus, the general problem of a construction algorithm with Nx*Ny
performance scaling is solved, when granted a HOMOGENEOUS Y-set of
boundary conditions.
If we have a problem S' where the Y-set is inhomogeneous, we proceed as
follows (algorithm B):
Let A' be the "complementary" solution algorithm: it solves the
problem with a linear X-set by splining first in y, then splining
the resulting y 1d spline coefficients in x.
Algorithm B involves a combination of algorithms A and A'. To solve
problem S using B:
(a) use algorithm A to solve the related problem S* for a bicubic
spline that satisfies the X-set boundary conditions, and
"not-a-knot" homogeneous boundary condition at y(1) and y(Ny).
(b) use algorithm A' to solve the residual problem R=S'-S*. Note
that the data values and X-set boundary conditions for R are
all zero. Problem R's Y-set boundary conditions are the
difference between the S* boundary condition obtained, and the
desired (inhomogeneous) boundary condition of the original problem.
(c) The sum of the solutions A(S*) and A'(R) is the solution to
the original problem, S'.
The computational cost of this algorithm-- approximately 1.5x the
cost of algorithm A alone-- and scales like (Nx*Ny), and is
acceptable.
The spline coefficients in the non-compact representation are
closely related to the coefficients of the compact representation,
i.e. for any grid cell in a bicubic spline,
C(1,1) = f
2*C(3,1) = fxx
2*C(1,3) = fyy
4*C(3,3) = fxxyy
Consideration of these relations allows a compact bicubic (or
tricubic) spline construction method to be written in a manner
analogous to the non-compact method (but considerably more
efficient: C(2,2), etc. are not re-splined). The issue of
inhomogeneous boundary conditions can be handled in the same
way as it is for the non-compact construction method.
NAG comparison notes -- D. McCune 21 Feb. 1999 --
Summary: NAG bicubic B-splines are about 4x slower to evaluate,
but, 16x more compact in memory consumption then the "explicit"
representation, and 4x more compact than the "compact"
representation used in PSPLINE.
Tests comparing PSPLINE r8bcspeval and NAG e02def indicate that
r8bcspeval is 4x faster. The B-spline representation created in NAG
e01daf and evaluated in NAG e02def is far more compact -- consuming
sixteen times less memory. But, there are no options for controlling
boundary conditions in the nag bicubic b- spline routine. Also, there
are no 3d routines.
The comments about the extensibility of 1d splines to 2d applies equally
well to 3d rectilinear grids. So, too, do the comments about the method
of construction, the differentiability properties of the result, and
the choice between the (faster) "explicit" representation requiring
64*Nx*Ny*Nz memory words, and the "compact" representation requiring
8*Nx*Ny*Nz memory words.
The input data for setup of a 3d spline are
f(i,j,k) at x(i),i=1 to Nx, y(j), j=1 to Ny, z(k), i=1 to Nz.
with x(i+1).gt.x(i),
y(j+1).gt.y(j), and
z(k+1).gt.z(k) required.
General boundary conditions are 1st or second derivative information
across a grid boundary plane, e.g. df/dx at (x(1),y(j),z(k)), an Ny x Nz
array of data.
Periodic and "not a knot" boundary conditions are again available. For
data given on a toroidal coordinate system, f(r,theta,phi), periodic
boundary conditions in the theta and phi directions would be the
preferred choice.
The storage requirement for the explicit 3d spline representation
data set is 64*Nx*Ny*Nz. The cost of evaluation is 21 computations
of the form ((A*x+B)*x+C)*x+D (63 A*x+B type operations). The storage
requirement for the compact representation is 8*Nx*Ny*Nz, and timing
tests (on a DEC alpha workstation) show the cost of evaluation to be
approximately 1.4 times greater than for the explicit representation.
The PSPLINE library offers three routines, with different boundary
condition option sets, for the setting up of tricubic spline data sets.
In addition, there is a scalar bicubic spline evaluation routine and
a vectorized spline evaluation routine.
Remarkably, standard mathematics software libraries (NAG, netlib) do not
seem to include tricubic splines. Yet this representation can be highly
efficient in practice -- the author has used a tricubic spline to replace
an optimized hybrid Fourier series & 1d spline representation of an
actual physics model object, with a speedup in evaluation of a factor
of nearly 100 on a DEC alpha workstation (although a substantial memory
cost is incurred). See the "Case Study" (top level subtopic).
Certainly 4d and even 5d multi-cubic splines could be constructed
using the compact representation. However, the storage and evaluation
costs (algebra) are exponential in the dimensionality.
Higher-dimensional splines (beyond 3d) are not currently available in
the library.
Hermite cubic interpolation is sometimes used as an alternative to
splines. The Hermite method avoids the global tridiagonal matrix
inversion for coefficients evaluation required by splines, and the
local behaviour of the Hermite interpolant is controlled strictly
by local data -- sometimes an advantage. However, the Hermite
interpolant is continuously differentiable only once.
Hermite interpolation requires function values and derivatives
to be provided by the user at all grid points. The following
inputs are required:
the grid: x(i), i=1 to N x(i+1).gt.x(i) for all i.
the data: f(0,i), i=1 to N f(0,i)=f(x(i))
the derivative: f(1,i), i=1 to N, f(1,i) = [df/dx](x(i)).
In 1d, the construction of the Hermite interpolant is extremely simple.
In each interval [x(i),x(i+1)], there is only one cubic polynomial
that matches the data values f(0,i), f(0,i+1) and the derivative values
f(1,i), f(1,i+1) at the end points. This cubic polynomial is the
Hermite cubic for that interval. The value and derivative of the
Hermite cubic is continuous with the Hermite cubics in the neighbouring
intervals simply by virtue of the sharing the data value and derivative
information {f(0,i), f(1,i)} and {f(0,i+1),f(1,i+1)}.
In practice, it is not necessary to solve for the coefficients of the
cubic piece in each interval. Rather, it can be constructed on the
fly from the Hermite basis functions, which, on the unit interval [0,1]
are given by:
alpha(x)=x*x*(-2*x + 3)
alpha(0)=0 alpha(1)=1 alpha'(0)=0 alpha'(1)=0
alphabar(x)=1-alpha(x)
alphabar(0)=1 alphabar(1)=0 alphabar'(0)=0 alphabar'(1)=0
beta(x)=x*x*(x-1)
beta(0)=0 beta(1)=0 beta'(0)=0 beta'(1)=1
betabar(x)=x*(x*x-2*x+1)
betabar(0)=0 betabar(1)=0 betabar'(0)=1 betabar'(1)=0
Speed of evaluation of values and derivatives on a general grid requires
access to the information hx(i)=(x(i+1)-x(i)) and hxi(i)=1/hx(i). Given
this information, Hermite evaluation speeds are similar to spline
evaluation speeds.
The memory consumption of Hermite and Spline single and multivariate
representations are equivalent, as per the following table:
1d 2d 3d M-dimensional
Hermite 2*N 4*Nx*Ny 8*Nx*Ny*Nz 2**M*(Nx1*Nx2*...*NxM)
or
"compact" spline
The behaviour of the bicubic Hermite interpolant is analogous to the 1d
case. The required input information is the function value and
derivatives at each grid point:
x(i), i = 1 to Nx, x(i).lt.x(i+1) -- the x grid
y(j), j = 1 to Ny, y(j).lt.y(j+1) -- the y grid
f(0,i,j) = f(x(i),y(j))
f(1,i,j) = [df/dx](x(i),y(j))
f(2,i,j) = [df/dy](x(i),y(j))
f(3,i,j) = [d2f/dxdy](x(i),y(j))
Bicubic Hermite evaluation at (x,y) is done by direct construction,
as follows. Suppose that (x,y) is in the grid cell whose lower left
corner is (x(i),y(j)).
Let hx(i)=(x(i+1)-x(i)),
hxi(i)=1/hx(i),
hy(j)=(y(j+1)-y(j)),
hyi(j)=1/hy(j).
x~=(x-x(i))*hxi(i)
y~=(y-y(j))*hyi(j)
The normalized parameters (x~,y~) describe the position of (x,y) within
the grid cell.
The Hermite basis functions alpha, alphabar, beta, and betabar
(defined in the previous section) are then used:
Let
ax=alpha(x~)
axbar=alphabar(x~)
bx=beta(x~)
bxbar=betabar(x~)
ay=alpha(y~)
aybar=alphabar(y~)
by=beta(y~)
bybar=betabar(y~)
then
h(x,y) = aybar*(axbar*f(0,i,j)+ax*f(0,i+1,j))
+ ay*(axbar*f(0,i,j+1)+ax*f(0,i+1,j+1))
+ hx(i)*(
aybar*(bxbar*f(1,i,j)+bx*f(1,i+1,j))
+ ay*(bxbar*f(1,i,j+1)+bx*f(1,i+1,j+1)) )
+ hy(j)*(
bybar*(axbar*f(2,i,j)+ax*f(2,i+1,j))
+ by*(axbar*f(2,i,j+1)+ax*f(2,i+1,j+1)) )
+ hx(i)*hy(j)*
bybar*(bxbar*f(3,i,j)+bx*f(3,i+1,j))
+ by*(bxbar*f(3,i,j+1)+bx*f(3,i+1,j+1)) )
with similar expressions, involving derivatives axbar', ax',
bxbar', and bx', for d[h(x,y)]/dx, and involving derivatives aybar',
ay', bybar', and by', for d[h(x,y)]/dy, and all polynomial
derivatives for d2[h(x,y)]/dxdy. Full details can be seen
in the PSPLINE library source code: herm2ev.for.
From the definitions of the basis functions and their it follows
that all the necessary continuity and differentiability conditions
are met. For example, h(x(i),y(j))=f(0,i,j) (as required) because
x(i) -> x~=0, y(i) -> y~=0, => axbar=aybar=1, ax=ay=0, and
bx=bxbar=by=bybar=0.
Similarly, d[h(x(i),y(j+1)]/dy = f(2,i,j+1) (as required) because
x~=0, y~=1 and the known simple behaviour of the basis functions and
their derivatives.
Inside each cell, h(x,y) is locally a bicubic polynomial patch, and
is multiply continuously differentiable.
Along any cell boundary, the formulae for h(x,y) and its 1st derivatives
reduce to expressions which depend only on the data at the endpoints
of the line segment defining the boundary. These limiting expressions
are identical for cells on either side of the boundary segment. This is
what yields the continuity of h(x,y) and all its 1st derivatives across
these boundaries and hence over the entire gridded region.
The PSPLINE library provides routines for evaluation of 2d Hermite
interpolants and their first derivatives -- see herm2ev.for.
Hermite interpolation is different than spline interpolation in
that (a) no global matrix solution is required for evaluation of
coefficients, and (b) continuity properties are maintained for ANY
choice of derivatives at the grid points.
Thus, if derivative information is not available analytically, it
may be provided by other suitable (i.e. numerical) means.
The PSPLINE library Hermite software includes a choice of routines
for numerical evaluation of derivatives-- see the section on software.
The construction method using 1d Hermite cubic basis functions, described
for 2d bicubic Hermite, generalizes in a straightforward manner to the 3d
case as well. The details will not be spelled out here -- they can be
seen in the PSPLINE library source code: herm3ev.for.
The input data needed for tricubic Hermite evaluation are:
x(i), i = 1 to Nx, x(i).lt.x(i+1) -- the x grid
y(j), j = 1 to Ny, y(j).lt.y(j+1) -- the y grid
z(k), k = 1 to Nz, z(k).lt.z(k+1) -- the y grid
f(0,i,j,k) = f(x(i),y(j),z(k))
f(1,i,j,k) = [df/dx](x(i),y(j),z(k))
f(2,i,j,k) = [df/dy](x(i),y(j),z(k))
f(3,i,j,k) = [df/dz](x(i),y(j),z(k))
f(4,i,j,k) = [d2f/dxdy](x(i),y(j),z(k))
f(5,i,j,k) = [d2f/dxdz](x(i),y(j),z(k))
f(6,i,j,k) = [d2f/dydz](x(i),y(j),z(k))
f(7,i,j,k) = [d3f/dxdydz](x(i),y(j),z(k))
The PSPLINE library provides routines for evaluation of 3d Hermite
interpolants and their first derivatives -- see herm3ev.for.
Hermite interpolation is different than spline interpolation in
that (a) no global matrix solution is required for evaluation of
coefficients, and (b) continuity properties are maintained for ANY
choice of derivatives at the grid points.
Thus, if derivative information is not available analytically, it
may be provided by other suitable (i.e. numerical) means.
The PSPLINE library Hermite software includes a choice of routines
for numerical evaluation of derivatives-- see the section on software.
Interpolation software often represents a performance critical
component of real applications. Here, in decreasing order of
importance, are some tips for obtaining the best possible
performance from pspline:
* where possible, do "vectorized" spline evaluations
(including zone lookup), i.e. use routines which carry out
multiple spline evaluations in a single call.
"Zone lookup" plays a critical role in the performance of
interpolation systems.
* where possible, use evenly spaced grids. Otherwise...
* using custom zone lookup software may help when special
properties of the grid are known in advance. Otherwise...
* use "(r8)genxpkg" to create a grid data object for each
grid dimension, and use the corresponding "genxpkg-
enabled" set of evaluation routines.
-> On all workstations tested so far, genxpkg algorithm #3,
piecewise-linear indexing function, is the fastest, out-
performing traditional binary search methods by almost a
factor of two.
-> if the target point locations for interpolation are
slowly varying, selecting the genxpkg option to "use
previous result" can be a win (~30% performance gain
on typical workstations). On the other hand, if the
target point locations vary rapidly, this option will
exact a small penalty (~5% performance loss on typical
workstations). Caution: in an evaluation loop, "use
previous result" creates a dependency on the previous
loop iteration, which may break vectorization and
exact a heavy performance penalty, on some machines.
There is a trade-off between memory and speed in this software...
* use "explicit" splines rather than "compact" splines, when
the memory cost is not too high. Prefer "explicit" splines
to Hermite methods, unless the spline setup cost is a
significant factor, or data noise is a consideration.
The following step will involve too much labor, except where the
interpolation performance is extraordinarily critical
* for the ultimate in performance, hand inline the spline
evaluation computation into your application (but note this
gets tricky, especially at higher dimensionalities)!
When moving to a finer grid, C2 cubic splines will converge much
more rapidly than C1 Hermite to an accurate reproduction of the
underlying model function, if the function is smooth and the
representative data points chosen sample the model function
with machine precision accuracy.
But: when fitting "noisy experimental data", splines can cause
amplification and spreading of the noise, called "ringing". The
same comment applies if the underlying model function is not
really C2 on the resolution of the given numerical grid. If
noisy data is expected, Hermite methods should be given preference.
The Akima Hermite method, in particular, has well demonstrated
variation minimizing properties.
[Hiroshi Akima, Communications of the ACM, Jan 1974, Vol. 17 No. 1]
Most users should follow this strategy:
* use EZspline
But, if the f77 interface must be used, then,
* use genxpkg to control zone lookup algorithms and precompute
grid-specific data which is useful to the lookup/evaluation
software, and, use vectorized spline lookup/evaluation routines,
or...
* provide your own vectorized zone lookup algorithm, and use the
vectorized spline evaluation routines.
fortran-90 interface with good performance and ease of use features
by Alex Pletzer...
for details, see the separate EZspline help documents. Links to
these documents can be found at the pspline web page:
https://w3.pppl.gov/NTCC/PSPLINE/pspline.html
The subtopics describe the pspline "standard routines" for Spline
and Hermite interpolation. See also the "EZspline_Software"
section for the description of Alex Pletzer's f90 interface to
pspline, which will probably be easier to use for new code.
Note: the routines described are coded for single precision
(floating point arguments and arrays are declared REAL ...),
but...
The PSPLINE library also contains a complete set of double
precision routines (floating point arguments and arrays are
declared REAL*8 ...). These routines have the same names as
the single precision routines except that the prefix "R8"
prepended.
For example:
standard name REAL*8 name
one-d spline setup routine CSPLINE R8CSPLINE
two-d spline setup routine BCSPLINE R8BCSPLINE
The declaration "REAL*8" instead of "DOUBLE PRECISION" is used
for reasons of portability to typical Cray environments.
The pspline software includes a "public" module with a definition
of all the (non-EZspline) pspline subroutine calls. Fortran-90
users can add the line
use pspline_calls
to their source, for compile-time checking of argument lists
of pspline subroutine calls.
Caution: as with any f90 interface module, some calls that are
functionally correct may be flagged. This happens e.g. if an
element of an array is passed to a dummy argument as an "address"
indicating that that array element and subsequent elements in the
caller's array contain the necessary data. Unfortunately, f90
sees this as an attempt to pass a scalar where an array is expected.
Note that the use of the pspline_calls modules is entirely optional:
it is essentially a debugging aid.
This routine takes a grid array x(1:nx) and prepares a dataset
xpkg(1:nx,4) imbedding precomputed data specific to the x
axis, and option specifications for pspline zone lookup software.
Usage:
real x(nx) ! the grid (in) -- strict ascending
real xpkg(nx,4) ! the "package" (out)
integer iper ! periodicity flag (in)
integer imsg ! flag for out-of-range warnings (in)
integer itol ! range tolerance option (in)
real ztol ! range tolerance, if itol is set (in)
integer ialg ! algorithm selection code (in)
integer ier ! completion code, 0=OK (out)
call genxpkg(nx,x,xpkg,iper,imsg,itol,ztol,ialg,ier)
c!--------------------
c! set iper=1 if x is a periodic coordinate, e.g. poloidal or
c! toroidal angle coordinate in a toroidal system.
c! if x is periodic, a mod function is applied to
c! normalize input x locations that are out of the
c! given range. x(1:nx) is presumed to cover 1
c! whole period.
c! set iper=0 if x is not a periodic coordinate
c!
c! recommendation: set correctly for the given x coordinate
c!--------------------
c! set imsg=1 for error messages when a given input x location
c! is beyond the range of x(1:nx).
c! set imsg=0 for no error message, but, regardless of imsg, a
c! warning flag will be set by the lookup routine
c!
c! recommendation: set imsg=0 but check warning flags in
c! routines which call for zone lookup / spline
c! evaluation.
c!--------------------
c! set itol=1 to specify relative tolerance for determining
c! even spacing, and for out-of-range errors
c! set itol=0 to use the default relative tolerance, (5.0e-7).
c!
c! set ztol= <the relative tolerance> if itol=1
c!
c! absolute tolerance is then
c! ztola = ztol*max(abs(x(1)),abs(x(nx)));
c! a data point within the range x(1)-ztola to x(nx)+ztola
c! is considered to be "in range".
c!
c! recommendation: the default is OK for most purposes.
c!--------------------
c! to set the search algorithm, in case the grid is unevenly
c! spaced:
c!
c! set ialg= +/- 1: <1/h(j)> "Newton" method
c! +/- 2: binary search
c! +/- 3: piecewise linear indexing function
c!
c! choose the negative value to allow use of the results of
c! the search for the preceding target point, as a start point
c! for the search for the subsequent target point.
c!
c! recommendation: for workstations: ialg = -3 ... is fastest.
c! if the grid is determined to be evenly spaced, this fact
c! will be recorded in xpkg, and the lookup software will
c! take advantage of this, regardless of the setting of "ialg".
c!--------------------
pspline comes with test programs "lookup_test" and "r8lookup_test"
to exercise the even-spaced and uneven-spaced grid zone lookup
algorithms, measuring cpu time consumption.
Tests were run (12 Apr 2000) on a COMPAQ/DEC Alpha 500au using
the DEC FORTRAN optimizing compiler (version 5.1) under
DIGITAL UNIX v4.0.
(a) "vectorization" has a strong performance impact: for example,
(500000 lookups per test)
...even spacing: NO
...algorithm: linear indexing fcn
...DO NOT use previous search result.
%genpkg(ngrid,grid,gpkg, 1,1,0,zdum, 3,ier): ier=
0
$$$$$> cpu time (secs) = 1.312721 vector size = 1
$$$$$> cpu time (secs) = 0.2918243 vector size = 10
$$$$$> cpu time (secs) = 0.1756783 vector size = 100
$$$$$> cpu time (secs) = 0.1668968 vector size = 1000
(this compares 500000 lookups with varying levels of vectorization,
ranging from 500000 single lookup calls to 500 calls with 1000
lookups each).
(b) for vector size = 10, testing 500000 lookups with on an
unevenly spaced grid. Here, numbers in parentheses are from a
Linux P-II 300MHz system running Red Hat 6.1 and using the fujitsu
fortran compiler v1.0.
cpu time, seconds
slowly varying target vector:
-----------------------------
algorithm init search from each iteration
previous iteration independent
Pseudo-Newton .304 (.73) .469 (1.07)
Binary-Search .381 (.85) .508 (.99)
Index-Function .208 (.49) .297 (.67)
rapidly varying target vector:
------------------------------
algorithm init search from each iteration
previous iteration independent
Pseudo-Newton .534 (1.19) .515 (1.20)
Binary-Search .566 (1.14) .537 (1.06)
Index-Function .321 (.71) .286 (.67)
(date of test: 12-Apr 2000)
Pseudo-Newton (genxpkg ialg = +/- 1) -- start in zone j. If that
is not the right zone, use delta(x)/h(j) to estimate which
zone to look in next. delta(x) = distance of target point
from zone j boundary; h(j) = x(j+1)-x(j).
Binary-Search (genxpkg ialg = +/- 2) -- start in zone nx/2. If
that is not the right zone, then, if the zone is too low try
3*nx/4 else try 1*nx/4. Keep subdividing the possible set of
zones by a factor of 2 per step, until the right zone is reached.
We expect log[2](nx) performance scaling of this algorithm.
Index-Function (genxpkg ialg = +/- 3) -- we have pretabulated
on an evenly spaced grid covering x(1) to x(nx) the index
(integer part) and offset with/in zone (fractional part) into
the original, unevenly spaced grid. To do a lookup, we carry out
a piecewise linear interpolation of this indexing function and
extract the integer part of the result. If this yields the
correct zone in the original grid, we are done. If not, we
correct by shifting right or left one zone at a time.
This routine takes a vector of target locations and determines
the corresponding vectors of zone indices and displacements.
It uses "xpkg" (output of subroutine genxpkg-- see above) to
acquire the grid data and lookup algorithm options.
The "imode" switch controls whether to generate location
information required by the explicit spline evaluation
routines (imode=1) or in the form required by compact
spline, Hermite, and piecewise linear routines (imode=2).
Usage:
**call genxpkg first**
real xvec(ivec) ! target vector (in)
real xpkg(nx,4) ! grid xpkg (in)
integer imode ! mode switch (in)
integer iv(ivec) ! list of zone indices (out)
real dxn(ivec) ! displacements w/in zones (out)
real hv(ivec) ! zone widths (out)
real hiv(ivec) ! inverse zone widths (out)
integer iwarn ! warning flag, 0=normal (out)
call xlookup(ivec,xvec,nx,xpkg,imode,iv,dxn,hv,hiv,iwarn)
c! xpkg: must be set up by prior call to genxpkg.
c! xpkg(1:nx,1) = x(1:nx), where x(1:nx) is the original
c! grid.
c!
c! imode switch:
c!
c! if imode=1, then, on output:
c! iv(i) contains the zone in the original grid x(1:nx)
c! in which xvec(i) is contained;
c! i.e. x(i).le.xvec(i).le.x(i+1)
c! iv(i) is always in the range [1,nx-1].
c!
c! dxn(i) = xvec(i) - x(i)
c!
c! (hv(i) and hiv(i) are never set)
c!
c! *** use this for non-compact spline evaluation ONLY ***
c!
c! if imode=2, then, on output:
c! iv(i) contains the zone in the original grid x(1:nx)
c! in which xvec(i) is contained;
c! i.e. x(i).le.xvec(i).le.x(i+1)
c! iv(i) is always in the range [1,nx-1].
c!
c! dxn(i) = (xvec(i) - x(i))/hv(i)
c!
c! hv(i) = (x(i+1)-x(i))
c!
c! hiv(i) = 1/hv(i)
c!
c! *** use this for everything except non-compact splines ***
c!
c! out-of-range warnings (iwarn):
c!
c! iwarn=0 on output: no warnings
c! iwarn=M on output: M points in xvec(1:ivec) were OUT OF RANGE
c!
c! If x(1:nx) is specified periodic (cf genxpkg call), range errors
c! cannot occur; a mod function is used instead to normalize the
c! range of the periodic coordinate; iwarn=0 on output.
c!
c! points xvec(j) out of range low are treated as if xvec(j)=x(1).
c! points xvec(k) out of range high are treated as if xvec(k)=x(nx).
c!
c! If out of range points are within a tolerance of the boundary,
c! the warning flag is not set. The tolerance was set by the
c! genxpkg call.
c!
Routines for setting up spline coefficients:
compact splines memory requirement
---------------------------------------------
1d: mkspline 2*Nx
2d: mkbicub 4*Nx*Ny
3d: mktricub 8*Nx*Ny*Nz
==> these routines create the spline coefficients and any other
information needed by the lookup / evaluation routines --
with various boundary condition control options.
Vectorized Evaluation Routines (recommended):
=============================================
given a location vector, do lookup & evaluation:
vecspline vecbicub vectricub
given a grid (different than the original grid), to
map the function from the old grid to the new grid:
gridspline gridbicub gridtricub
if zone location data is already computed, to do just
the vectorized spline evaluation:
fvspline fvbicub fvtricub
Scalar Routines:
================
Routines for compact spline function evaluation (scalar, grid
lookup included):
evspline evbicub evtricub
subroutine mkspline(x,nx,fspl,ibcxmin,bcxmin,ibcxmax,bcxmax,
> ilinx,ier)
C
C make a 2-coefficient 1d spline
C
C only 2 coefficients, the data and its 2nd derivative, are needed to
C fully specify a spline. See e.g. Numerical Recipies in Fortran-77
C (2nd edition) chapter 3, section on cubic splines.
C
C input:
integer nx ! no. of data points
real x(nx) ! x axis data, strict ascending order
C
C input/output:
real fspl(2,nx) ! f(1,*): data in; f(2,*): coeffs out
C
C f(1,j) = f(x(j)) on input (unchanged on output)
C f(2,j) = f''(x(j)) (of interpolating spline) (on output).
C
C ...boundary conditions...
C
C input:
C
integer ibcxmin ! b.c. flag @ x=xmin=x(1)
real bcxmin ! b.c. data @xmin
C
integer ibcxmax ! b.c. flag @ x=xmax=x(nx)
real bcxcmax ! b.c. data @xmax
C
C ibcxmin=-1 -- periodic boundary condition
C (bcxmin,ibcxmax,bcxmax are ignored)
C
C the output spline s satisfies
C s'(x(1))=s'(x(nx)) ..and.. s''(x(1))=s''(x(nx))
C
C if non-periodic boundary conditions are used, then the xmin and xmax
C boundary conditions can be specified independently:
C
C ibcxmin (ibcxmax) = 0 -- this specifies a "not a knot" boundary
C condition, see "cubsplb.for". This is a common way
C for inferring a "good" spline boundary condition
C automatically from data in the vicinity of the
C boundary. ... bcxmin (bcxmax) are ignored.
C
C ibcxmin (ibcxmax) = 1 -- boundary condition is to have s'(x(1))
C ( s'(x(nx)) ) match the passed value bcxmin (bcxmax).
C
C ibcxmin (ibcxmax) = 2 -- boundary condition is to have s''(x(1))
C ( s''(x(nx)) ) match the passed value bcxmin (bcxmax).
C
C ibcxmin (ibcxmax) = 3 -- boundary condition is to have s'(x(1))=0.0
C ( s'(x(nx))=0.0 )
C
C ibcxmin (ibcxmax) = 4 -- boundary condition is to have s''(x(1))=0.0
C ( s''(x(nx))=0.0 )
C
C ibcxmin (ibcxmax) = 5 -- boundary condition is to have s'(x(1))
C ( s'(x(nx)) ) match the 1st divided difference
C e.g. at x(1): d(1)/h(1), where
C d(j)=f(1,j+1)-f(1,j)
C h(j)=x(j+1)-x(j)
C
C ibcxmin (ibcxmax) = 6 -- BC is to have s''(x(1)) ( s''(x(nx)) )
C match the 2nd divided difference
C e.g. at x(1):
C e(1) = [d(2)/h(2) - d(1)/h(1)]/(0.5*(h(1)+h(2)))
C
C ibcxmin (ibcxmax) = 7 -- BC is to have s'''(x(1)) ( s'''(x(nx)) )
C match the 3rd divided difference
C e.g. at x(1): [e(2)-e(1)]/(0.33333*(h(1)+h(2)+h(3)))
C
C output:
C
integer ilinx ! =1: hint, x axis is ~evenly spaced
C
C let dx[avg] = (x(nx)-x(1))/(nx-1)
C let dx[j] = x(j+1)-x(j), for all j satisfying 1.le.j.lt.nx
C
C if for all such j, abs(dx[j]-dx[avg]).le.(1.0e-3*dx[avg]) then
C ilinx=1 is returned, indicating the data is (at least nearly)
C evenly spaced. Even spacing is useful, for speed of zone lookup,
C when evaluating a spline.
C
C if the even spacing condition is not satisfied, ilinx=2 is returned.
C
integer ier ! exit code, 0=OK
C
C an error code is returned if the x axis is not strict ascending,
C or if nx.lt.4, or if an invalid boundary condition specification was
C input.
subroutine mkbicub(x,nx,y,ny,f,nf2,
> ibcxmin,bcxmin,ibcxmax,bcxmax,
> ibcymin,bcymin,ibcymax,bcymax,
> ilinx,ilinth,ier)
From mkbicub.for:
subroutine mkbicub(x,nx,y,ny,f,nf2,
> ibcxmin,bcxmin,ibcxmax,bcxmax,
> ibcymin,bcymin,ibcymax,bcymax,
> ilinx,iliny,ier)
C
C setup a bicubic spline; store coefficients in compact form
C (as per suggestion of L. Zakharov, PPPL, Feb. 1999)
C 4 coeffs per grid point: f,fxx,fyy,fxxyy
C
C
C input:
integer nx ! length of x vector
integer ny ! length of y vector
real x(nx) ! x vector, strict ascending
real y(ny) ! y vector, strict ascending
C
integer nf2 ! 2nd dimension of f, nf2.ge.nx
C input/output:
C
real f(4,nf2,ny) ! data & spline coefficients
C
C on input: f(1,i,j) = f(x(i),y(j))
C on output: f(1,i,j) unchanged
C f(2,i,j) = d2f/dx2(x(i),y(j))
C f(3,i,j) = d2f/dy2(x(i),y(j))
C f(4,i,j) = d4f/dx2dy2(x(i),y(j))
C
C and the interpolation formula for (x,y) in
C (x(i),x(i+1))^(y(j),y(j+1)) is:
C hx = x(i+1)-x(i) hy = y(j+1)-y(j)
C dxp= (x-x(i))/hx dxm= 1-dxp dxp,dxm in (0,1)
C dyp= (x-x(i))/hx dym= 1-dyp dyp,dym in (0,1)
C dx3p = dxp**3-dxp dx3m = dxm**3-dxm dxp3,dxm3 in (0,1)
C
C finterp = dxm*(dym*f(1,i,j)+dyp*f(1,i,j+1))
C +dxp*(dym*f(1,i+1,j)+dyp*f(1,i+1,j+1))
C +1/6*hx**2*
C dx3m*(dym*f(2,i,j)+dyp*f(2,i,j+1))
C +dx3p*(dym*f(2,i+1,j)+dyp*f(2,i+1,j+1))
C +1/6*hy**2*
C dxm*(dy3m*f(3,i,j)+dy3p*f(3,i,j+1))
C +dxp*(dy3m*f(3,i+1,j)+dy3p*f(3,i+1,j+1))
C +1/36*hx**2*hy**2*
C dx3m*(dym*f(4,i,j)+dyp*f(4,i,j+1))
C +dx3p*(dym*f(4,i+1,j)+dyp*f(4,i+1,j+1))
C
C where the f(2:4,*,*) are cleverly evaluated to assure
C (a) finterp is continuous and twice differentiable across all
C grid cell boundaries, and
C (b) all boundary conditions are satisfied.
C
C the boundary conditions specification options are:
C inputs:
C
integer ibcxmin ! bc flag for x=xmin
real bcxmin(ny) ! bc data vs. y at x=xmin
integer ibcxmax ! bc flag for x=xmax
real bcxmax(ny) ! bc data vs. y at x=xmax
C
integer ibcymin ! bc flag for y=ymin
real bcymin(nx) ! bc data vs. x at y=ymin
integer ibcymax ! bc flag for y=ymax
real bcymax(nx) ! bc data vs. x at y=ymax
C
C with interpretation:
c ibcxmin -- indicator for boundary condition at x(1):
c bcxmin(...) -- boundary condition data
c =-1 -- periodic boundary condition
c =0 -- use "not a knot"
c =1 -- match slope, specified at x(1),th(ith) by bcxmin(ith)
c =2 -- match 2nd deriv., specified at x(1),th(ith) by bcxmin(ith)
c =3 -- boundary condition is slope=0 (df/dx=0) at x(1), all th(j)
c =4 -- boundary condition is d2f/dx2=0 at x(1), all th(j)
c =5 -- match 1st derivative to 1st divided difference
c =6 -- match 2nd derivative to 2nd divided difference
c =7 -- match 3rd derivative to 3rd divided difference
c (for more detailed definition of BCs 5-7, see the
c comments of subroutine mkspline)
c NOTE bcxmin(...) referenced ONLY if ibcxmin=1 or ibcxmin=2
c
c ibcxmax -- indicator for boundary condition at x(nx):
c bcxmax(...) -- boundary condition data
c (interpretation as with ibcxmin, bcxmin)
c NOTE: if ibcxmin=-1, ibcxmax is ignored! ...the BC is periodic.
c
C and analogous interpretation for ibcymin,bcymin,ibcymax,bcymax
C (df/dy or d2f/dy2 boundary conditions at y=ymin and y=ymax).
C
C
C and output arguments
C
integer ilinx ! =1: x grid is "nearly" equally spaced
integer iliny ! =1: y grid is "nearly" equally spaced
C
C ilinx and iliny are set to zero if corresponding grids are not equally
C spaced
C
C and completion code
C
integer ier ! =0 on exit if there is no error.
C
C if there is an error, ier is set and a message is output on unit 6.
C these are considered programming errors in the calling routine.
C
C possible errors:
C x(...) not strict ascending
C y(...) not strict ascending
C nx .lt. 4
C ny .lt. 4
C invalid boundary condition flag
C workspace too small
subroutine mktricub(x,nx,y,ny,z,nz,f,nf2,nf3,
> ibcxmin,bcxmin,ibcxmax,bcxmax,inb1x,
> ibcymin,bcymin,ibcymax,bcymax,inb1y,
> ibczmin,bczmin,ibczmax,bczmax,inb1z,
> ilinx,iliny,ilinz,ier)
c
c setup a tricubic spline; store coefficients in compatct form
c (as per suggestion of L. Zakharov, PPPL, Feb. 1999)
C 8 coeffs per (x,y,z) grid point:
C f,fxx,fyy,fzz,fxxyy,fxxzz,fyyzz,fxxyyzz
C
C input:
integer nx ! length of x vector
integer ny ! length of y vector
integer nz ! length of z vector
real x(nx) ! x vector, strict ascending
real y(ny) ! y vector, strict ascending
real z(nz) ! z vector, strict ascending
c
integer nf2 ! 2nd dim. of f array, nf2.ge.nx
integer nf3 ! 3rd dim. of f array, nf3.ge.ny
c
c input/output:
c
real f(8,nf2,nf3,nz) ! data and spline coefficients
c
C on input: f(1,i,j,k) = f(x(i),y(j),z(k))
C on output: f(1,i,j,k) unchanged
C f(2,i,j,k) = d2f/dx2(x(i),y(j),z(k))
C f(3,i,j,k) = d2f/dy2(x(i),y(j),z(k))
C f(4,i,j,k) = d2f/dz2(x(i),y(j),z(k))
C f(5,i,j,k) = d4f/dx2dy2(x(i),y(j),z(k))
C f(6,i,j,k) = d4f/dx2dz2(x(i),y(j),z(k))
C f(7,i,j,k) = d4f/dy2dz2(x(i),y(j),z(k))
C f(8,i,j,k) = d6f/dx2dy2dz2(x(i),y(j),z(k))
C
C there is a rather Hermite like interpolation formula to go with
C this-- see evtricub.for. Also the bicubic formula is given in
C mkbicubw.for; the tricubic formula is precisely analogous.
C
C boundary condition data
C inputs:
integer inb1x ! 1st dim of xmin & xmax bc arrays
integer inb1y ! 1st dim of ymin & ymax bc arrays
integer inb1z ! 1st dim of zmin & zmax bc arrays
C
integer ibcxmin,ibcxmax ! BC type flag @xmin, xmax
integer ibcymin,ibcymax ! BC type flag @ymin, ymax
integer ibczmin,ibczmax ! BC type flag @zmin, zmax
C
real bcxmin(inb1x,nz)
real bcxmax(inb1x,nz) ! xmin & xmax BC data, ny x nz
real bcymin(inb1y,nz)
real bcymax(inb1y,nz) ! ymin & ymax BC data, nx x nz
real bczmin(inb1z,ny)
real bczmax(inb1z,ny) ! zmin & zmax BC data, nx x ny
c
c where BC data is not required, dummy scalars may be passed.
C the ibc* flags determine whether BC data isneeded.
c
c BC data: bcxmin & bcxmax: BC vs. y,z @xmin,xmax
C bcymin & bcymax: BC vs. x,z @ymin,ymax
C bczmin & bczmax: BC vs. x,y @zmin,zmax
C
c ibcxmin -- indicator for boundary condition at xmin=x(1):
c bcxmin(...) -- boundary condition data
c =-1 -- use periodic boundary condition
c =0 -- use "not a knot"
c =1 -- match slope, specified at x(1),y(iy),z(iz) by bcxmin(iy,iz)
c =2 -- match 2nd derivative, specified at x(1),y(iy),z(iz)
c by bcxmin(iy,iz
c =3 -- boundary condition is slope=0 (df/dx=0) at x(1), all y(j)
c =4 -- boundary condition is d2f/dx2=0 at x(1), all y(j)
c =5 -- match 1st derivative to 1st divided difference
c =6 -- match 2nd derivative to 2nd divided difference
c =7 -- match 3rd derivative to 3rd divided difference
c (for more detailed definition of BCs 5-7, see the
c comments of subroutine mkspline)
c ***NOTE bcxmin(...) referenced ONLY if ibcxmin=1 or ibcxmin=2
c
c ibcxmax -- indicator for boundary condition at x(nx):
c bcxmax(...) -- boundary condition data
c (interpretation as with ibcxmin, bcxmin)
c NOTE: if ibcxmin=-1 then the periodic BC applies on both sides
c and ibcxmax, bcxmax are ignored.
c inb1x -- 1st dimension of bcxmin, bcxmax: if ibcxmin or ibcxmax
c .gt. 0 this must be .ge. ny.
c
c interpretation of ibcymin,bcymin,ibcymax,bcymax,inb1y
c is same as with ibcxmin,... ** inb1y.ge.nx required **
c
c interpretation of ibczmin,bczmin,ibczmax,bczmax,inb1z
c is same as with ibcxmin,... ** inb1z.ge.nx required **
c
c the explicit bdy condition arrays are referenced only if the
c corresponding "ibc" flag values are set to 1 or 2.
c
c output:
integer ilinx ! x vector equal spacing flag
integer iliny ! y vector equal spacing flag
integer ilinz ! z vector equal spacing flag
c
c ilinx -- =1 on output if x(nx) nearly evenly spaced (tol=1e-3)
c iliny -- =1 on output if y(ny) evenly spaced (tol=1e-3)
c ilinz -- =1 on output if z(nz) evenly spaced (tol=1e-3)
c
integer ier ! exit code
c ier -- completion code, 0 for normal
subroutine vecspline(ict,ivec,xvec,ivd,fval,nx,xpkg,fspl,
> iwarn,ier)
c
c vectorized spline evaluation routine -- 1d *compact* spline
c 1. call vectorized zone lookup routine
c 2. call vectorized spline evaluation routine
c
c--------------------------
c input:
integer ict(3) ! selector:
c ict(1)=1 for f (don't evaluate if ict(1)=0)
c ict(2)=1 for df/dx (don't evaluate if ict(2)=0)
c ict(3)=1 for d2f/dx2 (don't evaluate if ict(3)=0)
c
c set ict(1)=3 to get d3f/dx3 (only)
c
integer ivec ! vector dimensioning
c
c ivec-- number of vector pts (spline values to look up)
c
real xvec(ivec) ! x-locations at which to evaluate
c
integer ivd ! 1st dimension of output array
c
c ivd -- 1st dimension of fval, .ge.ivec
c
c output:
real fval(ivd,*) ! output array
c
c fval(1:ivec,1) -- values as per 1st non-zero ict(...) element
c fval(1:ivec,2) -- values as per 2nd non-zero ict(...) element
c --etc--
c
c input:
integer nx ! dimension of spline x grid
real xpkg(nx,4) ! x grid "package" (cf genxpkg)
real fspl(2,nx) ! (compact) spline coefficients
c
c output:
c condition codes, 0 = normal return
integer iwarn ! =1 if an x value was out of range
integer ier ! =1 if argument error detected
subroutine vecbicub(ict,ivec,xvec,yvec,ivd,fval,
> nx,xpkg,ny,ypkg,fspl,inf2,
> iwarn,ier)
c
c vectorized spline evaluation routine -- 2d *compact* spline
c 1. call vectorized zone lookup routine
c 2. call vectorized spline evaluation routine
c
c--------------------------
c input:
integer ict(6) ! selector:
c ict(1)=1 for f (don't evaluate if ict(1)=0)
c ict(2)=1 for df/dx (don't evaluate if ict(2)=0)
c ict(3)=1 for df/dy (don't evaluate if ict(3)=0)
c ict(4)=1 for d2f/dx2 (don't evaluate if ict(4)=0)
c ict(5)=1 for d2f/dy2 (don't evaluate if ict(5)=0)
c ict(6)=1 for d2f/dxdy (don't evaluate if ict(6)=0)
c the number of non zero values ict(1:6)
c determines the number of vector outputs...
c
c new dmc December 2005 -- access to higher derivatives (even if not
c continuous-- but can only go up to 3rd derivatives on any one coordinate.
c if ict(1)=3 -- want 3rd derivatives
c ict(2)=1 for d3f/dx3
c ict(3)=1 for d3f/dx2dy
c ict(4)=1 for d3f/dxdy2
c ict(5)=1 for d3f/dy3
c number of non-zero values ict(2:5) gives no. of outputs
c if ict(1)=4 -- want 4th derivatives
c ict(2)=1 for d4f/dx3dy
c ict(3)=1 for d4f/dx2dy2
c ict(4)=1 for d4f/dxdy3
c number of non-zero values ict(2:4) gives no. of outputs
c if ict(1)=5 -- want 5th derivatives
c ict(2)=1 for d5f/dx3dy2
c ict(3)=1 for d5f/dx2dy3
c number of non-zero values ict(2:3) gives no. of outputs
c if ict(1)=6 -- want 6th derivatives
c d6f/dx3dy3 -- one value is returned.
c
integer ivec ! vector dimensioning
c
c ivec-- number of vector pts (spline values to look up)
c
c list of (x,y) pairs:
c
real xvec(ivec) ! x-locations at which to evaluate
real yvec(ivec) ! y-locations at which to evaluate
c
integer ivd ! 1st dimension of output array
c
c ivd -- 1st dimension of fval, .ge.ivec
c
c output:
real fval(ivd,*) ! output array
c
c fval(1:ivec,1) -- values as per 1st non-zero ict(...) element
c fval(1:ivec,2) -- values as per 2nd non-zero ict(...) element
c --etc--
c
c input:
integer nx,ny ! dimension of spline grids
real xpkg(nx,4) ! x grid "package" (cf genxpkg)
real ypkg(ny,4) ! y grid "package" (cf genxpkg)
integer inf2 ! fspl 3rd array dimension, .ge.nx
real fspl(0:3,inf2,ny) ! (compact) spline coefficients
c
c output:
c condition codes, 0 = normal return
integer iwarn ! =1 if an x value was out of range
integer ier ! =1 if argument error detected
c
subroutine vectricub(ict,ivec,xvec,yvec,zvec,ivd,fval,
> nx,xpkg,ny,ypkg,nz,zpkg,fspl,inf2,inf3,
> iwarn,ier)
c
c vectorized spline evaluation routine -- 3d *compact* spline
c 1. call vectorized zone lookup routine
c 2. call vectorized spline evaluation routine
c
c--------------------------
c input:
integer ict(10) ! selector:
c ict(1)=1 for f (don't evaluate if ict(1)=0)
c ict(2)=1 for df/dx (don't evaluate if ict(2)=0)
c ict(3)=1 for df/dy (don't evaluate if ict(3)=0)
c ict(4)=1 for df/dy (don't evaluate if ict(4)=0)
c ict(5)=1 for d2f/dx2 (don't evaluate if ict(5)=0)
c ict(6)=1 for d2f/dy2 (don't evaluate if ict(6)=0)
c ict(7)=1 for d2f/dz2 (don't evaluate if ict(7)=0)
c ict(8)=1 for d2f/dxdy (don't evaluate if ict(8)=0)
c ict(9)=1 for d2f/dxdz (don't evaluate if ict(9)=0)
c ict(10)=1 for d2f/dydz (don't evaluate if ict(10)=0)
c
c (new dmc Dec 2005 -- higher derivatives available)
c ict(1)=3 --> 3rd derivative, .le.2 diff. in any coordinate
c ict(2:8) select: fxxy fxxz fxyy fxyz fxzz fyyz fyzz
c ->note ict(1)=3, ict(5)=1 gives fxyz = d3f/dxdydz
c ict(1)=-3 --> 3rd derivative, 3 in one coordinate
c ict(2:4) select: fxxx fyyy fzzz
c ict(1)=4 --> 3rd derivative, .le.2 diff. in any coordinate
c ict(2:7) select: fxxyy fxxyz fxxzz fxyyz fxyzz fyyzz
c ict(1)=-4 --> 3rd derivative, 3 in one coordinate
c ict(2:7) select: fxxxy fxxxz fxyyy fxzzz fyyyz fyzzz
c ict(1)=5 --> 3rd derivative, .le.2 diff. in any coordinate
c ict(2:4) select: fxxyyz fxxyzz fxyyzz
c ict(1)=-5 --> 3rd derivative, 3 in one coordinate
c ict(2:10) select: fxxxyy fxxxyz fxxxzz fxxyyy fxxzzz
c fxyyyz fxyzzz fyyyzz fzzzyy
c ict(1)=6 --> 3rd derivative, .le.2 diff. in any coordinate
c fxxyyzz
c ict(1)=-6 --> 3rd derivative, 3 in one coordinate
c ict(2:10) select: fxxxyyy fxxxyyz fxxxyzz fxxxyyz
c fxxyyyz fxxyzzz fxyyyzz fxyyzzz fyyyzzz
c ict(1)=-7 --> 7th derivative
c ict(2:7) select: fxxxyyyz fxxxyyzz fxxxyzzz
c fxxyyyzz fxxyyzzz fxyyyzzz
c ict(1)=-8 --> 8th derivative
c ict(2:4) select: fxxxyyyzz fxxxyyzzz fxxyyyzzz
c ict(1)=-9 --> 9th derivative: fxxxyyyzzz
c
integer ivec ! vector dimensioning
c
c ivec-- number of vector pts (spline values to look up)
c
c list of (x,y,z) triples:
c
real xvec(ivec) ! x-locations at which to evaluate
real yvec(ivec) ! y-locations at which to evaluate
real zvec(ivec) ! z-locations at which to evaluate
c
integer ivd ! 1st dimension of output array
c
c ivd -- 1st dimension of fval, .ge.ivec
c
c output:
real fval(ivd,*) ! output array
c
c fval(1:ivec,1) -- values as per 1st non-zero ict(...) element
c fval(1:ivec,2) -- values as per 2nd non-zero ict(...) element
c --etc--
c
c input:
integer nx,ny,nz ! dimension of spline grids
real xpkg(nx,4) ! x grid "package" (cf genxpkg)
real ypkg(ny,4) ! y grid "package" (cf genxpkg)
real zpkg(nz,4) ! z grid "package" (cf genxpkg)
integer inf2 ! fspl 4th array dimension, .ge.nx
integer inf3 ! fspl 5th array dimension, .ge.ny
real fspl(0:7,inf2,inf3,nz) ! (compact) spline coefficients
c
c output:
c condition codes, 0 = normal return
integer iwarn ! =1 if an x value was out of range
integer ier ! =1 if argument error detected
subroutine gridspline(x_newgrid,nx_new,f_new,nx,xpkg,fspl,
> iwarn,ier)
c
c regrid a spline function f defined vs. x as in xpkg
c to a new grid, given by x_newgrid.
c
c set warning flag if the range x_newgrid exceeds the range of the
c original xpkg.
c
c (xpkg -- see genxpkg subroutine)
c
c input:
c
real x_newgrid(nx_new) ! new grid
c
c output:
c
real f_new(nx_new) ! f evaluated on this grid
c
c input:
c
integer nx ! size of old grid
real xpkg(nx,4) ! old grid "package"
real fspl(nx,2) ! compact spline coefficients of f
c
c output:
c condition codes, =0 for normal exit
c
integer iwarn ! =1 if new grid points out of range
integer ier ! =1 if there is an argument error
c
c--------------------------------------------
c local
c
integer ict(3)
c
data ict/1,0,0/
c
c--------------------------------------------
c
call vecspline(ict,nx_new,x_newgrid,nx_new,f_new,nx,xpkg,fspl,
> iwarn,ier)
c
return
end
subroutine gridbicub(
> x_newgrid,nx_new,y_newgrid,ny_new,f_new,if1,
> nx,xpkg,ny,ypkg,fspl,inf3,iwarn,ier)
c
c regrid a spline function f defined vs. x,y as in xpkg,ypkg
c to a new grid, given by x_newgrid, y_newgrid
c
c set warning flag if the range x_newgrid, y_newgrid
c exceeds the range of the original xpkg,ypkg.
c
c (xpkg,ypkg -- axis data, see genxpkg subroutine)
c
c input:
c
real x_newgrid(nx_new) ! new x grid
real y_newgrid(ny_new) ! new y grid
c
c output:
c
integer if1 ! 1st dimension of f_new
real f_new(if1,ny_new) ! f evaluated on this grid
c
c input:
c
integer nx ! size of old grid
real xpkg(nx,4) ! old grid "package"
integer ny ! size of old grid
real ypkg(ny,4) ! old grid "package"
c
integer inf3 ! array dimension
real fspl(2,2,inf3,ny) ! spline coefficients of f
c
c output:
c condition codes, =0 for normal exit
c
integer iwarn ! =1 if new grid points out of range
integer ier ! =1 if there is an argument error
c
subroutine gridtricub(
> x_newgrid,nx_new,y_newgrid,ny_new,z_newgrid,nz_new,
> f_new,if1,if2,
> nx,xpkg,ny,ypkg,nz,zpkg,fspl,inf4,inf5,iwarn,ier)
c
c regrid a spline function f defined vs. x,y,z as in xpkg, etc.
c to a new grid, given by x_newgrid, y_newgrid, z_newgrid
c
c set warning flag if the range exceeds the range of the
c original x/y/zpkg's.
c
c (xpkg/ypkg/zpkg -- axis data, see genxpkg subroutine)
c
c input:
c
real x_newgrid(nx_new) ! new x grid
real y_newgrid(ny_new) ! new y grid
real z_newgrid(nz_new) ! new z grid
c
c output:
c
integer if1,if2 ! 1st dimensions of f_new
real f_new(if1,if2,nz_new) ! f evaluated on this grid
c
c input:
c
integer nx ! size of old grid
real xpkg(nx,4) ! old grid "package"
integer ny ! size of old grid
real ypkg(ny,4) ! old grid "package"
integer nz ! size of old grid
real zpkg(nz,4) ! old grid "package"
c
integer inf4,inf5 ! array dimensions
real fspl(2,2,2,inf4,inf5,nz) ! spline coefficients of f
c
c output:
c condition codes, =0 for normal exit
c
integer iwarn ! =1 if new grid points out of range
integer ier ! =1 if there is an argument error
c
subroutine fvspline(ict,ivec,ivecd,
> fval,ii,xparam,hx,hxi,fin,nx)
C
C ** use mkspline to set up the fin array **
C
integer ict(3) ! requested output control
integer ivec ! vector length
integer ivecd ! vector dimension (1st dim of fval)
C
integer ii(ivec) ! target cells (i)
real xparam(ivec)
! normalized displacements from (i) corners
C
real hx(ivec) ! grid spacing, and
real hxi(ivec) ! inverse grid spacing 1/(x(i+1)-x(i))
C
real fin(0:1,nx) ! interpolant data (cf "evspline")
C
real fval(ivecd,3) ! output returned
C
C for detailed description of fin, ict and fval see subroutine
C evspline comments. Note ict is not vectorized; the same output
C is expected to be returned for all input vector data points.
C
C note that the index inputs ii and parameter inputs
C xparam,hx,hxi, are vectorized, and the
C output array fval has a vector ** 1st dimension ** whose
C size must be given as a separate argument
C
C to use this routine in scalar mode, pass in ivec=ivecd=1
C
subroutine fvbicub(ict,ivec,ivecd,
> fval,ii,jj,xparam,yparam,hx,hxi,hy,hyi,
> fin,inf2,ny)
C
C ** use mkbicub to set up the fin array **
C
integer ict(6) ! requested output control
integer ivec ! vector length
integer ivecd ! vector dimension (1st dim of fval)
C
integer ii(ivec),jj(ivec) ! target cells (i,j)
real xparam(ivec),yparam(ivec)
! normalized displacements from (i,j) corners
C
real hx(ivec),hy(ivec) ! grid spacing, and
real hxi(ivec),hyi(ivec) ! inverse grid spacing 1/(x(i+1)-x(i))
! & 1/(y(j+1)-y(j))
C
real fin(0:3,inf2,ny) ! interpolant data (cf "evbicub")
C
real fval(ivecd,6) ! output returned
C
C for detailed description of fin, ict and fval see subroutine
C fvbicub comments. Note ict is not vectorized; the same output
C is expected to be returned for all input vector data points.
C
C note that the index inputs ii,jj and parameter inputs
C xparam,yparam,hx,hxi,hy,hyi are vectorized, and the
C output array fval has a vector ** 1st dimension ** whose
C size must be given as a separate argument
C
C to use this routine in scalar mode, pass in ivec=ivecd=1
subroutine fvtricub(ict,ivec,ivecd,
> fval,ii,jj,kk,xparam,yparam,zparam,
> hx,hxi,hy,hyi,hz,hzi,
> fin,inf2,inf3,nz)
C
C use mktricub to set up spline coefficients...
C
integer ict(10) ! requested output control
integer ivec ! vector length
integer ivecd ! vector dimension (1st dim of fval)
C
integer ii(ivec),jj(ivec),kk(ivec) ! target cells (i,j,k)
real xparam(ivec),yparam(ivec),zparam(ivec)
! normalized displacements from (i,j,k) corners
C
real hx(ivec),hy(ivec),hz(ivec) ! grid spacing, and
real hxi(ivec),hyi(ivec),hzi(ivec) ! inverse grid spacing
! 1/(x(i+1)-x(i)) & 1/(y(j+1)-y(j)) & 1/(z(k+1)-z(i))
C
real fin(0:7,inf2,inf3,nz) ! interpolant data (cf "evtricub")
C
real fval(ivecd,10) ! output returned
C
C for detailed description of fin, ict and fval see subroutine
C evtricub comments. Note ict is not vectorized; the same output
C is expected to be returned for all input vector data points.
C
C note that the index inputs ii,jj,kk and parameter inputs
C xparam,yparam,zparam,hx,hxi,hy,hyi,hz,hzi are vectorized, and
C the output array fval has a vector ** 1st dimension ** whose
C size must be given as a separate argument
C
C to use this routine in scalar mode, pass in ivec=ivecd=1
C
subroutine evspline(xget,x,nx,ilinx,f,ict,fval,ier)
C
C evaluate a 1d cubic Spline interpolant -- this is C2
C
C ** use mkspline to set up the f array **
C
C this subroutine calls two subroutines:
C herm1x -- find cell containing (xget)
C fvspline -- evaluate interpolant function and (optionally)
C derivatives
C
C input arguments:
C ================
C
integer nx ! grid size
real xget ! target of this interpolation
real x(nx) ! ordered x grid
integer ilinx ! ilinx=1 => assume x evenly spaced
C
real f(0:1,nx) ! function data
C
C f(0,i) = f @ x(i)
C f(1,i) = d2f/dx2 @ x(i)
C
C (these are spline coefficients selected for continuous 2-
C diffentiability, see mkspline.for)
C
integer ict(3) ! code specifying output desired
C
C ict(1)=1 -- return f (0, don't)
C ict(2)=1 -- return df/dx (0, don't)
C ict(3)=1 -- return d2f/dx2 (0, don't)
C
C set ict(1)=3 to get d3f/dx3 (only)
C
C output arguments:
C =================
C
real fval(3) ! output data
integer ier ! error code =0 ==> no error
C
C fval(1) receives the first output (depends on ict(...) spec)
C fval(2) receives the second output (depends on ict(...) spec)
C fval(3) receives the third output (depends on ict(...) spec)
C
C examples:
C on input ict = [1,1,1]
C on output fval= [f,df/dx,d2f/dx2]
C
C on input ict = [1,0,0]
C on output fval= [f] ... elements 2 -- 3 never referenced.
C
C on input ict = [0,0,1]
C on output fval= [d2f/dx2] ... elements 2 -- 3 never referenced.
C
C ier -- completion code: 0 means OK
subroutine evbicub(xget,yget,x,nx,y,ny,ilinx,iliny,
> f,inf2,ict,fval,ier)
C
C ** use mkbicub to set up the f array **
C
C evaluate a 2d cubic Spline interpolant on a rectilinear
C grid -- this is C2 in both directions.
C
C this subroutine calls two subroutines:
C herm2xy -- find cell containing (xget,yget)
C fvbicub -- evaluate interpolant function and (optionally)
C derivatives
C
C input arguments:
C ================
C
integer nx,ny ! grid sizes
real xget,yget ! target of this interpolation
real x(nx) ! ordered x grid
real y(ny) ! ordered y grid
integer ilinx ! ilinx=1 => assume x evenly spaced
integer iliny ! iliny=1 => assume y evenly spaced
C
real f(0:3,inf2,ny) ! function data
C
C f 2nd dimension inf2 must be .ge. nx
C contents of f:
C
C f(0,i,j) = f @ x(i),y(j)
C f(1,i,j) = d2f/dx2 @ x(i),y(j)
C f(2,i,j) = d2f/dy2 @ x(i),y(j)
C f(3,i,j) = d4f/dx2dy2 @ x(i),y(j)
C
C (these are spline coefficients selected for continuous 2-
C diffentiability, see mkbicub[w].for)
C
integer ict(4) ! code specifying output desired
C
C ict(1)=1 -- return f (0, don't)
C ict(2)=1 -- return df/dx (0, don't)
C ict(3)=1 -- return df/dy (0, don't)
C ict(4)=1 -- return d2f/dx2 (0, don't)
C ict(5)=1 -- return d2f/dy2 (0, don't)
C ict(6)=1 -- return d2f/dxdy (0, don't)
c the number of non zero values ict(1:6)
c determines the number of outputs...
c
c new dmc December 2005 -- access to higher derivatives (even if not
c continuous-- but can only go up to 3rd derivatives on any one coordinate.
c if ict(1)=3 -- want 3rd derivatives
c ict(2)=1 for d3f/dx3
c ict(3)=1 for d3f/dx2dy
c ict(4)=1 for d3f/dxdy2
c ict(5)=1 for d3f/dy3
c number of non-zero values ict(2:5) gives no. of outputs
c if ict(1)=4 -- want 4th derivatives
c ict(2)=1 for d4f/dx3dy
c ict(3)=1 for d4f/dx2dy2
c ict(4)=1 for d4f/dxdy3
c number of non-zero values ict(2:4) gives no. of outputs
c if ict(1)=5 -- want 5th derivatives
c ict(2)=1 for d5f/dx3dy2
c ict(3)=1 for d5f/dx2dy3
c number of non-zero values ict(2:3) gives no. of outputs
c if ict(1)=6 -- want 6th derivatives
c d6f/dx3dy3 -- one value is returned.
C
C output arguments:
C =================
C
real fval(6) ! output data
integer ier ! error code =0 ==> no error
C
C fval(1) receives the first output (depends on ict(...) spec)
C fval(2) receives the second output (depends on ict(...) spec)
C fval(3) receives the third output (depends on ict(...) spec)
C fval(4) receives the fourth output (depends on ict(...) spec)
C fval(5) receives the fourth output (depends on ict(...) spec)
C fval(6) receives the fourth output (depends on ict(...) spec)
C
C examples:
C on input ict = [1,1,1,0,0,1]
C on output fval= [f,df/dx,df/dy,d2f/dxdy], elements 5 & 6 not referenced.
C
C on input ict = [1,0,0,0,0,0]
C on output fval= [f] ... elements 2 -- 6 never referenced.
C
C on input ict = [0,0,0,1,1,0]
C on output fval= [d2f/dx2,d2f/dy2] ... elements 3 -- 6 never referenced.
C
C on input ict = [0,0,1,0,0,0]
C on output fval= [df/dy] ... elements 2 -- 6 never referenced.
C
C ier -- completion code: 0 means OK
subroutine evtricub(xget,yget,zget,x,nx,y,ny,z,nz,
> ilinx,iliny,ilinz,
> f,inf2,inf3,ict,fval,ier)
C
C use mktricub to set up spline coefficients...
C
C evaluate a 3d cubic Spline interpolant on a rectilinear
C grid -- this is C2 in all directions.
C
C this subroutine calls two subroutines:
C herm3xyz -- find cell containing (xget,yget,zget)
C fvtricub -- evaluate the spline function (w/derivatives if
C requested).
C
C input arguments:
C ================
C
real xget,yget,zget ! target of this interpolation
real x(nx) ! ordered x grid
real y(ny) ! ordered y grid
real z(nz) ! ordered z grid
integer ilinx ! ilinx=1 => assume x evenly spaced
integer iliny ! iliny=1 => assume y evenly spaced
integer ilinz ! ilinz=1 => assume z evenly spaced
C
real f(0:7,inf2,inf3,nz) ! function data
C
C f 2nd dimension inf2 must be .ge. nx; 3rd dim inf3 .ge. ny
C contents of f:
C
C f(0,i,j,k) = f @ x(i),y(j),z(k)
C f(1,i,j,k) = d2f/dx2 @ x(i),y(j),z(k)
C f(2,i,j,k) = d2f/dy2 @ x(i),y(j),z(k)
C f(3,i,j,k) = d2f/dz2 @ x(i),y(j),z(k)
C f(4,i,j,k) = d4f/dx2dy2 @ x(i),y(j),z(k)
C f(5,i,j,k) = d4f/dx2dz2 @ x(i),y(j),z(k)
C f(6,i,j,k) = d4f/dy2dz2 @ x(i),y(j),z(k)
C f(7,i,j,k) = d6f/dx2dy2dz2 @ x(i),y(j),z(k)
C
integer ict(10) ! code specifying output desired
C
C ict(1)=1 -- return f (0, don't)
C ict(2)=1 -- return df/dx (0, don't)
C ict(3)=1 -- return df/dy (0, don't)
C ict(4)=1 -- return df/dz (0, don't)
C ict(5)=1 -- return d2f/dx2 (0, don't)
C ict(6)=1 -- return d2f/dy2 (0, don't)
C ict(7)=1 -- return d2f/dz2 (0, don't)
C ict(8)=1 -- return d2f/dxdy (0, don't)
C ict(9)=1 -- return d2f/dxdz (0, don't)
C ict(10)=1-- return d2f/dydz (0, don't)
c
c (new dmc Dec 2005 -- higher derivatives available)
c ict(1)=3 --> 3rd derivative, .le.2 diff. in any coordinate
c ict(2:8) select: fxxy fxxz fxyy fxyz fxzz fyyz fyzz
c ->note ict(1)=3, ict(5)=1 gives fxyz = d3f/dxdydz
c ict(1)=-3 --> 3rd derivative, 3 in one coordinate
c ict(2:4) select: fxxx fyyy fzzz
c ict(1)=4 --> 3rd derivative, .le.2 diff. in any coordinate
c ict(2:7) select: fxxyy fxxyz fxxzz fxyyz fxyzz fyyzz
c ict(1)=-4 --> 3rd derivative, 3 in one coordinate
c ict(2:7) select: fxxxy fxxxz fxyyy fxzzz fyyyz fyzzz
c ict(1)=5 --> 3rd derivative, .le.2 diff. in any coordinate
c ict(2:4) select: fxxyyz fxxyzz fxyyzz
c ict(1)=-5 --> 3rd derivative, 3 in one coordinate
c ict(2:10) select: fxxxyy fxxxyz fxxxzz fxxyyy fxxzzz
c fxyyyz fxyzzz fyyyzz fzzzyy
c ict(1)=6 --> 3rd derivative, .le.2 diff. in any coordinate
c fxxyyzz
c ict(1)=-6 --> 3rd derivative, 3 in one coordinate
c ict(2:10) select: fxxxyyy fxxxyyz fxxxyzz fxxxyyz
c fxxyyyz fxxyzzz fxyyyzz fxyyzzz fyyyzzz
c ict(1)=-7 --> 7th derivative
c ict(2:7) select: fxxxyyyz fxxxyyzz fxxxyzzz
c fxxyyyzz fxxyyzzz fxyyyzzz
c ict(1)=-8 --> 8th derivative
c ict(2:4) select: fxxxyyyzz fxxxyyzzz fxxyyyzzz
c ict(1)=-9 --> 9th derivative: fxxxyyyzzz
C
C output arguments:
C =================
C
real fval(10) ! output data
integer ier ! error code =0 ==> no error
C
C fval(1) receives the first output (depends on ict(...) spec)
C fval(2) receives the second output (depends on ict(...) spec)
C fval(3) receives the third output (depends on ict(...) spec)
C fval(4) receives the 4th output (depends on ict(...) spec)
C fval(5-10) receive 5th thru 10th outputs (if required by ict(...) spec)
C
C examples:
C on input ict = [1,1,1,1,0,0,0,0,0,0,0]
C on output fval= [f,df/dx,df/dy,df/dz]
C
C on input ict = [1,0,0,0,0,0,0,0,0,0,0]
C on output fval= [f] ... elements 2-10 never referenced
C
C on input ict = [0,1,1,0,0,0,0,0,0,0,0]
C on output fval= [df/dx,df/dy] ... elements 3-10 never referenced
C
C on input ict = [0,0,0,0,1,0,0,0,0,0,0]
C on output fval= [d2f/dx2] ... elements 2-10 never referenced.
C
C ier -- completion code: 0 means OK
PSPLINE's "explicit splines" are somewhat faster to evaluate than
PSPLINE's "compact splines", but they require far more memory.
Spline setup routines:
explicit spline routines memory requirement
------------------------------------------------------------
1d: cspline 4*Nx
2d: bcspline 16*Nx*Ny
3d: tcspline 64*Nx*Ny*Nz
==> these routines create the spline coefficients and any other
information needed by the lookup / evaluation routines --
with various boundary condition control options.
Vectorized Evaluation Routines (recommended):
=============================================
given a location vector, do lookup & evaluation:
1d 2d 3d
spvec bcspvec tcspvec
given a grid (different than the original grid), to
map the function from the old grid to the new grid:
spgrid bcspgrid tcspgrid
if zone location data is already computed, to do just
the vectorized spline evaluation:
cspevfn bcspevfn tcspevfn
Scalar Routines:
================
Routines for compact spline function evaluation (scalar, grid
lookup included):
cspeval bcspeval tcspeval
Generally, the spline setup routines require workspaces. As these
routines are written in highly portable fortran-77 code (and thus
lack dynamic memory allocation capability), calling programs need
to declare and provide writable workspace arrays.
The following table summarizes the sizes of workspaces needed.
In this table, "x" denotes the 1st dimension, "y" the 2nd
dimension, and "z" the 3rd dimension of the input data.
=> 1d spline routines:
cspline workspace size: 1*nx (periodic spline option only; for
other options the workspace is unused and a dummy
argument can be provided).
=> 2d spline routines:
bcspline: ws size: 5*max(nx,ny) --or--
5*max(nx,ny)+4*nx*ny [note#1]
=> 3d spline routines:
tcspline: ws size: 10*max(nx,ny,nz) + 20*nx*ny --or--
16*nx*ny*nz [note#2]
[note#1] the larger number is required if the "y" boundary
condition set is "non-linear" (see below).
[note#2] the larger number is required if the "z" boundary
condition set is "non-linear".
A boundary condition set BC is considered "linear" if and only
if, given two splines s1 and s2 satisfying BC, the spline given
by the formula
s = a*s1 + s2 (a is a real number)
also satisfies BC.
Examples of linear boundary condition sets:
* periodic: s'(x(1))=s'(x(nx)) & s''(x(1))=s''(x(nx))
or any combination of two of the following:
* not a knot: s'''(x) continuous at x(2) or x(nx-1).
* zero derivative: s'(x(1))=0 and/or s'(x(nx))=0
* zero 2nd derivative: s''(x(1))=0 and/or s''(x(nx))=0
* a derivative estimated by a divided difference formula
Non-linear boundary condition sets are those for which an
explicit non-zero derivative or 2nd derivative value is
specified. Splines satisfying such boundary conditions
cannot be added without alteration of the boundary condition.
c cspline -- dmc 15 Feb 1999
c
c a standard interface to various 1d spline setup routines
c
subroutine cspline(x,nx,fspl,ibcxmin,bcxmin,ibcxmax,bcxmax,
> wk,iwk,ilinx,ier)
c
real x(nx) ! x axis (in)
real fspl(4,nx) ! spline data (in/out)
integer ibcxmin ! x(1) BC flag (in, see comments)
real bcxmin ! x(1) BC data (in, see comments)
integer ibcxmax ! x(nx) BC flag (in, see comments)
real bcxmax ! x(nx) BC data (in, see comments)
real wk(iwk) ! workspace of size at least nx
c
c workspace only used if "periodic" boundary condition is selected.
c
integer ilinx ! even spacing flag (out)
integer ier ! output status, =0 means OK
c
c this routine computes spline coefficients for a 1d spline --
c evaluation of the spline can be done by cspeval.for subroutines
c or directly by inline code.
c
c the input x axis x(1...nx) must be strictly ascending, i.e.
c x(i+1).gt.x(i) is required for i=1 to nx-1. This is checked and
c ier=1 is set and the routine exits if the test is not satisfied.
c
c on output, ilinx=1 is set if, to a reasonably close tolerance,
c all grid spacings x(i+1)-x(i) are equal. This allows a speedier
c grid lookup algorithm on evaluation of the spline. If on output
c ilinx=2, this means the spline x axis is not evenly spaced.
c
c the input data for the spline are given in f[j] = fspl(1,j). The
c output data are the spline coefficients fspl(2,j),fspl(3,j), and
c fspl(4,j), j=1 to nx. The result is a spline s(x) satisfying the
c boundary conditions and with the properties
c
c s(x(j)) = f(1,j)
c s'(x) is continuous even at the grid points x(j)
c s''(x) is continuous even at the grid points x(j)
c
c the formula for evaluation of s(x) is:
c
c let dx = x-x(i), where x(i).le.x.le.x(i+1). Then,
c s(x)=f(1,i) + dx*(f(2,i) +dx*(f(3,i) + dx*f(4,i)))
c
c ==>boundary conditions. Complete specification of a 1d spline
c requires specification of boundary conditions at x(1) and x(nx).
c
c this routine provides 4 options:
c
c ibcxmin=-1 -- periodic boundary condition. This means the
c boundary conditions s'(x(1))=s'(x(nx)) and s''(x(1))=s''(x(nx))
c are imposed. Note that s(x(1))=s(x(nx)) (i.e. f(1,1)=f(1,nx))
c is not required -- that is determined by the f array input data.
c The periodic boundary condition is to be preferred for periodic
c data. When splining periodic data f(x) with period P, the relation
c x(nx)=x(1)+n*P, n = the number of periods (usually 1), should hold.
c (ibcxmax, bcxmin, bcxmax are ignored).
c
c if a periodic boundary condition is set, this covers both boundaries.
c for the other types of boundary conditions, the type of condition
c chosen for the x(1) boundary need not be the same as the type chosen
c for the x(nx) boundary.
c
c ibcxmin=0 | ibcxmax=0 -- this specifies a "not a knot" boundary
c condition -- see cubsplb.for. This is a common way for inferring
c a "good" spline boundary condition automatically from data in the
c vicinity of the boundary. (bcxmin | bcxmax are ignored).
c
c ibcxmin=1 | ibcxmax=1 -- boundary condition is to have s'(x(1)) |
c s'(x(nx)) match the passed value (bcxmin | bcxmax).
c
c ibcxmin=2 | ibcxmax=2 -- boundary condition is to have s''(x(1)) |
c s''(x(nx)) match the passed value (bcxmin | bcxmax).
c
c ibcxmin=3 | ibcxmax=3 -- boundary condition is to have s'(x(1)) |
c s'(x(nx)) = 0.0 (ZERO)
c
c ibcxmin=4 | ibcxmax=4 -- boundary condition is to have s''(x(1)) |
c s''(x(nx)) = 0.0 (ZERO)
c
c ibcxmin=5 | ibcxmax=5 == BC is to have s' match the 1st divided
c difference at the endpoint.
c
c ibcxmin=6 | ibcxmax=6 == BC is to have s'' match the 2nd divided
c difference at the endpoint.
c
c ibcxmin=7 | ibcxmax=7 == BC is to have s''' match the 3rd divided
c difference at the endpoint.
c
c set up coefficients for bicubic spline with following BC's:
c FULL BC CONTROL at all bdys
c
subroutine bcspline(x,inx,th,inth,fspl,inf3,
> ibcxmin,bcxmin,ibcxmax,bcxmax,
> ibcthmin,bcthmin,ibcthmax,bcthmax,
> wk,nwk,ilinx,ilinth,ier)
c
real x(inx),th(inth),fspl(4,4,inf3,inth),wk(nwk)
real bcxmin(inth),bcxmax(inth)
real bcthmin(inx),bcthmax(inx)
c
c input:
c x(1...inx) -- abscissae, first dimension of data
c th(1...inth) -- abscissae, second (periodic) dimension of data
c fspl(1,1,1..inx,1..inth) -- function values
c inf3 -- fspl dimensioning, inf3.ge.inx required.
c
c boundary conditions input:
c ibcxmin -- indicator for boundary condition at x(1):
c bcxmin(...) -- boundary condition data
c =-1 -- periodic boundary condition
c =0 -- use "not a knot", bcxmin(...) ignored
c =1 -- match slope, specified at x(1),th(ith) by bcxmin(ith)
c =2 -- match 2nd derivative, specified at x(1),th(ith) by
c bcxmin(ith)
c =3 -- boundary condition is slope=0 (df/dx=0) at x(1), all th(j)
c =4 -- boundary condition is d2f/dx2=0 at x(1), all th(j)
c =5 -- match 1st derivative to 1st divided difference
c =6 -- match 2nd derivative to 2nd divided difference
c =7 -- match 3rd derivative to 3rd divided difference
c (for more detailed definition of BCs 5-7, see the
c comments of subroutine mkspline)
c the bcxmin(...) array is referenced only if ibcxmin=1 or
c ibcxmin=2.
c
c ibcxmax -- indicator for boundary condition at x(nx):
c bcxmax(...) -- boundary condition data
c (interpretation as with ibcxmin, bcxmin)
c NOTE: if ibcxmin=-1, ibcxmax is ignored! ...and the BC is periodic.
c
c ibcthmin -- indicator for boundary condition at th(1):
c bcthmin(...) -- boundary condition data
c (interpretation as with ibcxmin, bcxmin)
c ibcthmax -- indicator for boundary condition at th(inth):
c bcthmax(...) -- boundary condition data
c (interpretation as with ibcxmin, bcxmin)
c NOTE: if ibcthmin=-1, ibcthmax is ignored! ...the BC is periodic.
c
c NOTE the bcxmin,bcxmax,bcthmin,bcthmax arrays are only used if the
c corresponding boundary condition flags are set to 1 or 2.
c Carefully note the dimensioning of these arrays!
c
c output:
c fspl(*,*,1..inx,1..inth) -- bicubic spline coeffs (4x4)
c ...fspl(1,1,*,*) is not replaced.
c
c ilinx -- =1 on output if x(inx) pts are nearly evenly spaced
c (tol=1e-3)
c ilinth-- =1 on output if th(inth) are nearly evenly spaced
c (tol=1e-3)
c
c ier -- completion code, 0 for normal
c
c workspace:
c wk -- must be at least 5*max(inx,inth) large *** OR ***
c 4*inx*inth + 5*max(inx,inth) if ibcthmin=1 or 2
c or ibcthmax=1 or 2 and non-zero explicit derivative
c boundary information is provided.
c nwk -- size of workspace provided.
c
c tcspline -- dmc 20 Jan 1999
c
c set up coefficients for bicubic spline with following BC's:
c * LHS and RHS handled as in cubspl.for for 1st coordinate
c * derivatives periodic in second coordinate (use pspline.for)
c
subroutine tcspline(x,inx,th,inth,ph,inph,fspl,inf4,inf5,
> ibcxmin,bcxmin,ibcxmax,bcxmax,inb1x,
> ibcthmin,bcthmin,ibcthmax,bcthmax,inb1th,
> ibcphmin,bcphmin,ibcphmax,bcphmax,inb1ph,
> wk,nwk,ilinx,ilinth,ilinph,ier)
c
real x(inx),th(inth),ph(inph)
real fspl(4,4,4,inf4,inf5,inph),wk(nwk)
real bcxmin(inb1x,inph),bcxmax(inb1x,inph) ! inth x inph used
real bcthmin(inb1th,inph),bcthmax(inb1th,inph) ! inx x inph used
real bcphmin(inb1ph,inth),bcphmax(inb1ph,inth) ! inx x inth used
c
c input:
c x(1...inx) -- abscissae, first dimension of data
c th(1...inth) -- abscissae, second (periodic) dimension of data
c ph(1...inph) -- abscissae, third (periodic) dimension of data
c fspl(1,1,1,1..inx,1..inth,1..inph) -- function values
c inf4 -- fspl dimensioning, inf4.ge.inx required.
c inf5 -- fspl dimensioning, inf5.ge.inth required.
c
c boundary conditions input:
c
c bc data at xmin, xmax vs. theta,phi
c bc data at thmin, thmax vs. x,phi
c bc data at phmin, phmax vs. x,theta
c
c ibcxmin -- indicator for boundary condition at x(1):
c bcxmin(...) -- boundary condition data
c =-1 -- use periodic boundary condition
c =0 -- use "not a knot", bcxmin(...) ignored
c =1 -- match slope, specified at x(1),th(ith),ph(iph)
c by bcxmin(ith,iph)
c =2 -- match 2nd derivative, specified at x(1),th(ith),ph(iph)
c by bcxmin(ith,iph)
c =3 -- boundary condition is slope=0 (df/dx=0) at x(1), all pts
c =4 -- boundary condition is d2f/dx2=0 at x(1), all (th,ph) pts
c =5 -- match 1st derivative to 1st divided difference
c =6 -- match 2nd derivative to 2nd divided difference
c =7 -- match 3rd derivative to 3rd divided difference
c (for more detailed definition of BCs 5-7, see the
c comments of subroutine mkspline)
c bcxmin(...) is referenced ONLY if ibxcmin=1 or ibxcmin=2.
c
c ibcxmax -- indicator for boundary condition at x(nx):
c bcxmax(...) -- boundary condition data
c (interpretation as with ibcxmin, bcxmin)
c
c NOTE: if ibcxmin=-1 then the periodic BC applies on both sides
c and ibcxmax is ignored.
c inb1x -- 1st dimension of bcxmin, bcxmax: if ibcxmin .gt. 0
c or ibcxmax .gt. 0, this must be .ge. inth.
c
c interpretation of ibcthmin,bcthmin,ibcthmax,bcthmax,inb1th
c is same as with ibcxmin,...
c
c interpretation of ibcphmin,bcphmin,ibcphmax,bcphmax,inb1ph
c is same as with ibcxmin,...
c
c output:
c fspl(*,*,*,1..inx,1..inth,1..inph) -- bicubic spline coeffs (4x4)
c ...fspl(1,1,1,*,*,*) is not replaced.
c
c ilinx -- =1 on output if x(inx) pts are nearly evenly spaced
c (tol=1e-3)
c ilinth-- =1 on output if th(inth) evenly spaced (tol=1e-3)
c ilinph-- =1 on output if ph(inph) evenly spaced (tol=1e-3)
c
c ier -- completion code, 0 for normal
c
c workspace:
c wk -- because tcspline uses bcspline as a subroutine, a large
c workspace is needed, depending on boundary conditions:
c
c if ibcphmin = 1 or 2, or, ibcphmax = 1 or 2, and non-zero
c explicit d/dphi or d2/dphi2 boundary condition data is
c provided, the workspace size must be at least
c
c 16*inx*inth*inph.
c
c if ibcphmin and ibcphmax .ne. 1 or 2, and no explicit
c non-zero d/dphi or d2/dphi2 boundary condition data is
c provided, then, the workspace size must be at least:
c
c 20*inx*inth + 10*max(inx,inth,inph).
c
c nwk -- actual size of workspace
c
There is a small collection of routines which have fewer arguments,
but also less control over boundary conditions. These are shown here.
c
c set up coefficients for bicubic spline with following BC's:
c * LHS and RHS handled as in cubspl.for for 1st coordinate
c * derivatives periodic in second coordinate (use pspline.for)
c
subroutine bpsplinb(x,inx,th,inth,fspl,inf3,
> ibcxmin,bcxmin,ibcxmax,bcxmax,
> wk,nwk,ilinx,ilinth,ier)
c
real x(inx),th(inth),fspl(4,4,inf3,inth),wk(nwk)
real bcxmin(inth),bcxmax(inth)
c
c input:
c x(1...inx) -- abscissae, first dimension of data
c th(1...inth) -- abscissae, second (periodic) dimension of data
c fspl(1,1,1..inx,1..inth) -- function values
c inf3 -- fspl dimensioning, inf3.ge.inx required.
c
c boundary conditions input:
c ibcxmin -- indicator for boundary condition at x(1):
c bcxmin(...) -- boundary condition data
c =0 -- use "not a knot", bcxmin(...) ignored
c =1 -- match slope, specified at x(1),th(ith) by bcxmin(ith)
c =2 -- match 2nd derivative, specified at x(1),th(ith) by
c bcxmin(ith)
c =3 -- boundary condition is slope=0 (df/dx=0) at x(1), all th(j)
c =4 -- boundary condition is d2f/dx2=0 at x(1), all th(j)
c =5 -- match 1st derivative to 1st divided difference
c =6 -- match 2nd derivative to 2nd divided difference
c =7 -- match 3rd derivative to 3rd divided difference
c (for more detailed definition of BCs 5-7, see the
c comments of subroutine mkspline)
c ibcxmax -- indicator for boundary condition at x(nx):
c bcxmax(...) -- boundary condition data
c (interpolation as with ibcxmin, bcxmin)
c
c output:
c fspl(*,*,1..inx,1..inth) -- bicubic spline coeffs (4x4)
c ...fspl(1,1,*,*) is not replaced.
c
c ilinx -- =1 on output if x(inx) pts are nearly evenly spaced
c (tol=1e-3)
c ilinth-- =1 on output if th(inth) evenly spaced (tol=1e-3)
c
c ier -- completion code, 0 for normal
c
c workspace:
c wk -- must be at least 5*max(inx,inth) large
c nwk -- size of workspace
c
c set up coefficients for bicubic spline with following BC's:
c * LHS and RHS handled as in spline.for for 1st coordinate
c * derivatives periodic in second coordinate (use pspline.for)
c
c see similar routine, bpsplinb, for more control over x boundary
c conditions...
c
subroutine bpspline(x,inx,th,inth,fspl,inf3,
> wk,nwk,ilinx,ilinth,ier)
c
real x(inx),th(inth),fspl(4,4,inf3,inth),wk(nwk)
integer ilinx,ilinth
c
c input:
c x(1...inx) -- abscissae, first dimension of data
c th(1...inth) -- abscissae, second (periodic) dimension of data
c fspl(1,1,1..inx,1..inth) -- function values
c inf3 -- fspl dimensioning, inf3.ge.inx required.
c
c output:
c fspl(*,*,1..inx,1..inth) -- bicubic spline coeffs (4x4)
c ...fspl(1,1,*,*) is not replaced.
c
c ilinx -- =1 on output if x(inx) pts are nearly evenly spaced
c (tol=1e-3)
c ilinth-- =1 on output if th(inth) evenly spaced (tol=1e-3)
c
c ier -- completion code, 0 for normal
c
c workspace:
c wk -- must be at least 5*max(inx,inth) large
c nwk -- size of workspace
c
c set up coefficients for bicubic spline with following BC's:
c * LHS and RHS handled as in cubspl.for for 1st coordinate
c * derivatives periodic in second coordinate (use pspline.for)
c
subroutine tpsplinb(x,inx,th,inth,ph,inph,fspl,inf4,inf5,
> ibcxmin,bcxmin,ibcxmax,bcxmax,inb1,
> wk,nwk,ilinx,ilinth,ilinph,ier)
c
real x(inx),th(inth),ph(inph)
real fspl(4,4,4,inf4,inf5,inph),wk(nwk)
real bcxmin(inb1,inph),bcxmax(inb1,inph)
c
c input:
c x(1...inx) -- abscissae, first dimension of data
c th(1...inth) -- abscissae, second (periodic) dimension of data
c ph(1...inph) -- abscissae, third (periodic) dimension of data
c fspl(1,1,1,1..inx,1..inth,1..inph) -- function values
c inf4 -- fspl dimensioning, inf4.ge.inx required.
c inf5 -- fspl dimensioning, inf5.ge.inth required.
c
c boundary conditions input:
c ibcxmin -- indicator for boundary condition at x(1):
c bcxmin(...) -- boundary condition data
c =0 -- use "not a knot", bcxmin(...) ignored
c =1 -- match slope, specified at x(1),th(ith),ph(iph)
c by bcxmin(ith,iph)
c =2 -- match 2nd derivative, specified at x(1),th(ith),ph(iph)
c by bcxmin(ith,iph)
c =3 -- boundary condition is slope=0 (df/dx=0) at x(1), all th(j)
c =4 -- boundary condition is d2f/dx2=0 at x(1), all th(j)
c =5 -- match 1st derivative to 1st divided difference
c =6 -- match 2nd derivative to 2nd divided difference
c =7 -- match 3rd derivative to 3rd divided difference
c (for more detailed definition of BCs 5-7, see the
c comments of subroutine mkspline)
c ibcxmax -- indicator for boundary condition at x(nx):
c bcxmax(...) -- boundary condition data
c (interpolation as with ibcxmin, bcxmin)
c
c output:
c fspl(*,*,*,1..inx,1..inth,1..inph) -- bicubic spline coeffs (4x4)
c ...fspl(1,1,1,*,*,*) is not replaced.
c
c ilinx -- =1 on output if x(inx) pts are nearly evenly spaced
c (tol=1e-3)
c ilinth-- =1 on output if th(inth) evenly spaced (tol=1e-3)
c ilinph-- =1 on output if ph(inph) evenly spaced (tol=1e-3)
c
c ier -- completion code, 0 for normal
c
c workspace:
c wk -- must be at least 5*max(inx,inth,inph) large
c nwk -- size of workspace
c set up coefficients for bicubic spline with following BC's:
c * 1st dimension: bounding cubic (divided difference)
c (use spline.for)
c * derivatives periodic in second coordinate (use pspline.for)
c * derivatives periodic in third coordinate (use pspline.for)
c
c for more control over boundary conditions, use tpspline or tcspline.
c
c for evaluation of interpolant, see tcspeval.for
c
subroutine tpspline(x,inx,th,inth,ph,inph,fspl,inf4,inf5,
> wk,nwk,ilinx,ilinth,ilinph,ier)
c
real x(inx),th(inth),ph(inph)
real fspl(4,4,4,inf4,inf5,inph),wk(nwk)
c
c input:
c x(1...inx) -- abscissae, first dimension of data
c th(1...inth) -- abscissae, second (periodic) dimension of data
c ph(1...inph) -- abscissae, third (periodic) dimension of data
c fspl(1,1,1,1..inx,1..inth,1..inph) -- function values
c inf4 -- fspl dimensioning, inf4.ge.inx required.
c inf5 -- fspl dimensioning, inf5.ge.inth required.
c
c output:
c fspl(*,*,*,1..inx,1..inth,1..inph) -- bicubic spline coeffs (4x4)
c ...fspl(1,1,1,*,*,*) is not replaced.
c
c ilinx -- =1 on output if x(inx) pts are nearly evenly spaced
c (tol=1e-3)
c ilinth-- =1 on output if th(inth) evenly spaced (tol=1e-3)
c ilinph-- =1 on output if ph(inph) evenly spaced (tol=1e-3)
c
c ier -- completion code, 0 for normal
c
c workspace:
c wk -- must be at least 5*max(inx,inth,inph) large
c nwk -- size of workspace
c
subroutine spvec(ict,ivec,xvec,ivd,fval,nx,xpkg,fspl,iwarn,ier)
c
c vectorized spline evaluation routine -- 1d spline
c 1. call vectorized zone lookup routine
c 2. call vectorized spline evaluation routine
c
c--------------------------
c input:
integer ict(3) ! selector:
c ict(1)=1 for f (don't evaluate if ict(1)=0)
c ict(2)=1 for df/dx (don't evaluate if ict(2)=0)
c ict(3)=1 for d2f/dx2 (don't evaluate if ict(3)=0)
c
c set ict(1)=3 to get d3f/dx3 (only)
c
integer ivec ! vector dimensioning
c
c ivec-- number of vector pts (spline values to look up)
c
real xvec(ivec) ! x-locations at which to evaluate
c
integer ivd ! 1st dimension of output array
c
c ivd -- 1st dimension of fval, .ge.ivec
c
c output:
real fval(ivd,*) ! output array
c
c fval(1:ivec,1) -- values as per 1st non-zero ict(...) element
c fval(1:ivec,2) -- values as per 2nd non-zero ict(...) element
c --etc--
c
c input:
integer nx ! dimension of spline x grid
real xpkg(nx,4) ! x grid "package" (cf genxpkg)
real fspl(4,nx) ! (non-compact) spline coefficients
c
c output:
c condition codes, 0 = normal return
integer iwarn ! =1 if an x value was out of range
integer ier ! =1 if argument error detected
subroutine bcspvec(ict,ivec,xvec,yvec,ivd,fval,
> nx,xpkg,ny,ypkg,fspl,inf3,
> iwarn,ier)
c
c vectorized spline evaluation routine -- 2d spline
c 1. call vectorized zone lookup routine
c 2. call vectorized spline evaluation routine
c
c--------------------------
c input:
integer ict(6) ! selector:
c ict(1)=1 for f (don't evaluate if ict(1)=0)
c ict(2)=1 for df/dx (don't evaluate if ict(2)=0)
c ict(3)=1 for df/dy (don't evaluate if ict(3)=0)
c ict(4)=1 for d2f/dx2 (don't evaluate if ict(4)=0)
c ict(5)=1 for d2f/dy2 (don't evaluate if ict(5)=0)
c ict(6)=1 for d2f/dxdy (don't evaluate if ict(6)=0)
c the number of non zero values ict(1:6)
c determines the number of vector outputs...
c
c new dmc December 2005 -- access to higher derivatives (even if not
c continuous-- but can only go up to 3rd derivatives on any one coordinate.
c if ict(1)=3 -- want 3rd derivatives
c ict(2)=1 for d3f/dx3
c ict(3)=1 for d3f/dx2dy
c ict(4)=1 for d3f/dxdy2
c ict(5)=1 for d3f/dy3
c number of non-zero values ict(2:5) gives no. of outputs
c if ict(1)=4 -- want 4th derivatives
c ict(2)=1 for d4f/dx3dy
c ict(3)=1 for d4f/dx2dy2
c ict(4)=1 for d4f/dxdy3
c number of non-zero values ict(2:4) gives no. of outputs
c if ict(1)=5 -- want 5th derivatives
c ict(2)=1 for d5f/dx3dy2
c ict(3)=1 for d5f/dx2dy3
c number of non-zero values ict(2:3) gives no. of outputs
c if ict(1)=6 -- want 6th derivatives
c d6f/dx3dy3 -- one value is returned.
c
integer ivec ! vector dimensioning
c
c ivec-- number of vector pts (spline values to look up)
c
c list of (x,y) pairs:
c
real xvec(ivec) ! x-locations at which to evaluate
real yvec(ivec) ! y-locations at which to evaluate
c
integer ivd ! 1st dimension of output array
c
c ivd -- 1st dimension of fval, .ge.ivec
c
c output:
real fval(ivd,*) ! output array
c
c fval(1:ivec,1) -- values as per 1st non-zero ict(...) element
c fval(1:ivec,2) -- values as per 2nd non-zero ict(...) element
c --etc--
c
c input:
integer nx,ny ! dimension of spline grids
real xpkg(nx,4) ! x grid "package" (cf genxpkg)
real ypkg(ny,4) ! y grid "package" (cf genxpkg)
integer inf3 ! fspl 3rd array dimension, .ge.nx
real fspl(4,4,inf3,ny) ! (non-compact) spline coefficients
c
c output:
c condition codes, 0 = normal return
integer iwarn ! =1 if an x value was out of range
integer ier ! =1 if argument error detected
subroutine tcspvec(ict,ivec,xvec,yvec,zvec,ivd,fval,
> nx,xpkg,ny,ypkg,nz,zpkg,fspl,inf4,inf5,
> iwarn,ier)
c
c vectorized spline evaluation routine -- 3d spline
c 1. call vectorized zone lookup routine
c 2. call vectorized spline evaluation routine
c
c--------------------------
c input:
integer ict(10) ! selector:
c ict(1)=1 for f (don't evaluate if ict(1)=0)
c ict(2)=1 for df/dx (don't evaluate if ict(2)=0)
c ict(3)=1 for df/dy (don't evaluate if ict(3)=0)
c ict(4)=1 for df/dy (don't evaluate if ict(4)=0)
c ict(5)=1 for d2f/dx2 (don't evaluate if ict(5)=0)
c ict(6)=1 for d2f/dy2 (don't evaluate if ict(6)=0)
c ict(7)=1 for d2f/dz2 (don't evaluate if ict(7)=0)
c ict(8)=1 for d2f/dxdy (don't evaluate if ict(8)=0)
c ict(9)=1 for d2f/dxdz (don't evaluate if ict(9)=0)
c ict(10)=1 for d2f/dydz (don't evaluate if ict(10)=0)
c
c (new dmc Dec 2005 -- higher derivatives available)
c ict(1)=3 --> 3rd derivative, .le.2 diff. in any coordinate
c ict(2:8) select: fxxy fxxz fxyy fxyz fxzz fyyz fyzz
c ->note ict(1)=3, ict(5)=1 gives fxyz = d3f/dxdydz
c ict(1)=-3 --> 3rd derivative, 3 in one coordinate
c ict(2:4) select: fxxx fyyy fzzz
c ict(1)=4 --> 3rd derivative, .le.2 diff. in any coordinate
c ict(2:7) select: fxxyy fxxyz fxxzz fxyyz fxyzz fyyzz
c ict(1)=-4 --> 3rd derivative, 3 in one coordinate
c ict(2:7) select: fxxxy fxxxz fxyyy fxzzz fyyyz fyzzz
c ict(1)=5 --> 3rd derivative, .le.2 diff. in any coordinate
c ict(2:4) select: fxxyyz fxxyzz fxyyzz
c ict(1)=-5 --> 3rd derivative, 3 in one coordinate
c ict(2:10) select: fxxxyy fxxxyz fxxxzz fxxyyy fxxzzz
c fxyyyz fxyzzz fyyyzz fzzzyy
c ict(1)=6 --> 3rd derivative, .le.2 diff. in any coordinate
c fxxyyzz
c ict(1)=-6 --> 3rd derivative, 3 in one coordinate
c ict(2:10) select: fxxxyyy fxxxyyz fxxxyzz fxxxyyz
c fxxyyyz fxxyzzz fxyyyzz fxyyzzz fyyyzzz
c ict(1)=-7 --> 7th derivative
c ict(2:7) select: fxxxyyyz fxxxyyzz fxxxyzzz
c fxxyyyzz fxxyyzzz fxyyyzzz
c ict(1)=-8 --> 8th derivative
c ict(2:4) select: fxxxyyyzz fxxxyyzzz fxxyyyzzz
c ict(1)=-9 --> 9th derivative: fxxxyyyzzz
c
integer ivec ! vector dimensioning
c
c ivec-- number of vector pts (spline values to look up)
c
c list of (x,y,z) triples:
c
real xvec(ivec) ! x-locations at which to evaluate
real yvec(ivec) ! y-locations at which to evaluate
real zvec(ivec) ! z-locations at which to evaluate
c
integer ivd ! 1st dimension of output array
c
c ivd -- 1st dimension of fval, .ge.ivec
c
c output:
real fval(ivd,*) ! output array
c
c fval(1:ivec,1) -- values as per 1st non-zero ict(...) element
c fval(1:ivec,2) -- values as per 2nd non-zero ict(...) element
c --etc--
c
c input:
integer nx,ny,nz ! dimension of spline grids
real xpkg(nx,4) ! x grid "package" (cf genxpkg)
real ypkg(ny,4) ! y grid "package" (cf genxpkg)
real zpkg(nz,4) ! z grid "package" (cf genxpkg)
integer inf4 ! fspl 4th array dimension, .ge.nx
integer inf5 ! fspl 5th array dimension, .ge.ny
real fspl(4,4,4,inf4,inf5,nz) ! (non-compact) spline coefficients
c
c output:
c condition codes, 0 = normal return
integer iwarn ! =1 if an x value was out of range
integer ier ! =1 if argument error detected
c
subroutine spgrid(x_newgrid,nx_new,f_new,nx,xpkg,fspl,iwarn,ier)
c
c regrid a spline function f defined vs. x as in xpkg
c to a new grid, given by x_newgrid.
c
c set warning flag if the range x_newgrid exceeds the range of the
c original xpkg.
c
c (xpkg -- see genxpkg subroutine)
c
c input:
c
real x_newgrid(nx_new) ! new grid
c
c output:
c
real f_new(nx_new) ! f evaluated on this grid
c
c input:
c
integer nx ! size of old grid
real xpkg(nx,4) ! old grid "package"
real fspl(4,nx) ! spline coefficients of f
c
c output:
c condition codes, =0 for normal exit
c
integer iwarn ! =1 if new grid points out of range
integer ier ! =1 if there is an argument error
c
subroutine bcspgrid(
> x_newgrid,nx_new,y_newgrid,ny_new,f_new,if1,
> nx,xpkg,ny,ypkg,fspl,inf3,iwarn,ier)
c
c regrid a spline function f defined vs. x,y as in xpkg,ypkg
c to a new grid, given by x_newgrid, y_newgrid
c
c set warning flag if the range of x_newgrid or y_newgrid exceeds
c the range of the original xpkg or ypkg.
c
c (xpkg,ypkg arrays -- see genxpkg subroutine)
c
c input:
c
real x_newgrid(nx_new) ! new x grid
real y_newgrid(ny_new) ! new y grid
c
c output:
c
integer if1 ! 1st dimension of f_new
real f_new(if1,ny_new) ! f evaluated on this grid
c
c input:
c
integer nx ! size of old grid
real xpkg(nx,4) ! old grid "package"
integer ny ! size of old grid
real ypkg(ny,4) ! old grid "package"
c
integer inf3 ! array dimension
real fspl(4,4,inf3,ny) ! spline coefficients of f
c
c output:
c condition codes, =0 for normal exit
c
integer iwarn ! =1 if new grid points out of range
integer ier ! =1 if there is an argument error
c
subroutine tcspgrid(
> x_newgrid,nx_new,y_newgrid,ny_new,z_newgrid,nz_new,
> f_new,if1,if2,
> nx,xpkg,ny,ypkg,nz,zpkg,fspl,inf4,inf5,iwarn,ier)
c
c regrid a spline function f defined vs. x,y,z as in xpkg,ypkg,zpkg
c to a new grid, given by x_newgrid, y_newgrid, z_newgrid
c
c set warning flag if the range of x_newgrid or y_newgrid or z_newgrid
c exceeds the range of the original xpkg or ypkg or zpkg.
c
c (xpkg,ypkg,zpkg arrays -- axis data, see genxpkg subroutine)
c
c input:
c
real x_newgrid(nx_new) ! new x grid
real y_newgrid(ny_new) ! new y grid
real z_newgrid(nz_new) ! new z grid
c
c output:
c
integer if1,if2 ! 1st dimensions of f_new
real f_new(if1,if2,nz_new) ! f evaluated on this grid
c
c input:
c
integer nx ! size of old grid
real xpkg(nx,4) ! old grid "package"
integer ny ! size of old grid
real ypkg(ny,4) ! old grid "package"
integer nz ! size of old grid
real zpkg(nz,4) ! old grid "package"
c
integer inf4 ! array dimension
integer inf5 ! array dimension
real fspl(4,4,4,inf4,inf5,nz) ! spline coefficients of f
c
c output:
c condition codes, =0 for normal exit
c
integer iwarn ! =1 if new grid points out of range
integer ier ! =1 if there is an argument error
c
==> note, in real life, if speed is a concern, the caller should
evaluate the spline directly, using the coefficients array.
c cspevfn -- OK now evaluate the cubic spline
c
subroutine cspevfn(ict,ivec,ivd,fval,iv,dxv,f,nx)
c
c input:
integer ict(3) ! selector:
c ict(1)=1 for f (don't evaluate f if ict(1)=0)
c ict(2)=1 for df/dx ""
c ict(3)=1 for d2f/dx2
c
c set ict(1)=3 to get d3f/dx3 (only)
c
integer ivec,ivd ! vector dimensioning
c
c ivec-- number of vector pts (spline values to look up)
c ivd -- 1st dimension of fval, .ge.ivec
c
c output:
real fval(ivd,6) ! output array
c
c v = index to element in vector;
c fval(v,1) = first item requested by ict(...),
c fval(v,2) = 2nd item requested, ...etc...
c
c input:
integer iv(ivec) ! grid cell indices -- vectors
real dxv(ivec) ! displacements w/in cell -- vectors
c
real f(4,nx) ! cubic fcn spline coeffs array
c
c usage example:
c
c 1. for each element xx(v) in a vector of x values:
c find the x zone index and displacement with respect to the
c lower end point of the zone; store thes in vectors iv and dxv.
c
c 2. set ict(1)=0, ict(2)=1, ict(3)=0 -- get only the 1st derivative
c
c 3. ivec is the length of the vector; ivd is the 1st dimension of
c the array fval to receive the output
c
c real fval(ivd,3)
c real xv(ivd)
c real dxv(ivd)
c integer iv(ivd)
c integer ict(3)
c
c real fspline(4,nx) ! spline coeffs
c data ict/0,1,0/ ! this call: want 1st derivative only
c ! this will be output to fval(*,1).
c ...
c do iv=1,ivec
c ... ! compute indices and displacements
c enddo
c call cspevfn(ict,ivec,ivd,fval,iv,dxv,fspline,nx)
c
c bcspevfn -- OK now evaluate the bicubic spline
c
subroutine bcspevfn(ict,ivec,ivd,fval,iv,jv,dxv,dyv,f,inf3,ny)
c
c input:
integer ict(6) ! selector:
c ict(1)=1 for f (don't evaluate f if ict(1)=0)
c ict(2)=1 for df/dx ""
c ict(3)=1 for df/dy ""
c ict(4)=1 for d2f/dx2
c ict(5)=1 for d2f/dy2
c ict(6)=1 for d2f/dxdy
c
c the number of non zero values ict(1:6)
c determines the number of outputs...
c
c new dmc December 2005 -- access to higher derivatives (even if not
c continuous-- but can only go up to 3rd derivatives on any one coordinate.
c if ict(1)=3 -- want 3rd derivatives
c ict(2)=1 for d3f/dx3
c ict(3)=1 for d3f/dx2dy
c ict(4)=1 for d3f/dxdy2
c ict(5)=1 for d3f/dy3
c number of non-zero values ict(2:5) gives no. of outputs
c if ict(1)=4 -- want 4th derivatives
c ict(2)=1 for d4f/dx3dy
c ict(3)=1 for d4f/dx2dy2
c ict(4)=1 for d4f/dxdy3
c number of non-zero values ict(2:4) gives no. of outputs
c if ict(1)=5 -- want 5th derivatives
c ict(2)=1 for d5f/dx3dy2
c ict(3)=1 for d5f/dx2dy3
c number of non-zero values ict(2:3) gives no. of outputs
c if ict(1)=6 -- want 6th derivatives
c d6f/dx3dy3 -- one value is returned.
c
integer ivec,ivd ! vector dimensioning
c
c ivec-- number of vector pts (spline values to look up)
c ivd -- 1st dimension of fval, .ge.ivec
c
c output:
real fval(ivd,6) ! output array
c
c v = index to element in vector;
c fval(v,1) = first item requested by ict(...),
c fval(v,2) = 2nd item requested, ...etc...
c
c input:
integer iv(ivec),jv(ivec) ! grid cell indices -- vectors
real dxv(ivec),dyv(ivec) ! displacements w/in cell -- vectors
c
real f(4,4,inf3,ny) ! bicubic fcn spline coeffs array
integer inf3 ! 3rd dimension of f -- .ge. nx
c
c usage example:
c
c 1. for each element (xx(v),yy(v)) in a vector of (x,y) pairs,
c find the x and y zone indices and displacements with respect
c to the "lower left corner" of the zone; store these in vectors
c iv,jv and dxv,dyv.
c
c 2. set ict(1)=0, ict(2)=1, ict(3)=1, the rest zero -- get only
c the 1st derivatives.
c
c 3. ivec is the length of the vector; ivd is the 1st dimension
c of the array fval to receive the output
c
c real fval(ivd,6)
c real xv(ivd),yv(ivd)
c integer iv(ivd),jv(ivd)
c real dxv(ivd),dyv(ivd)
c integer ict(6)
c
c real fspline(4,4,nx,ny) ! spline coeffs
c data ict/0,1,1,0,0,0/ ! this call: want 1st derivatives
c ! only ... these will be output to
c ! fval(*,1) fval(*,2)
c ...
c do iv=1,ivec
c ... ! find indices and displacements
c enddo
c call bcspevfn(ict,ivec,ivd,fval,iv,jv,dxv,dyv,fspline,nx,ny)
c
c tcspevfn -- OK now evaluate the tricubic spline
c evaluate at a vector set of target locations as specified by
c input vectors (iv,jv,kv), (dxv,dyv,dzv)
c
subroutine tcspevfn(ict,ivec,ivd,fval,iv,jv,kv,dxv,dyv,dzv,
> f,inf4,inf5,nz)
c
c input:
integer ict(10) ! selector:
c ict(1)=1 for f (don't evaluate f if ict(1)=0)
c ict(2)=1 for df/dx ""
c ict(3)=1 for df/dy ""
c ict(4)=1 for df/dz ""
c ict(5)=1 for d2f/dx2
c ict(6)=1 for d2f/dy2
c ict(7)=1 for d2f/dz2
c ict(8)=1 for d2f/dxdy
c ict(9)=1 for d2f/dxdz
c ict(10)=1 for d2f/dydz
c
c (new dmc Dec 2005 -- higher derivatives available)
c ict(1)=3 --> 3rd derivative, .le.2 diff. in any coordinate
c ict(2:8) select: fxxy fxxz fxyy fxyz fxzz fyyz fyzz
c ->note ict(1)=3, ict(5)=1 gives fxyz = d3f/dxdydz
c ict(1)=-3 --> 3rd derivative, 3 in one coordinate
c ict(2:4) select: fxxx fyyy fzzz
c ict(1)=4 --> 3rd derivative, .le.2 diff. in any coordinate
c ict(2:7) select: fxxyy fxxyz fxxzz fxyyz fxyzz fyyzz
c ict(1)=-4 --> 3rd derivative, 3 in one coordinate
c ict(2:7) select: fxxxy fxxxz fxyyy fxzzz fyyyz fyzzz
c ict(1)=5 --> 3rd derivative, .le.2 diff. in any coordinate
c ict(2:4) select: fxxyyz fxxyzz fxyyzz
c ict(1)=-5 --> 3rd derivative, 3 in one coordinate
c ict(2:10) select: fxxxyy fxxxyz fxxxzz fxxyyy fxxzzz
c fxyyyz fxyzzz fyyyzz fzzzyy
c ict(1)=6 --> 3rd derivative, .le.2 diff. in any coordinate
c fxxyyzz
c ict(1)=-6 --> 3rd derivative, 3 in one coordinate
c ict(2:10) select: fxxxyyy fxxxyyz fxxxyzz fxxxyyz
c fxxyyyz fxxyzzz fxyyyzz fxyyzzz fyyyzzz
c ict(1)=-7 --> 7th derivative
c ict(2:7) select: fxxxyyyz fxxxyyzz fxxxyzzz
c fxxyyyzz fxxyyzzz fxyyyzzz
c ict(1)=-8 --> 8th derivative
c ict(2:4) select: fxxxyyyzz fxxxyyzzz fxxyyyzzz
c ict(1)=-9 --> 9th derivative: fxxxyyyzzz
c
integer ivec,ivd ! vector dimensioning
c
c ivec-- number of vector pts (spline values to look up)
c ivd -- 1st dimension of fval, .ge. ivec
c
c output:
real fval(ivd,10) ! output array
c
c for vector elements v, (iv(v),jv(v),kv(v),dxv(v),dyv(v),dzv(v))
c
c fval(v,1) = first item requested by ict(...),
c fval(v,2) = 2nd item requested, ...etc...
c
c input:
integer iv(ivec),jv(ivec),kv(ivec)
! grid cell indices
real dxv(ivec),dyv(ivec),dzv(ivec)
! displacements w/in cell
c
real f(4,4,4,inf4,inf5,nz) ! tricubic fcn spline coeffs array
integer inf4 ! 4th dimension of f, .le.nx
integer inf5 ! 5th dimension of f, .le.ny
c
c usage example:
c
c 1. for each element (xx(v),yy(v),zz(v)) in a vector of (x,y,z)
c triples, find the x,y,z zone indices and displacements with
c to the "lower left corner" of the zone; store these in vectors
c iv,jv,kv and dxv,dyv,dzv
c
c 2. set ict(1)=0, ict(2)=1, ict(3)=1, ict(4)=1 & the rest zero --
c to get only the 1st derivatives.
c
c 3. ivec is the length of the vector; ivd is the 1st dimension
c of the array fval to receive the output
c
c real fval(ivd,10)
c real xv(ivd),yv(ivd),zv(ivd)
c integer iv(ivd),jv(ivd),kv(ivd)
c real dxv(ivd),dyv(ivd),dzv(ivd)
c integer ict(10)
c
c real fspline(4,4,4,nx,ny,nz) ! spline coeffs
c data ict/0,1,1,1,0,0,0,0,0,0/ ! this call: want 1st
c ! derivatives only ... these will
c ! be output to
c ! fval(*,1) fval(*,2) fval(*,3)
c ...
c do iv=1,ivec
c ... ! find indices and displacements
c enddo
c call tcspevfn(ict,ivec,ivd,fval,iv,jv,kv,dxv,dyv,dzv,
c > fspline,nx,ny,nz)
c
c cspeval -- eval cubic spline function and/or derivatives
c
subroutine cspeval(xget,iselect,fval,x,nx,ilinx,f,ier)
c
real xget ! interpolation target
real fval(3) ! output values
real x(nx),f(4,nx) ! spline data
c
integer iselect(3) ! output selector
integer ilinx ! =1 if x(...) is evenly spaced
c
c input:
c (xget) location where interpolated value is desired
c x(1).le.xget.le.x(nx) expected
c
c iselect select desired output
c
c iselect(1)=1 -- want function value (f) itself
c iselect(2)=1 -- want df/dx
c iselect(3)=1 -- want d2f/dx2
c
c example: iselect(1)=iselect(2)=iselect(3)=1
c f, df/dx, and d2f/dx2 are all evaluated
c and returned in order in fval(1), fval(2),
c and fval(3)
c iselect(1)=0, iselect(2)=1, iselect(3)=0
c only the 1st derivative is evaluated
c and returned in fval(1).
c
c set iselect(1)=3 to get d3f/dx3, 1 value only.
c
c see fval (output) description.
c
c x(1...nx) independent coordinate x, strict ascending
c
c ilinx -- =1: flag that x is linearly spaced
c (find zone at once, avoid search for speed)
c
c **CAUTION** actual even spacing of x, is NOT CHECKED HERE!
c
c f the function values (at grid points) and spline
c coefficients
c
c evaluation formula: for point x btw x(i) and x(i+1), dx=x-x(i)
c
c spline value =
c f(1,i) + dx*f(2,i) + dx**2*f(3,i) + dx**3*f(4,i)
c
c output:
c up to 3 elements of fval, ordered as follows:
c fval(1)=function value or lowest order derivative requested
c fval(2)=next order derivative
c etc
c the ordering is a subset of the sequence given under the
c "iselect" description.
c
c ier = 0 -- successful completion; = 1 -- an error occurred.
c bcspeval -- eval bicubic spline function and/or derivatives
c
subroutine bcspeval(xget,yget,iselect,fval,
> x,nx,y,ny,ilinx,iliny,f,inf3,ier)
c
real xget,yget
real fval(6)
real x(nx),y(ny),f(4,4,inf3,ny)
c
integer iselect(6)
integer ilinx,iliny,nx,ny,inf3,ier
c
c input:
c (xget,yget) location where interpolated value is desired
c x(1).le.xget.le.x(nx) expected
c y(1).le.yget.le.y(ny) expected
c
c iselect select desired output
c
c iselect(1)=1 -- want function value (f) itself
c iselect(2)=1 -- want df/dx
c iselect(3)=1 -- want df/dy
c iselect(4)=1 -- want d2f/dx2
c iselect(5)=1 -- want d2f/dy2
c iselect(6)=1 -- want d2f/dxdy
c
c example: iselect(1)=iselect(2)=iselect(3)=1
c f, df/dx, and df/dy all evaluated
c iselect(4)=iselect(5)=iselect(6)=0
c 2nd derivatives not evaluated.
c
c the number of non zero values iselect(1:6)
c determines the number of outputs...
c see fval (output) description.
c
c new dmc December 2005 -- access to higher derivatives (even if not
c continuous-- but can only go up to 3rd derivatives on any one coordinate.
c if iselect(1)=3 -- want 3rd derivatives
c iselect(2)=1 for d3f/dx3
c iselect(3)=1 for d3f/dx2dy
c iselect(4)=1 for d3f/dxdy2
c iselect(5)=1 for d3f/dy3
c number of non-zero values iselect(2:5) gives no. of outputs
c if iselect(1)=4 -- want 4th derivatives
c iselect(2)=1 for d4f/dx3dy
c iselect(3)=1 for d4f/dx2dy2
c iselect(4)=1 for d4f/dxdy3
c number of non-zero values iselect(2:4) gives no. of outputs
c if iselect(1)=5 -- want 5th derivatives
c iselect(2)=1 for d5f/dx3dy2
c iselect(3)=1 for d5f/dx2dy3
c number of non-zero values iselect(2:3) gives no. of outputs
c if iselect(1)=6 -- want 6th derivatives
c d6f/dx3dy3 -- one value is returned.
c
c x(1...nx) independent coordinate x, strict ascending
c y(1...ny) independent coordinate y, strict ascending
c
c ilinx -- =1: flag that x is linearly spaced
c iliny -- =1: flag that y is linearly spaced
c
c **CAUTION** actual even spacing of x, y is NOT CHECKED HERE!
c
c
c f the function values (at grid points)
c and spline coefficients
c
c evaluation formula: for point x btw x(i) and x(i+1), dx=x-x(i)
c and y btw y(j) and y(j+1), dy=y-y(j),
c
c spline value =
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)
c +dy*(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))
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))
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))
c
c where d2=dy**2 and d3=dy**3.
c
c output:
c up to 6 elements of fval, ordered as follows:
c fval(1)=function value or lowest order derivative requested
c fval(2)=next order derivative
c etc
c the ordering is a subset of the sequence given under the
c "iselect" input variable description.
c
c ier = 0 -- successful completion; = 1 -- an error occurred.
c tcspeval -- eval tricubic spline function and/or derivatives
c
subroutine tcspeval(xget,yget,zget,iselect,fval,
> x,nx,y,ny,z,nz,ilinx,iliny,ilinz,f,inf4,inf5,
> ier)
c
real xget,yget,zget
real fval(10)
real x(nx),y(ny),z(nz),f(4,4,4,inf4,inf5,nz)
c
integer iselect(10)
integer ilinx,iliny,ilinz,nx,ny,nz,inf4,inf5,ier
c
c input:
c (xget,yget,zget) location where interpolated value is desired
c x(1).le.xget.le.x(nx) expected
c y(1).le.yget.le.y(ny) expected
c z(1).le.zget.le.z(nz) expected
c
c iselect select desired output
c
c iselect(1)=1 -- want function value (f) itself
c iselect(2)=1 -- want df/dx
c iselect(3)=1 -- want df/dy
c iselect(4)=1 -- want df/dz
c iselect(5)=1 -- want d2f/dx2
c iselect(6)=1 -- want d2f/dy2
c iselect(7)=1 -- want d2f/dz2
c iselect(8)=1 -- want d2f/dxdy
c iselect(9)=1 -- want d2f/dxdz
c iselect(10)=1 -- want d2f/dydz
c
c
c example: iselect(1)=iselect(2)=iselect(3)=iselect(4)=1
c f, df/dx, df/dy, and df/dz all evaluated
c iselect(5)=iselect(6)=iselect(7)=0
c iselect(8)=iselect(9)=iselect(10)=0
c 2nd derivatives not evaluated.
c
c see fval (output) description.
c
c (new dmc Dec 2005 -- higher derivatives available)
c iselect(1)=3 --> 3rd derivative, .le.2 diff. in any coordinate
c iselect(2:8) select: fxxy fxxz fxyy fxyz fxzz fyyz fyzz
c ->note iselect(1)=3, iselect(5)=1 gives fxyz = d3f/dxdydz
c iselect(1)=-3 --> 3rd derivative, 3 in one coordinate
c iselect(2:4) select: fxxx fyyy fzzz
c iselect(1)=4 --> 3rd derivative, .le.2 diff. in any coordinate
c iselect(2:7) select: fxxyy fxxyz fxxzz fxyyz fxyzz fyyzz
c iselect(1)=-4 --> 3rd derivative, 3 in one coordinate
c iselect(2:7) select: fxxxy fxxxz fxyyy fxzzz fyyyz fyzzz
c iselect(1)=5 --> 3rd derivative, .le.2 diff. in any coordinate
c iselect(2:4) select: fxxyyz fxxyzz fxyyzz
c iselect(1)=-5 --> 3rd derivative, 3 in one coordinate
c iselect(2:10) select: fxxxyy fxxxyz fxxxzz fxxyyy fxxzzz
c fxyyyz fxyzzz fyyyzz fzzzyy
c iselect(1)=6 --> 3rd derivative, .le.2 diff. in any coordinate
c fxxyyzz
c iselect(1)=-6 --> 3rd derivative, 3 in one coordinate
c iselect(2:10) select: fxxxyyy fxxxyyz fxxxyzz fxxxyyz
c fxxyyyz fxxyzzz fxyyyzz fxyyzzz fyyyzzz
c iselect(1)=-7 --> 7th derivative
c iselect(2:7) select: fxxxyyyz fxxxyyzz fxxxyzzz
c fxxyyyzz fxxyyzzz fxyyyzzz
c iselect(1)=-8 --> 8th derivative
c iselect(2:4) select: fxxxyyyzz fxxxyyzzz fxxyyyzzz
c iselect(1)=-9 --> 9th derivative: fxxxyyyzzz
c
c x(1...nx) independent coordinate x, strict ascending
c y(1...ny) independent coordinate y, strict ascending
c z(1...nz) independent coordinate y, strict ascending
c
c ilinx -- =1: flag that x is linearly spaced
c iliny -- =1: flag that y is linearly spaced
c ilinz -- =1: flat that z is linearly spaced
c
c **CAUTION** actual even spacing of x, y, z is NOT CHECKED HERE!
c
c
c f the function values (at grid points)
c and spline coefficients
c
c evaluation formula: for point x btw x(i) and x(i+1), dx=x-x(i)
c and y btw y(j) and y(j+1), dy=y-y(j),
c and z btw z(k) and z(k+1), dz=z-z(k)
c
c do m=1,4
c p(m) =
c f(1,1,m,i,j,k)+dx*f(2,1,m,i,j,k)+dx**2*f(3,1,m,i,j,k)+dx**3*f(4,1,m,i,j,k)
c +dy*(
c f(1,2,m,i,j,k)+dx*f(2,2,m,i,j,k)+dx**2*f(3,2,m,i,j,k)+dx**3*f(4,2,m,i,j,k))
c +dy**2*(
c f(1,3,m,i,j,k)+dx*f(2,3,m,i,j,k)+dx**2*f(3,3,m,i,j,k)+dx**3*f(4,3,m,i,j,k))
c +dy**3*(
c f(1,4,m,i,j,k)+dx*f(2,4,m,i,j,k)+dx**2*f(3,4,m,i,j,k)+dx**3*f(4,4,m,i,j,k))
c enddo
c answer = p(1)+dz*p(2)+dz**2*p(3)+dz**3*p(4)
c
c output:
c up to 10 elements of fval, ordered as follows:
c fval(1)=function value or lowest order derivative requested
c fval(2)=next order derivative
c etc
c the ordering is a subset of the sequence given under the
c "iselect" input variable description; the first M elements
c of fval are used, where M = the number of non-zero elements
c of iselect.
c
c ier = 0 -- successful completion; = 1 -- an error occurred.
c
Coefficient setup: spline, pspline
(Forsythe/Malcolm/Moler style spline coefficients computation)
Evaluation: seval, speval
(Forsythe/Malcolm/Moler style spline function evaluation)
C THE CODES (SPLINE & SEVAL) ARE TAKEN FROM:
C FORSYTHE,MALCOLM AND MOLER, "COMPUTER METHODS FOR
C MATHEMATICAL COMPUTATIONS",PRENTICE-HALL, 1977.
C
SUBROUTINE SPLINE (N, X, Y, B, C, D)
INTEGER N
REAL X(N), Y(N), B(N), C(N), D(N)
C
C THE COEFFICIENTS B(I), C(I), AND D(I), I=1,2,...,N ARE COMPUTED
C FOR A CUBIC INTERPOLATING SPLINE
C
C S(X) = Y(I) + B(I)*(X-X(I)) + C(I)*(X-X(I))**2 + D(I)*(X-X(I))**3
C
C FOR X(I) .LE. X .LE. X(I+1)
C
C INPUT..
C
C N = THE NUMBER OF DATA POINTS OR KNOTS (N.GE.2)
C X = THE ABSCISSAS OF THE KNOTS IN STRICTLY INCREASING ORDER
C Y = THE ORDINATES OF THE KNOTS
C
C OUTPUT..
C
C B, C, D = ARRAYS OF SPLINE COEFFICIENTS AS DEFINED ABOVE.
C
C PSPLINE -- modified from SPLINE.FOR -- dmc 3 Oct 1995
C
C THE CODES (SPLINE & SEVAL) ARE TAKEN FROM:
C FORSYTHE,MALCOLM AND MOLER, "COMPUTER METHODS FOR
C MATHEMATICAL COMPUTATIONS",PRENTICE-HALL, 1977.
c
c PSPLINE -- periodic spline -- adaptation of SPLINE.FOR by D. McCune
c October 1995: The function to be splined is taken to have periodic
C derivatives (it need not be periodic itself), so that
c the interpolating function satisfies
c y'(x(1))=y'(x(N))
c y''(x(1))=y''(x(N))
c this results in a fully specified system (no addl BC assumptions
c required); it is not a tridiagonal system but a modified NM1xNM1
c tridiagonal system with non-zero offdiagonal corner (1,NM1),(NM1,1)
c elements. It remains symmetric & diagonally dominant, and is
c solved by Gaussian elimination w/o pivoting. NM1=N-1.
c
SUBROUTINE PSPLINE (N, X, Y, B, C, D, WK)
INTEGER N
REAL X(N), Y(N), B(N), C(N), D(N), WK(N)
C
C THE COEFFICIENTS B(I), C(I), AND D(I), I=1,2,...,N ARE COMPUTED
C FOR A periodic CUBIC INTERPOLATING SPLINE
C
C S(X) = Y(I) + B(I)*(X-X(I)) + C(I)*(X-X(I))**2 + D(I)*(X-X(I))**3
C
C FOR X(I) .LE. X .LE. X(I+1)
C
C WK(...) is used as a workspace during the Guassian elimination.
C it represents the rightmost column of the array
C
C note B(N),C(N),D(N) give coeffs for a periodic continuation into
C the next period, up to X(N)+(X(2)-X(1)) only.
C
C SEVAL can be used to evaluate the spline, but, it is up to the
C SEVAL caller to bring the argument into the normal range
C X(1) to X(N).
C
C INPUT..
C
C N = THE NUMBER OF DATA POINTS OR KNOTS (N.GE.2)
C includes a complete period of the data, the last point being
C a repeat of the first point.
C X = THE ABSCISSAS OF THE KNOTS IN STRICTLY INCREASING ORDER
C X(N)-X(1) is the period of the periodic function represented.
C Y = THE ORDINATES OF THE KNOTS
C
C OUTPUT..
C
C B, C, D = ARRAYS OF SPLINE COEFFICIENTS AS DEFINED ABOVE.
C
(not recommended for MPP use)
REAL FUNCTION SEVAL(N, U, X, Y, B, C, D)
INTEGER N
REAL U, X(N), Y(N), B(N), C(N), D(N)
C
C THIS SUBROUTINE EVALUATES THE CUBIC SPLINE FUNCTION
C
C SEVAL = Y(I) + B(I)*(U-X(I)) + C(I)*(U-X(I))**2 + D(I)*(U-X(I))**3
C
C WHERE X(I) .LT. U .LT. X(I+1), USING HORNER'S RULE
C
C IF U .LT. X(1) THEN I = 1 IS USED.
C IF U .GE. X(N) THEN I = N IS USED.
C
C INPUT..
C
C N = THE NUMBER OF DATA POINTS
C U = THE ABSCISSA AT WHICH THE SPLINE IS TO BE EVALUATED
C X,Y = THE ARRAYS OF DATA ABSCISSAS AND ORDINATES
C B,C,D = ARRAYS OF SPLINE COEFFICIENTS COMPUTED BY SPLINE
C
C IF U IS NOT IN THE SAME INTERVAL AS THE PREVIOUS CALL, THEN A
C BINARY SEARCH IS PERFORMED TO DETERMINE THE PROPER INTERVAL.
C
INTEGER I, J, K
REAL DX
DATA I/1/
IF ( I .GE. N ) I = 1
IF ( U .LT. X(I) ) GO TO 10
IF ( U .LE. X(I+1) ) GO TO 30
C
C BINARY SEARCH
C
10 I = 1
J = N+1
20 K = (I+J)/2
IF ( U .LT. X(K) ) J = K
IF ( U .GE. X(K) ) I = K
IF ( J .GT. I+1 ) GO TO 20
C
C EVALUATE SPLINE
C
30 DX = U - X(I)
SEVAL = Y(I) + DX*(B(I) + DX*(C(I) + DX*D(I)))
C
RETURN
END
(not recommended for MPP use)
REAL FUNCTION SPEVAL(N, U, X, Y, B, C, D)
INTEGER N
REAL U, X(N), Y(N), B(N), C(N), D(N)
C
C THIS SUBROUTINE EVALUATES THE DERIVATIVE OF THE CUBIC SPLINE
C FUNCTION
C
C
C WHERE X(I) .LT. U .LT. X(I+1), USING HORNER'S RULE
C
C IF U .LT. X(1) THEN I = 1 IS USED.
C IF U .GE. X(N) THEN I = N IS USED.
C
C INPUT..
C
C N = THE NUMBER OF DATA POINTS
C U = THE ABSCISSA AT WHICH THE SPLINE IS TO BE EVALUATED
C X,Y = THE ARRAYS OF DATA ABSCISSAS AND ORDINATES
C B,C,D = ARRAYS OF SPLINE COEFFICIENTS COMPUTED BY SPLINE
C
C IF U IS NOT IN THE SAME INTERVAL AS THE PREVIOUS CALL, THEN A
C BINARY SEARCH IS PERFORMED TO DETERMINE THE PROPER INTERVAL.
C
INTEGER I, J, K
REAL DX
DATA I/1/
IF ( I .GE. N ) I = 1
IF ( U .LT. X(I) ) GO TO 10
IF ( U .LE. X(I+1) ) GO TO 30
C
C BINARY SEARCH
C
10 I = 1
J = N+1
20 K = (I+J)/2
IF ( U .LT. X(K) ) J = K
IF ( U .GE. X(K) ) I = K
IF ( J .GT. I+1 ) GO TO 20
C
C EVALUATE SPLINE
C
30 DX = U - X(I)
SPEVAL = B(I) + DX*(2.*C(I) + 3.*DX*D(I))
RETURN
END
Routines for generating Hermite data sets (i.e. derivative
evaluation at grid points):
*using Akima variation minimization method for derivative evaluation:
akherm1, akherm2, akherm3
or Akima with periodic boundary conditions:
akherm1p, akherm2p, akherm3p
(these routines may be treated as examples of ways to create
Hermite data sets)
*using calls to user supplied function for derivative evaluation:
mkherm1, mkherm2, mkherm3
*using simple numerical differentiation on the given grid:
dnherm1, dnherm2, dnherm3
(note: usually better to use the Akima routines)
Vectorized Evaluation Routines (recommended):
=============================================
given a location vector, do lookup & evaluation:
vecherm1 vecherm2 vecherm3
given a grid (different than the original grid), to
map the function from the old grid to the new grid:
gridherm1 gridherm2 gridherm3
if zone location data is already computed, to do just
the vectorized spline evaluation:
herm1fcn herm2fcn herm3fcn
Scalar Evaluation Routines
==========================
herm1ev herm2ev herm3ev
subroutine akherm1p(x,nx,fherm,ilinx,ipx,ier)
C
C Akima subroutine with boundary condition option:
C
C =>ipx=1 for periodic boundary condition
C
C =>ipx=2 for user specified boundary condition:
C in which case fherm(1,1) and fherm(1,nx) must contain
C the derivatives fx(x(1)) and fx(x(nx)), respectively, and these
C will not be modified on output.
C
C =>ipx=0: default boundary conditions
C
C other arguments as with akherm1.
C
C input:
C
integer nx ! array dimensions
real x(nx) ! x coordinate array
real fherm(0:1,nx) ! data/Hermite array
integer ipx ! =1: f' periodic in x
C
C----------------------------
(or call with default boundary condition...)
subroutine akherm1(x,nx,fherm,ilinx,ier)
C
C create a data set for Hermite interpolation, based on Akima's method
C [Hiroshi Akima, Communications of the ACM, Jan 1974, Vol. 17 No. 1]
C
C 1d routine
C
C input:
C
integer nx ! array dimensions
real x(nx) ! x coordinate array
real fherm(0:1,nx) ! data/Hermite array
C
C fherm(0,i) = function value f at x(i) **on input**
C
C fherm(1,i) = derivative df/dx at x(i) **on output**
C
C addl output:
C ilinx=1 if x is "evenly spaced" ier=0 if no errors
C
C ** x must be strict ascending **
C
C ier =0 on out if there is no error.
C
subroutine akherm2p(x,nx,y,ny,fherm,nf2,ilinx,iliny,ipx,ipy,ier)
C
C create a data set for Hermite interpolation, based on Akima's method
C [Hiroshi Akima, Communications of the ACM, Jan 1974, Vol. 17 No. 1]
C
C with independently settable boundary condition options:
C ipx or ipy =0 -- default, boundary conditions from divided diffs
C ipx or ipy =1 -- periodic boundary condition
C ipx or ipy =2 -- user supplied df/dx or df/dy
C input:
C
integer nx,ny,nf1 ! array dimensions
real x(nx) ! x coordinate array
real y(ny) ! y coordinate array
real fherm(0:3,nf2,ny) ! data/Hermite array
C
integer ipx ! =1 if df/dx periodic in x
integer ipy ! =1 if df/dy periodic in y
C
C fherm(0,1:nx,1:ny) supplied by user; this routine computes the
C rest of the elements, but note:
C
C if ipx=2: fherm(1,1,1:ny) & fherm(1,nx,1:ny) = INPUT df/dx BCs
C if ipy=2: fherm(2,1:nx,1) & fherm(2,1:nx,ny) = INPUT df/dy BCs
C
C on output, at all grid points (i,j) covering (1:nx,1:ny):
C fherm1(1,i,j) -- df/dx at the grid point
C fherm1(2,i,j) -- df/dy at the grid point
C fherm1(3,i,j) -- d2f/dxdy at the grid point
C
C---------------------
(or call with default boundary condition...)
subroutine akherm2(x,nx,y,ny,fherm,nf2,ilinx,iliny,ier)
C
C create a data set for Hermite interpolation, based on Akima's method
C [Hiroshi Akima, Communications of the ACM, Jan 1974, Vol. 17 No. 1]
C
C input:
C
integer nx,ny,nf1 ! array dimensions
real x(nx) ! x coordinate array
real y(ny) ! y coordinate array
real fherm(0:3,nf2,ny) ! data/Hermite array
C
C fherm(0,i,j) = function value f at x(i),y(j) **on input**
C
C fherm(1,i,j) = derivative df/dx at x(i),y(j) **on output**
C fherm(2,i,j) = derivative df/dy at x(i),y(j) **on output**
C fherm(3,i,j) = derivative d2f/dxdy at x(i),y(j) **on output**
C
C addl output:
C ilinx=1 if x axis is evenly spaced
C iliny=1 if y axis is evenly spaced
C ier=0 if no error:
C x, y must both be strict ascending
C nf2.ge.nx is required.
C
C----------------------------
subroutine akherm3p(x,nx,y,ny,z,nz,fherm,nf2,nf3,
> ilinx,iliny,ilinz,ipx,ipy,ipz,ier)
C
C create a data set for Hermite interpolation, based on Akima's method
C [Hiroshi Akima, Communications of the ACM, Jan 1974, Vol. 17 No. 1]
C
C 3d routine --
C with independently settable boundary condition options:
C ipx|ipy|ipz =0 -- default, boundary conditions from divided diffs
C ipx|ipy|ipz =1 -- periodic boundary condition
C ipx|ipy|ipz =2 -- user supplied df/dx|df/dy|df/dz
C
C input:
C
C nf2.get.nx and nf3.ge.ny expected!
C
integer nx,ny,nz,nf2,nf3 ! array dimensions
real x(nx) ! x coordinate array
real y(ny) ! y coordinate array
real z(nz) ! z coordinate array
real fherm(0:7,nf2,nf3,nz) ! data/Hermite array
C
integer ipx ! =0: dflt; =1: periodic df/dx
integer ipy ! =0: dflt; =1: periodic df/dy
integer ipz ! =0: dflt; =1: periodic df/dz
C
C if ipx=2: fherm(1,1,1:ny,1:nz) & fherm(1,nx,1:ny,1:nz) = INPUT df/dx BCs
C if ipy=2: fherm(2,1:nx,1,1:nz) & fherm(2,1:nx,ny,1:nz) = INPUT df/dy BCs
C if ipz=2: fherm(3,1:nx,1:ny,1) & fherm(3,1:nx,1:ny,nz) = INPUT df/dz BCs
C
C for the rest of the grid points...
C
C fherm(0,i,j,k) = function value f at x(i),y(j),z(k) **on input**
C
C fherm(1,i,j,k) = derivative df/dx at x(i),y(j),z(k) **on output**
C fherm(2,i,j,k) = derivative df/dy at x(i),y(j),z(k) **on output**
C fherm(3,i,j,k) = derivative df/dz at x(i),y(j),z(k) **on output**
C fherm(4,i,j,k) = derivative d2f/dxdy at x(i),y(j),z(k) **on output**
C fherm(5,i,j,k) = derivative d2f/dxdz at x(i),y(j),z(k) **on output**
C fherm(6,i,j,k) = derivative d2f/dydz at x(i),y(j),z(k) **on output**
C fherm(7,i,j,k) = derivative d3f/dxdydz at x(i),y(j),z(k) **on output**
C
C additional outputs:
C
C ilinx = 1 if x is evenly spaced (approximately)
C iliny = 1 if y is evenly spaced (approximately)
C ilinz = 1 if z is evenly spaced (approximately)
C
C---------------------
(or call with default boundary condition...)
subroutine akherm3(x,nx,y,ny,z,nz,fherm,nf2,nf3,
> ilinx,iliny,ilinz,ier)
C
C create a data set for Hermite interpolation, based on Akima's method
C [Hiroshi Akima, Communications of the ACM, Jan 1974, Vol. 17 No. 1]
C
C 3d routine
C
C input:
C
C nf2.get.nx and nf3.ge.ny expected!
C
integer nx,ny,nz,nf2,nf3 ! array dimensions
real x(nx) ! x coordinate array
real y(ny) ! y coordinate array
real z(nz) ! z coordinate array
real fherm(0:7,nf2,nf3,nz) ! data/Hermite array
C
C fherm(0,i,j,k) = function value f at x(i),y(j),z(k) **on input**
C
C fherm(1,i,j,k) = derivative df/dx at x(i),y(j),z(k) **on output**
C fherm(2,i,j,k) = derivative df/dy at x(i),y(j),z(k) **on output**
C fherm(3,i,j,k) = derivative df/dz at x(i),y(j),z(k) **on output**
C fherm(4,i,j,k) = derivative d2f/dxdy at x(i),y(j),z(k) **on output**
C fherm(5,i,j,k) = derivative d2f/dxdz at x(i),y(j),z(k) **on output**
C fherm(6,i,j,k) = derivative d2f/dydz at x(i),y(j),z(k) **on output**
C fherm(7,i,j,k) = derivative d3f/dxdydz at x(i),y(j),z(k) **on output**
C
C additional outputs:
C
C ilinx = 1 if x is evenly spaced (approximately)
C iliny = 1 if y is evenly spaced (approximately)
C ilinz = 1 if z is evenly spaced (approximately)
C
C ier = 0 if there are no errors
C
C note possible errors: x y or z NOT strict ascending
C nf2.lt.nx .or. nf3.lt.ny
C
subroutine mkherm1(fun,x,nx,fherm)
C
C create a data set for Hermite interpolation, from evaluation of
C function and derivatives. function of 1 indep. coordinate.
C
C input:
C
external fun ! passed real function(x,dfdx)
real x(nx) ! x coordinate array
C
C the passed function fun must have the interface:
C
C real function <name>(x,dfdx)
C where x is input, the function returns the function value and the
C argument dfdx returns as output the function derivative at x.
C
C output:
C
real fherm(0:1,nx) ! function data & derivative
C
C fherm(0,j) = function value f at x(j)
C fherm(1,j) = derivative df/dx at x(j)
C
C----------------------------
C
do ix=1,nx
fherm(0,ix)=fun(x(ix),dfdx)
fherm(1,ix)=dfdx
enddo
C
return
end
subroutine mkherm2(fun,x,nx,y,ny,fherm)
C
C create a data set for Hermite interpolation, from evaluation of
C function and derivatives. function of 2 indep. coordinates.
C
C input:
C
external fun ! passed real function(x,y)
real x(nx) ! x coordinate array
real y(ny) ! y coordinate array
C
C the passed function fun must have the interface:
C
C real function <name>(x,y,dfdx,dfdy,d2fdxdy)
C where x,y are input, the function returns the function value,
C and the arguments dfdx and dfdy return as output the function
C derivative at the point (x,y).
C
C output:
C
real fherm(0:3,nx,ny) ! function data & derivatives
C
C fherm(0,i,j) = function value f at x(i),y(j)
C fherm(1,i,j) = derivative df/dx at x(i),y(j)
C fherm(2,i,j) = derivative df/dy at x(i),y(j)
C fherm(3,i,j) = derivative d2f/dxdy at x(i),y(j)
C
C----------------------------
C
do ix=1,nx
do iy=1,ny
fherm(0,ix,iy)=fun(x(ix),y(iy),dfdx,dfdy,d2fdxdy)
fherm(1,ix,iy)=dfdx
fherm(2,ix,iy)=dfdy
fherm(3,ix,iy)=d2fdxdy
enddo
enddo
C
return
end
subroutine mkherm3(fun,x,nx,y,ny,z,nz,fherm)
C
C create a data set for Hermite interpolation, from evaluation of
C function and derivatives. Function of 3 indep. coordinates.
C
C input:
C
external fun ! passed real function(x,y,z)
real x(nx) ! x coordinate array (1st dim)
real y(ny) ! y coordinate array (2nd dim)
real z(nz) ! z coordinate array (3rd dim)
C
C the passed function fun must have the interface:
C
C real function <name>(x,y,z,dfdx,dfdy,dfdz,
C d2fdxdy,d2fdxdz,d2fdydz,d3fdxdydz)
C where x,y,z are input, the function returns the function value,
C and the arguments dfdx, dfdy and dfdz return as output the function
C derivative at the point (x,y,z).
C
C output:
C
real fherm(0:7,nx,ny,nz) ! function data & derivatives
C
C fherm(0,i,j,k) = function value f at x(i),y(j),z(k)
C fherm(1,i,j,k) = derivative df/dx at x(i),y(j),z(k)
C fherm(2,i,j,k) = derivative df/dy at x(i),y(j),z(k)
C fherm(3,i,j,k) = derivative df/dz at x(i),y(j),z(k)
C fherm(4,i,j,k) = derivative d2f/dxdy at x(i),y(j),z(k)
C fherm(5,i,j,k) = derivative d2f/dxdz at x(i),y(j),z(k)
C fherm(6,i,j,k) = derivative d2f/dydz at x(i),y(j),z(k)
C fherm(7,i,j,k) = derivative d3f/dxdydz at x(i),y(j),z(k)
C
C----------------------------
C
do iz=1,nz
do iy=1,ny
do ix=1,nx
fherm(0,ix,iy,iz)=
> fun(x(ix),y(iy),z(iz),dfdx,dfdy,dfdz,
> d2fdxdy,d2fdxdz,d2fdydz,d3fdxyz)
fherm(1,ix,iy,iz)=dfdx
fherm(2,ix,iy,iz)=dfdy
fherm(3,ix,iy,iz)=dfdz
fherm(4,ix,iy,iz)=d2fdxdy
fherm(5,ix,iy,iz)=d2fdxdz
fherm(6,ix,iy,iz)=d2fdydz
fherm(7,ix,iy,iz)=d3fdxyz
enddo
enddo
enddo
C
return
end
subroutine dnherm1(x,nx,fherm,ilinx,ier)
C
C create a data set for Hermite interpolation, based on simple
C numerical differentiation using the given grid.
C
C 1d routine
C
C input:
C
integer nx ! array dimensions
real x(nx) ! x coordinate array
real fherm(0:1,nx) ! data/Hermite array
C
C fherm(0,i) = function value f at x(i) **on input**
C
C fherm(1,i) = derivative df/dx at x(i) **on output**
C
C addl output:
C ilinx=1 if x is "evenly spaced" ier=0 if no errors
C
C ** x must be strict ascending **
C
C----------------------------
C
ztol=1.0e-3
call splinck(x,nx,ilinx,ztol,ier)
if(ier.ne.0) then
write(6,*) '?dnherm1: x axis not strict ascending.'
return
endif
C
do ix=1,nx
c
c x div. diffs in vicinity
c
ixp=min(nx,ix+1)
ixm=max(1,ix-1)
zd=(fherm(0,ixp)-fherm(0,ixm))/(x(ixp)-x(ixm))
c
fherm(1,ix)=zd
enddo
C
return
end
subroutine dnherm2(x,nx,y,ny,fherm,nf2,ilinx,iliny,ier)
C
C create a data set for Hermite interpolation, based on simple
C numerical differentiation using the given grid. The grid should
C be "approximately" evenly spaced.
C
C 2d routine
C
C input:
C
C nf2.gt.nx expected!
C
integer nx,ny,nf2 ! array dimensions
real x(nx) ! x coordinate array
real y(ny) ! y coordinate array
real fherm(0:3,nf2,ny) ! data/Hermite array
C
C fherm(0,i,j) = function value f at x(i),y(j) **on input**
C
C fherm(1,i,j) = derivative df/dx at x(i),y(j) **on output**
C fherm(2,i,j) = derivative df/dy at x(i),y(j) **on output**
C fherm(3,i,j) = derivative d2f/dxdy at x(i),y(j) **on output**
C
C addl output:
C ilinx=1 if x axis is evenly spaced
C iliny=1 if y axis is evenly spaced
C ier=0 if no error:
C x, y must both be strict ascending
C nf2.ge.nx is required.
C
C----------------------------
c
ier=0
c
call splinck(x,nx,ilinx,1.0e-3,ierx)
if(ierx.ne.0) ier=2
c
if(ier.eq.2) then
write(6,'('' ?dnherm2: x axis not strict ascending'')')
endif
c
call splinck(y,ny,iliny,1.0e-3,iery)
if(iery.ne.0) ier=3
c
if(ier.eq.3) then
write(6,'('' ?dnherm2: y axis not strict ascending'')')
endif
c
if(nf2.lt.nx) then
ier=4
write(6,*) '?dnherm2: fherm array dimension too small.'
endif
C
if(ier.ne.0) return
C
do iy=1,ny
c
iyp=min(ny,iy+1)
iym=max(1,iy-1)
c
do ix=1,nx
c
ixp=min(nx,ix+1)
ixm=max(1,ix-1)
c
c x diffs in vicinity
c
zd=(fherm(0,ixp,iy)-fherm(0,ixm,iy))/(x(ixp)-x(ixm))
c
fherm(1,ix,iy)=zd
c
c y diffs in vicinity
c
zd=(fherm(0,ix,iyp)-fherm(0,ix,iym))/(y(iyp)-y(iym))
c
fherm(2,ix,iy)=zd
c
c xy cross derivs (except at corners, iedge=2)
c
fherm(3,ix,iy)=(fherm(0,ixp,iyp)-fherm(0,ixm,iyp)
> -fherm(0,ixp,iym)+fherm(0,ixm,iym))/
> ((x(ixp)-x(ixm))*(y(iyp)-y(iym)))
c
enddo
enddo
C
return
end
subroutine dnherm3(x,nx,y,ny,z,nz,fherm,nf2,nf3,
> ilinx,iliny,ilinz,ier)
C
C create a data set for Hermite interpolation, based on simple
C numerical differentiation using the given grid. The grid should
C be "approximately" evenly spaced.
C
C 3d routine
C
C input:
C
C nf2.get.nx and nf3.ge.ny expected!
C
integer nx,ny,nz,nf2,nf3 ! array dimensions
real x(nx) ! x coordinate array
real y(ny) ! y coordinate array
real z(nz) ! z coordinate array
real fherm(0:7,nf2,nf3,nz) ! data/Hermite array
C
C fherm(0,i,j,k) = function value f at x(i),y(j),z(k) **on input**
C
C fherm(1,i,j,k) = derivative df/dx at x(i),y(j),z(k) **on output**
C fherm(2,i,j,k) = derivative df/dy at x(i),y(j),z(k) **on output**
C fherm(3,i,j,k) = derivative df/dz at x(i),y(j),z(k) **on output**
C fherm(4,i,j,k) = derivative d2f/dxdy at x(i),y(j),z(k) **on output**
C fherm(5,i,j,k) = derivative d2f/dxdz at x(i),y(j),z(k) **on output**
C fherm(6,i,j,k) = derivative d2f/dydz at x(i),y(j),z(k) **on output**
C fherm(7,i,j,k) = derivative d3f/dxdydz at x(i),y(j),z(k) **on output**
C
C additional outputs:
C
C ilinx = 1 if x is evenly spaced (approximately)
C iliny = 1 if y is evenly spaced (approximately)
C ilinz = 1 if z is evenly spaced (approximately)
C
C ier = 0 if there are no errors
C
C note possible errors: x y or z NOT strict ascending
C nf2.lt.nx .or. nf3.lt.ny
C
C----------------------------
C
C
C error check
C
ier=0
c
call splinck(x,nx,ilinx,1.0e-3,ierx)
if(ierx.ne.0) ier=2
c
if(ier.eq.2) then
write(6,'('' ?dnherm3: x axis not strict ascending'')')
endif
c
call splinck(y,ny,iliny,1.0e-3,iery)
if(iery.ne.0) ier=3
c
if(ier.eq.3) then
write(6,'('' ?dnherm3: y axis not strict ascending'')')
endif
c
call splinck(z,nz,ilinz,1.0e-3,ierz)
if(ierz.ne.0) ier=4
c
if(ier.eq.4) then
write(6,'('' ?dnherm3: z axis not strict ascending'')')
endif
c
if(nf2.lt.nx) then
ier=5
write(6,*) '?dnherm3: fherm (x) array dimension too small.'
endif
c
if(nf3.lt.ny) then
ier=6
write(6,*) '?dnherm3: fherm (y) array dimension too small.'
endif
C
if(ier.ne.0) return
C
C error check OK
C
do iz=1,nz
c
izp=min(nz,iz+1)
izm=max(1,iz-1)
c
do iy=1,ny
c
iyp=min(ny,iy+1)
iym=max(1,iy-1)
c
do ix=1,nx
c
ixp=min(nx,ix+1)
ixm=max(1,ix-1)
c
c x diffs in vicinity
c
zd=(fherm(0,ixp,iy,iz)-fherm(0,ixm,iy,iz))/
> (x(ixp)-x(ixm))
c
fherm(1,ix,iy,iz)=zd
c
c y diffs in vicinity
c
zd=(fherm(0,ix,iyp,iz)-fherm(0,ix,iym,iz))/
> (y(iyp)-y(iym))
c
fherm(2,ix,iy,iz)=zd
c
c z diffs in vicinity
c
zd=(fherm(0,ix,iy,izp)-fherm(0,ix,iy,izm))/
> (z(izp)-z(izm))
c
fherm(3,ix,iy,iz)=zd
c
c xy cross derivs
c
fherm(4,ix,iy,iz)=
> (fherm(0,ixp,iyp,iz)-fherm(0,ixm,iyp,iz)
> -fherm(0,ixp,iym,iz)+fherm(0,ixm,iym,iz))/
> ((x(ixp)-x(ixm))*(y(iyp)-y(iym)))
c
c xz cross derivs
c
fherm(5,ix,iy,iz)=
> (fherm(0,ixp,iy,izp)-fherm(0,ixm,iy,izp)
> -fherm(0,ixp,iy,izm)+fherm(0,ixm,iy,izm))/
> ((x(ixp)-x(ixm))*(z(izp)-z(izm)))
c
c yz cross derivs
c
fherm(6,ix,iy,iz)=
> (fherm(0,ix,iyp,izp)-fherm(0,ix,iym,izp)
> -fherm(0,ix,iyp,izm)+fherm(0,ix,iym,izm))/
> ((y(iyp)-y(iym))*(z(izp)-z(izm)))
c
c xyz cross deriv
c
fherm(7,ix,iy,iz)=
> ((fherm(0,ixp,iyp,izp)-fherm(0,ixp,iym,izp)
> -fherm(0,ixp,iyp,izm)+fherm(0,ixp,iym,izm))-
> (fherm(0,ixm,iyp,izp)-fherm(0,ixm,iym,izp)
> -fherm(0,ixm,iyp,izm)+fherm(0,ixm,iym,izm)))/
> ((x(ixp)-x(ixm))*(y(iyp)-y(iym))*(z(izp)-z(izm)))
c
enddo
enddo
enddo
C
return
end
subroutine vecherm1(ict,ivec,xvec,ivd,fval,nx,xpkg,fspl,
> iwarn,ier)
c
c vectorized hermite evaluation routine -- 1d
c 1. call vectorized zone lookup routine
c 2. call vectorized hermite evaluation routine
c
c--------------------------
c input:
integer ict(3) ! selector:
c ict(1)=1 for f (don't evaluate if ict(1)=0)
c ict(2)=1 for df/dx (don't evaluate if ict(2)=0)
c
integer ivec ! vector dimensioning
c
c ivec-- number of vector pts (hermite values to look up)
c
real xvec(ivec) ! x-locations at which to evaluate
c
integer ivd ! 1st dimension of output array
c
c ivd -- 1st dimension of fval, .ge.ivec
c
c output:
real fval(ivd,*) ! output array
c
c fval(1:ivec,1) -- values as per 1st non-zero ict(...) element
c fval(1:ivec,2) -- values as per 2nd non-zero ict(...) element
c --etc--
c
c input:
integer nx ! dimension of hermite x grid
real xpkg(nx,4) ! x grid "package" (cf genxpkg)
real fspl(2,nx) ! Hermite coefficients (cf herm1ev)
c
c output:
c condition codes, 0 = normal return
integer iwarn ! =1 if an x value was out of range
integer ier ! =1 if argument error detected
c
subroutine vecherm2(ict,ivec,xvec,yvec,ivd,fval,
> nx,xpkg,ny,ypkg,fspl,inf2,
> iwarn,ier)
c
c vectorized hermite evaluation routine -- 2d
c 1. call vectorized zone lookup routine
c 2. call vectorized hermite evaluation routine
c
c--------------------------
c input:
integer ict(6) ! selector:
c ict(1)=1 for f (don't evaluate if ict(1)=0)
c ict(2)=1 for df/dx (don't evaluate if ict(2)=0)
c ict(3)=1 for df/dy (don't evaluate if ict(3)=0)
c ict(4)=1 for d2f/dxdy (don't evaluate if ict(4)=0)
c
integer ivec ! vector dimensioning
c
c ivec-- number of vector pts (hermite values to look up)
c
c list of (x,y) pairs:
c
real xvec(ivec) ! x-locations at which to evaluate
real yvec(ivec) ! y-locations at which to evaluate
c
integer ivd ! 1st dimension of output array
c
c ivd -- 1st dimension of fval, .ge.ivec
c
c output:
real fval(ivd,*) ! output array
c
c fval(1:ivec,1) -- values as per 1st non-zero ict(...) element
c fval(1:ivec,2) -- values as per 2nd non-zero ict(...) element
c --etc--
c
c input:
integer nx,ny ! dimension of hermite grids
real xpkg(nx,4) ! x grid "package" (cf genxpkg)
real ypkg(ny,4) ! y grid "package" (cf genxpkg)
integer inf2 ! fspl 3rd array dimension, .ge.nx
real fspl(0:3,inf2,ny) ! Hermite coefficients, cf herm2ev
c
c output:
c condition codes, 0 = normal return
integer iwarn ! =1 if an x value was out of range
integer ier ! =1 if argument error detected
subroutine vecherm3(ict,ivec,xvec,yvec,zvec,ivd,fval,
> nx,xpkg,ny,ypkg,nz,zpkg,fspl,inf2,inf3,
> iwarn,ier)
c
c vectorized hermite evaluation routine --
c 1. call vectorized zone lookup routine
c 2. call vectorized hermite evaluation routine
c
c--------------------------
c input:
integer ict(8) ! selector:
c ict(1)=1 for f (don't evaluate if ict(1)=0)
c ict(2)=1 for df/dx (don't evaluate if ict(2)=0)
c ict(3)=1 for df/dy (don't evaluate if ict(3)=0)
c ict(4)=1 for df/dy (don't evaluate if ict(4)=0)
c ict(5)=1 for d2f/dxdy (don't evaluate if ict(5)=0)
c ict(6)=1 for d2f/dxdz (don't evaluate if ict(6)=0)
c ict(7)=1 for d2f/dydz (don't evaluate if ict(7)=0)
c ict(8)=1 for d3f/dxdydz (don't evaluate if ict(8)=0)
c
integer ivec ! vector dimensioning
c
c ivec-- number of vector pts (hermite values to look up)
c
c list of (x,y,z) triples:
c
real xvec(ivec) ! x-locations at which to evaluate
real yvec(ivec) ! y-locations at which to evaluate
real zvec(ivec) ! z-locations at which to evaluate
c
integer ivd ! 1st dimension of output array
c
c ivd -- 1st dimension of fval, .ge.ivec
c
c output:
real fval(ivd,*) ! output array
c
c fval(1:ivec,1) -- values as per 1st non-zero ict(...) element
c fval(1:ivec,2) -- values as per 2nd non-zero ict(...) element
c --etc--
c
c input:
integer nx,ny,nz ! dimension of hermite grids
real xpkg(nx,4) ! x grid "package" (cf genxpkg)
real ypkg(ny,4) ! y grid "package" (cf genxpkg)
real zpkg(nz,4) ! z grid "package" (cf genxpkg)
integer inf2 ! fspl 4th array dimension, .ge.nx
integer inf3 ! fspl 5th array dimension, .ge.ny
real fspl(0:7,inf2,inf3,nz) ! Hermite coefficients cf herm3ev
c
c output:
c condition codes, 0 = normal return
integer iwarn ! =1 if an x value was out of range
integer ier ! =1 if argument error detected
c
subroutine gridherm1(x_newgrid,nx_new,f_new,nx,xpkg,fspl,
> iwarn,ier)
c
c regrid a hermite function f defined vs. x as in xpkg
c to a new grid, given by x_newgrid.
c
c set warning flag if the range x_newgrid exceeds the range of the
c original xpkg.
c
c (xpkg -- see genxpkg subroutine)
c
c input:
c
real x_newgrid(nx_new) ! new grid
c
c output:
c
real f_new(nx_new) ! f evaluated on this grid
c
c input:
c
integer nx ! size of old grid
real xpkg(nx,4) ! old grid "package"
real fspl(nx,2) ! compact spline coefficients of f
c
c output:
c condition codes, =0 for normal exit
c
integer iwarn ! =1 if new grid points out of range
integer ier ! =1 if there is an argument error
c
c--------------------------------------------
c local
c
integer ict(2)
c
data ict/1,0/
c
c--------------------------------------------
c
call vecherm1(ict,nx_new,x_newgrid,nx_new,f_new,nx,xpkg,fspl,
> iwarn,ier)
c
return
end
subroutine gridherm2(
> x_newgrid,nx_new,y_newgrid,ny_new,f_new,if1,
> nx,xpkg,ny,ypkg,fspl,inf2,iwarn,ier)
c
c regrid a hermite function f defined vs. x,y as in xpkg,ypkg
c to a new grid, given by x_newgrid, y_newgrid
c
c set warning flag if the range x_newgrid, y_newgrid
c exceeds the range of the original xpkg,ypkg.
c
c (xpkg,ypkg -- axis data, see genxpkg subroutine)
c
c input:
c
real x_newgrid(nx_new) ! new x grid
real y_newgrid(ny_new) ! new y grid
c
c output:
c
integer if1 ! 1st dimension of f_new
real f_new(if1,ny_new) ! f evaluated on this grid
c
c input:
c
integer nx ! size of old grid
real xpkg(nx,4) ! old grid "package"
integer ny ! size of old grid
real ypkg(ny,4) ! old grid "package"
c
integer inf2 ! array dimension
real fspl(0:4,inf2,ny) ! hermite coefficients of f
c
c output:
c condition codes, =0 for normal exit
c
integer iwarn ! =1 if new grid points out of range
integer ier ! =1 if there is an argument error
c
subroutine gridherm3(
> x_newgrid,nx_new,y_newgrid,ny_new,z_newgrid,nz_new,
> f_new,if1,if2,
> nx,xpkg,ny,ypkg,nz,zpkg,fspl,inf2,inf3,iwarn,ier)
c
c regrid a hermite function f defined vs. x,y,z as in xpkg, etc.
c to a new grid, given by x_newgrid, y_newgrid, z_newgrid
c
c set warning flag if the range exceeds the range of the
c original x/y/zpkg's.
c
c (xpkg/ypkg/zpkg -- axis data, see genxpkg subroutine)
c
c input:
c
real x_newgrid(nx_new) ! new x grid
real y_newgrid(ny_new) ! new y grid
real z_newgrid(nz_new) ! new z grid
c
c output:
c
integer if1,if2 ! 1st dimensions of f_new
real f_new(if1,if2,nz_new) ! f evaluated on this grid
c
c input:
c
integer nx ! size of old grid
real xpkg(nx,4) ! old grid "package"
integer ny ! size of old grid
real ypkg(ny,4) ! old grid "package"
integer nz ! size of old grid
real zpkg(nz,4) ! old grid "package"
c
integer inf2,inf3 ! array dimensions
real fspl(0:7,inf2,inf3,nz) ! hermite coefficients of f
c
c output:
c condition codes, =0 for normal exit
c
integer iwarn ! =1 if new grid points out of range
integer ier ! =1 if there is an argument error
C evaluate C1 cubic Hermite function interpolation -- 1d fcn
C --vectorized-- dmc 10 Feb 1999
C
subroutine herm1fcn(ict,ivec,ivecd,
> fval,ii,xparam,hx,hxi,fin,nx)
C
C input:
C
integer ict(2) ! requested output control
integer ivec ! vector length
integer ivecd ! vector dimension (1st dim of fval)
C
integer ii(ivec) ! target cells
real xparam(ivec)
! normalized displacements in cells
C
real hx(ivec) ! grid spacing, and
real hxi(ivec) ! inverse grid spacing 1/(x(i+1)-x(i))
C
real fin(0:1,nx) ! Hermite dataset
C
C output:
C
real fval(ivecd,2) ! interpolation results
C
C for detailed description of fin, ict and fval see subroutine
C herm1ev comments. Note ict is not vectorized -- the same output
C is expected to be returned for all input vector data points.
C
C note that the index inputs ii,jj and parameter inputs
C xparam,yparam,hx,hxi,hy,hyi are vectorized, and the
C
C output array fval has a vector ** 1st dimension ** whose
C size must be given as a separate argument; ivecd.ge.ivec
C expected!
C
C to use this routine in scalar mode, pass in ivec=ivecd=1
C
C evaluate C1 cubic Hermite function interpolation -- 2d fcn
C --vectorized-- dmc 10 Feb 1999
C
subroutine herm2fcn(ict,ivec,ivecd,
> fval,ii,jj,xparam,yparam,hx,hxi,hy,hyi,
> fin,inf2,ny)
C
integer ict(4) ! requested output selector
integer ivec ! vector length
integer ivecd ! vector dimension (1st dim of fval)
C
integer ii(ivec),jj(ivec) ! target cells (i,j)
real xparam(ivec),yparam(ivec)
! normalized displacements from (i,j) corners
C
real hx(ivec),hy(ivec) ! grid spacing, and
real hxi(ivec),hyi(ivec) ! inverse grid spacing 1/(x(i+1)-x(i))
! & 1/(y(j+1)-y(j))
C
real fin(0:3,inf2,ny) ! interpolant data (cf "herm2ev")
C ! f,fx,fy,fxy at each grid point.
C
real fval(ivecd,4) ! output returned
C ! function value, derivatives...
C
C for detailed description of fin, ict and fval see subroutine
C herm2ev comments. Note ict is not vectorized; the same output
C is expected to be returned for all input vector data points.
C
C note that the index inputs ii,jj and parameter inputs
C xparam,yparam,hx,hxi,hy,hyi are vectorized, and the
C output array fval has a vector ** 1st dimension ** whose
C size must be given as a separate argument
C
C to use this routine in scalar mode, pass in ivec=ivecd=1
C
C evaluate C1 cubic Hermite function interpolation -- 3d fcn
C --vectorized-- dmc 10 Feb 1999
C
subroutine herm3fcn(ict,ivec,ivecd,
> fval,ii,jj,kk,xparam,yparam,zparam,
> hx,hxi,hy,hyi,hz,hzi,
> fin,inf2,inf3,nz)
C
integer ict(8) ! requested output selector
integer ivec ! vector length
integer ivecd ! vector dimension (1st dim of fval)
C
integer ii(ivec),jj(ivec),kk(ivec) ! target cells (i,j,k)
real xparam(ivec),yparam(ivec),zparam(ivec)
! normalized displacements from (i,j,k) corners
C
real hx(ivec),hy(ivec),hz(ivec) ! grid spacing, and
real hxi(ivec),hyi(ivec),hzi(ivec) ! inverse grid spacing
! 1/(x(i+1)-x(i)) & 1/(y(j+1)-y(j)) & 1/(z(k+1)-z(i))
C
real fin(0:7,inf2,inf3,nz) ! interpolant data (cf "herm3ev")
C ! f,fx,fy,fz,fxy,fxz,fyz,fxyz at each
C ! grid point
C
real fval(ivecd,8) ! output returned
C ! function value, derivatives...
C
C for detailed description of fin, ict and fval see subroutine herm3ev
C comments. Note ict is not vectorized; the same output
C is expected to be returned for all input vector data points.
C
C note that the index inputs ii,jj,kk and parameter inputs
C xparam,yparam,zparam,hx,hxi,hy,hyi,hz,hzi are vectorized, and
C the output array fval has a vector ** 1st dimension ** whose
C size must be given as a separate argument
C
C to use this routine in scalar mode, pass in ivec=ivecd=1
C
subroutine herm1ev(xget,x,nx,ilinx,f,ict,fval,ier)
C
C evaluate a 1d cubic Hermite interpolant -- this is C1.
C
C input arguments:
C ================
C
real xget ! target of this interpolation
real x(nx) ! ordered x grid
integer ilinx ! ilinx=1 => assume x evenly spaced
C
real f(0:1,nx) ! function data
C
C contents of f:
C
C f(0,i) = f @ x(i)
C f(1,i) = df/dx @ x(i)
C
integer ict(2) ! code specifying output desired
C
C ict(1)=1 -- return f (0, don't)
C ict(2)=1 -- return df/dx (0, don't)
C
C output arguments:
C =================
C
real fval(2) ! output data
integer ier ! error code =0 ==> no error
C
C fval(1) receives the first output (depends on ict(...) spec)
C fval(2) receives the second output (depends on ict(...) spec)
C
C examples:
C on input ict = [1,1]
C on output fval= [f,df/dx]
C
C on input ict = [1,0]
C on output fval= [f] ... element 2 never referenced
C
C on input ict = [0,1]
C on output fval= [df/dx] ... element 2 never referenced
C
C ier -- completion code: 0 means OK
subroutine herm2ev(xget,yget,x,nx,y,ny,ilinx,iliny,
> f,inf2,ict,fval,ier)
C
C evaluate a 2d cubic Hermite interpolant on a rectilinear
C grid -- this is C1 in both directions.
C
C this subroutine calls two subroutines:
C herm2xy -- find cell containing (xget,yget)
C herm2fcn -- evaluate interpolant function and (optionally)
C the derivatives
C
C input arguments:
C ================
C
real xget,yget ! target of this interpolation
real x(nx) ! ordered x grid
real y(ny) ! ordered y grid
integer ilinx ! ilinx=1 => assume x evenly spaced
integer iliny ! iliny=1 => assume y evenly spaced
C
real f(0:3,inf2,ny) ! function data
C
C f 2nd dimension inf2 must be .ge. nx
C contents of f:
C
C f(0,i,j) = f @ x(i),y(j)
C f(1,i,j) = df/dx @ x(i),y(j)
C f(2,i,j) = df/dy @ x(i),y(j)
C f(3,i,j) = d2f/dxdy @ x(i),y(j)
C
integer ict(4) ! code specifying output desired
C
C ict(1)=1 -- return f (0, don't)
C ict(2)=1 -- return df/dx (0, don't)
C ict(3)=1 -- return df/dy (0, don't)
C ict(4)=1 -- return d2f/dxdy (0, don't)
C
C output arguments:
C =================
C
real fval(4) ! output data
integer ier ! error code =0 ==> no error
C
C fval(1) receives the first output (depends on ict(...) spec)
C fval(2) receives the second output (depends on ict(...) spec)
C fval(3) receives the third output (depends on ict(...) spec)
C fval(4) receives the fourth output (depends on ict(...) spec)
C
C examples:
C on input ict = [1,1,1,1]
C on output fval= [f,df/dx,df/dy,d2f/dxdy]
C
C on input ict = [1,0,0,0]
C on output fval= [f] ... elements 2 & 3 & 4 never referenced
C
C on input ict = [0,1,1,0]
C on output fval= [df/dx,df/dy] ... element 3 & 4 never referenced
C
C on input ict = [0,0,1,0]
C on output fval= [df/dy] ... elements 2 & 3 & 4 never referenced.
C
C ier -- completion code: 0 means OK
subroutine herm3ev(xget,yget,zget,x,nx,y,ny,z,nz,
> ilinx,iliny,ilinz,
> f,inf2,inf3,ict,fval,ier)
C
C evaluate a 3d cubic Hermite interpolant on a rectilinear
C grid -- this is C1 in all directions.
C
C this subroutine calls two subroutines:
C herm3xyz -- find cell containing (xget,yget,zget)
C herm3fcn -- evaluate interpolant function and (optionally)
C the derivatives
C
C input arguments:
C ================
C
real xget,yget,zget ! target of this interpolation
real x(nx) ! ordered x grid
real y(ny) ! ordered y grid
real z(nz) ! ordered z grid
integer ilinx ! ilinx=1 => assume x evenly spaced
integer iliny ! iliny=1 => assume y evenly spaced
integer ilinz ! ilinz=1 => assume z evenly spaced
C
real f(0:7,inf2,inf3,nz) ! function data
C
C f 2nd dimension inf2 must be .ge. nx; 3rd dim inf3 .ge. ny
C contents of f:
C
C f(0,i,j,k) = f @ x(i),y(j),z(k)
C f(1,i,j,k) = df/dx @ x(i),y(j),z(k)
C f(2,i,j,k) = df/dy @ x(i),y(j),z(k)
C f(3,i,j,k) = df/dz @ x(i),y(j),z(k)
C f(4,i,j,k) = d2f/dxdy @ x(i),y(j),z(k)
C f(5,i,j,k) = d2f/dxdz @ x(i),y(j),z(k)
C f(6,i,j,k) = d2f/dydz @ x(i),y(j),z(k)
C f(7,i,j,k) = d3f/dxdydz @ x(i),y(j),z(k)
C
integer ict(8) ! code specifying output desired
C
C ict(1)=1 -- return f (0, don't)
C ict(2)=1 -- return df/dx (0, don't)
C ict(3)=1 -- return df/dy (0, don't)
C ict(4)=1 -- return df/dz (0, don't)
C ict(5)=1 -- return d2f/dxdy (0, don't)
C ict(6)=1 -- return d2f/dxdz (0, don't)
C ict(7)=1 -- return d2f/dydz (0, don't)
C ict(8)=1 -- return d3f/dxdydz (0, don't)
C
C output arguments:
C =================
C
real fval(8) ! output data
integer ier ! error code =0 ==> no error
C
C fval(1) receives the first output (depends on ict(...) spec)
C fval(2) receives the second output (depends on ict(...) spec)
C fval(3) receives the third output (depends on ict(...) spec)
C fval(4) receives the 4th output (depends on ict(...) spec)
C fval(5-8) receive 5th thru 8th outputs (if req. by ict(...) spec)
C
C examples:
C on input ict = [1,1,1,1,0,0,0,0]
C on output fval= [f,df/dx,df/dy,df/dz]
C
C on input ict = [1,0,0,0,0,0,0,0]
C on output fval= [f] ... elements 2-8 never referenced
C
C on input ict = [0,1,1,0,0,0,0,0]
C on output fval= [df/dx,df/dy] ... elements 3-8 never referenced
C
C on input ict = [0,0,0,0,1,0,0,0]
C on output fval= [d2f/dxdy] ... elements 2-8 never referenced.
C
C ier -- completion code: 0 means OK
The pspline library also contains routines for 1d, 2d, and 3d
piecewise linear interpolation. Piecewise linear interpolation
is fast (most of the work is in zone lookup), and the user
need only supply data at the grid points; no spline or Hermite
coefficients are required. But, the interpolant, while
continuous, is not differentiable across grid zone boundaries.
There are no setup routines for the data itself (use of genxpkg
is still recommended for the grids).
Vectorized Evaluation Routines (recommended):
=============================================
given a location vector, do lookup & evaluation:
vecpc1 vecpc2 vecpc3
if zone location data is already computed, to do just
the vectorized interpolation evaluation:
pc1fcn pc2fcn pc3fcn
Scalar Evaluation Routines
==========================
pc1ev pc2ev pc3ev
subroutine vecpc1(ict,ivec,xvec,ivd,fval,nx,xpkg,fspl,
> iwarn,ier)
c
c vectorized piecewise linear evaluation routine -- 1d
c 1. call vectorized zone lookup routine
c 2. call vectorized piecewise linear evaluation routine
c
c note -- derivatives are *not* continuous
c
c--------------------------
c input:
integer ict(3) ! selector:
c ict(1)=1 for f (don't evaluate if ict(1)=0)
c ict(2)=1 for df/dx (don't evaluate if ict(2)=0)
c
integer ivec ! vector dimensioning
c
c ivec-- number of vector pts (piecewise linear values to look up)
c
real xvec(ivec) ! x-locations at which to evaluate
c
integer ivd ! 1st dimension of output array
c
c ivd -- 1st dimension of fval, .ge.ivec
c
c output:
real fval(ivd,*) ! output array
c
c fval(1:ivec,1) -- values as per 1st non-zero ict(...) element
c fval(1:ivec,2) -- values as per 2nd non-zero ict(...) element
c --etc--
c
c input:
integer nx ! dimension of x grid
real xpkg(nx,4) ! x grid "package" (cf genxpkg)
real fspl(nx) ! Piecewise Linear function
c
c output:
c condition codes, 0 = normal return
integer iwarn ! =1 if an x value was out of range
integer ier ! =1 if argument error detected
subroutine vecpc2(ict,ivec,xvec,yvec,ivd,fval,
> nx,xpkg,ny,ypkg,fspl,inf1,
> iwarn,ier)
c
c vectorized piecewise linear evaluation routine -- 2d
c 1. call vectorized zone lookup routine
c 2. call vectorized piecewise linear evaluation routine
c
c--------------------------
c input:
integer ict(6) ! selector:
c ict(1)=1 for f (don't evaluate if ict(1)=0)
c ict(2)=1 for df/dx (don't evaluate if ict(2)=0)
c ict(3)=1 for df/dy (don't evaluate if ict(3)=0)
c ict(4)=1 for d2f/dxdy (don't evaluate if ict(4)=0)
c
c note -- derivatives are *not* continuous!
c
integer ivec ! vector dimensioning
c
c ivec-- number of vector pts (piecewise linear values to look up)
c
c list of (x,y) pairs:
c
real xvec(ivec) ! x-locations at which to evaluate
real yvec(ivec) ! y-locations at which to evaluate
c
integer ivd ! 1st dimension of output array
c
c ivd -- 1st dimension of fval, .ge.ivec
c
c output:
real fval(ivd,*) ! output array
c
c fval(1:ivec,1) -- values as per 1st non-zero ict(...) element
c fval(1:ivec,2) -- values as per 2nd non-zero ict(...) element
c --etc--
c
c input:
integer nx,ny ! dimension of axis grids
real xpkg(nx,4) ! x grid "package" (cf genxpkg)
real ypkg(ny,4) ! y grid "package" (cf genxpkg)
integer inf1 ! fspl 3rd array dimension, .ge.nx
real fspl(inf1,ny) ! Piecewise Linear function
c
c output:
c condition codes, 0 = normal return
integer iwarn ! =1 if an x value was out of range
integer ier ! =1 if argument error detected
subroutine vecpc3(ict,ivec,xvec,yvec,zvec,ivd,fval,
> nx,xpkg,ny,ypkg,nz,zpkg,fspl,inf1,inf2,
> iwarn,ier)
c
c vectorized piecewise linear evaluation routine --
c 1. call vectorized zone lookup routine
c 2. call vectorized piecewise linear evaluation routine
c
c--------------------------
c input:
integer ict(8) ! selector:
c ict(1)=1 for f (don't evaluate if ict(1)=0)
c ict(2)=1 for df/dx (don't evaluate if ict(2)=0)
c ict(3)=1 for df/dy (don't evaluate if ict(3)=0)
c ict(4)=1 for df/dy (don't evaluate if ict(4)=0)
c ict(5)=1 for d2f/dxdy (don't evaluate if ict(5)=0)
c ict(6)=1 for d2f/dxdz (don't evaluate if ict(6)=0)
c ict(7)=1 for d2f/dydz (don't evaluate if ict(7)=0)
c ict(8)=1 for d3f/dxdydz (don't evaluate if ict(8)=0)
c
c note -- derivatives are *not* continuous!
c
integer ivec ! vector dimensioning
c
c ivec-- number of vector pts (piecewise linear values to look up)
c
c list of (x,y,z) triples:
c
real xvec(ivec) ! x-locations at which to evaluate
real yvec(ivec) ! y-locations at which to evaluate
real zvec(ivec) ! z-locations at which to evaluate
c
integer ivd ! 1st dimension of output array
c
c ivd -- 1st dimension of fval, .ge.ivec
c
c output:
real fval(ivd,*) ! output array
c
c fval(1:ivec,1) -- values as per 1st non-zero ict(...) element
c fval(1:ivec,2) -- values as per 2nd non-zero ict(...) element
c --etc--
c
c input:
integer nx,ny,nz ! dimension of axis grids
real xpkg(nx,4) ! x grid "package" (cf genxpkg)
real ypkg(ny,4) ! y grid "package" (cf genxpkg)
real zpkg(nz,4) ! z grid "package" (cf genxpkg)
integer inf1 ! fspl 1st array dimension, .ge.nx
integer inf2 ! fspl 2nd array dimension, .ge.ny
real fspl(0:7,inf1,inf2,nz) ! Piecewise Linear function array
c
c output:
c condition codes, 0 = normal return
integer iwarn ! =1 if an x value was out of range
integer ier ! =1 if argument error detected
c
C---------------------------------------------------------------------
C evaluate C0 piecewise linear function interpolation -- 1d fcn
C --vectorized-- dmc 10 Feb 1999
C
subroutine pc1fcn(ict,ivec,ivecd,
> fval,ii,xparam,hx,hxi,fin,nx)
C
C input:
C
integer ict(2) ! requested output control
integer ivec ! vector length
integer ivecd ! vector dimension (1st dim of fval)
C
integer ii(ivec) ! target cells
real xparam(ivec)
! normalized displacements in cells
C
real hx(ivec) ! grid spacing, and
real hxi(ivec) ! inverse grid spacing 1/(x(i+1)-x(i))
C
real fin(nx) ! the data
C
C output:
C
real fval(ivecd,2) ! interpolation results
C
C for detailed description of fin, ict and fval see subroutine
C pc1ev comments. Note ict is not vectorized -- the same output
C is expected to be returned for all input vector data points.
C
C note that the index inputs ii, and parameter inputs
C xparam,hx,hxi,are vectorized, and the
C
C output array fval has a vector ** 1st dimension ** whose
C size must be given as a separate argument; ivecd.ge.ivec
C expected!
C
C to use this routine in scalar mode, pass in ivec=ivecd=1
C
C---------------------------------------------------------------------
C evaluate piecewise bilinear function interpolation -- 2d fcn
C --vectorized-- dmc 10 Feb 1999
C
subroutine pc2fcn(ict,ivec,ivecd,
> fval,ii,jj,xparam,yparam,hx,hxi,hy,hyi,
> fin,inf2,ny)
C
integer ict(4) ! requested output control
integer ivec ! vector length
integer ivecd ! vector dimension (1st dim of fval)
C
integer ii(ivec),jj(ivec) ! target cells (i,j)
real xparam(ivec),yparam(ivec)
! normalized displacements from (i,j) corners
C
real hx(ivec),hy(ivec) ! grid spacing, and
real hxi(ivec),hyi(ivec) ! inverse grid spacing 1/(x(i+1)-x(i))
! & 1/(y(j+1)-y(j))
C
real fin(inf2,ny) ! interpolant data (cf "pc2ev")
C
real fval(ivecd,4) ! output returned
C
C for detailed description of fin, ict and fval see subroutine
C pc2ev comments. Note ict is not vectorized; the same output
C is expected to be returned for all input vector data points.
C
C note that the index inputs ii,jj and parameter inputs
C xparam,yparam,hx,hxi,hy,hyi are vectorized, and the
C output array fval has a vector ** 1st dimension ** whose
C size must be given as a separate argument
C
C to use this routine in scalar mode, pass in ivec=ivecd=1
C---------------------------------------------------------------------
C evaluate trilinear function interpolation -- 3d fcn
C --vectorized-- dmc 10 Feb 1999
C
subroutine pc3fcn(ict,ivec,ivecd,
> fval,ii,jj,kk,xparam,yparam,zparam,
> hx,hxi,hy,hyi,hz,hzi,
> fin,inf2,inf3,nz)
C
integer ict(8) ! requested output control
integer ivec ! vector length
integer ivecd ! vector dimension (1st dim of fval)
C
integer ii(ivec),jj(ivec),kk(ivec) ! target cells (i,j,k)
real xparam(ivec),yparam(ivec),zparam(ivec)
! normalized displacements from (i,j,k) corners
C
real hx(ivec),hy(ivec),hz(ivec) ! grid spacing, and
real hxi(ivec),hyi(ivec),hzi(ivec) ! inverse grid spacing
! 1/(x(i+1)-x(i)) & 1/(y(j+1)-y(j)) & 1/(z(k+1)-z(i))
C
real fin(inf2,inf3,nz) ! interpolant data (cf "pc3ev")
C
real fval(ivecd,8) ! output returned
C
C for detailed description of fin, ict and fval see subroutine pc3ev
C comments. Note ict is not vectorized; the same output
C is expected to be returned for all input vector data points.
C
C note that the index inputs ii,jj,kk and parameter inputs
C xparam,yparam,zparam,hx,hxi,hy,hyi,hz,hzi are vectorized, and the
C output array fval has a vector ** 1st dimension ** whose
C size must be given as a separate argument
C
C to use this routine in scalar mode, pass in ivec=ivecd=1
subroutine pc1ev(xget,x,nx,ilinx,f,ict,fval,ier)
C
C evaluate a 1d piecewise linear interpolant -- this is C0.
C ...a derivative can be evaluated but it is not continuous.
C
C this subroutine calls two subroutines:
C herm1x -- find cell containing (xget,yget)
C pc1fcn -- evaluate interpolant function and (optionally) derivatives
C
C input arguments:
C ================
C
real xget ! target of this interpolation
real x(nx) ! ordered x grid
integer ilinx ! ilinx=1 => assume x evenly spaced
C
real f(nx) ! function data
C
C contents of f:
C
C f(i) = f @ x(i)
C
integer ict(2) ! code specifying output desired
C
C ict(1)=1 -- return f (0, don't)
C ict(2)=1 -- return df/dx (0, don't)
C
C output arguments:
C =================
C
real fval(2) ! output data
integer ier ! error code =0 ==> no error
C
C fval(1) receives the first output (depends on ict(...) spec)
C fval(2) receives the second output (depends on ict(...) spec)
C
C examples:
C on input ict = [1,1]
C on output fval= [f,df/dx]
C
C on input ict = [1,0]
C on output fval= [f] ... element 2 never referenced
C
C on input ict = [0,1]
C on output fval= [df/dx] ... element 2 never referenced
C
C ier -- completion code: 0 means OK
subroutine pc2ev(xget,yget,x,nx,y,ny,ilinx,iliny,
> f,inf2,ict,fval,ier)
C
C evaluate a piecewise bilinear interpolant on a rectilinear
C grid -- this is C1 in both directions.
C
C this subroutine calls two subroutines:
C herm2xy -- find cell containing (xget,yget)
C pc2fcn -- evaluate interpolant function and (optionally) derivatives
C
C input arguments:
C ================
C
real xget,yget ! target of this interpolation
real x(nx) ! ordered x grid
real y(ny) ! ordered y grid
integer ilinx ! ilinx=1 => assume x evenly spaced
integer iliny ! iliny=1 => assume y evenly spaced
C
real f(inf2,ny) ! function data
C
C f 2nd dimension inf2 must be .ge. nx
C contents of f:
C
C f(i,j) = f @ x(i),y(j)
C
integer ict(4) ! code specifying output desired
C
C ict(1)=1 -- return f (0, don't)
C ict(2)=1 -- return df/dx (0, don't)
C ict(3)=1 -- return df/dy (0, don't)
C ict(4)=1 -- return d2f/dxdy (0, don't)
C
C output arguments:
C =================
C
real fval(4) ! output data
integer ier ! error code =0 ==> no error
C
C fval(1) receives the first output (depends on ict(...) spec)
C fval(2) receives the second output (depends on ict(...) spec)
C fval(3) receives the third output (depends on ict(...) spec)
C fval(4) receives the fourth output (depends on ict(...) spec)
C
C examples:
C on input ict = [1,1,1,1]
C on output fval= [f,df/dx,df/dy,d2f/dxdy]
C
C on input ict = [1,0,0,0]
C on output fval= [f] ... elements 2 & 3 & 4 never referenced
C
C on input ict = [0,1,1,0]
C on output fval= [df/dx,df/dy] ... element 3 & 4 never referenced
C
C on input ict = [0,0,1,0]
C on output fval= [df/dy] ... elements 2 & 3 & 4 never referenced.
C
C ier -- completion code: 0 means OK
subroutine pc3ev(xget,yget,zget,x,nx,y,ny,z,nz,
> ilinx,iliny,ilinz,
> f,inf2,inf3,ict,fval,ier)
C
C evaluate a trilinear interpolant on a rectilinear grid
C derivatives are available, but, not continuous across grid planes.
C
C this subroutine calls two subroutines:
C herm3xyz -- find cell containing (xget,yget,zget)
C pc3fcn -- evaluate interpolant function and (optionally) derivatives
C
C input arguments:
C ================
C
real xget,yget,zget ! target of this interpolation
real x(nx) ! ordered x grid
real y(ny) ! ordered y grid
real z(nz) ! ordered z grid
integer ilinx ! ilinx=1 => assume x evenly spaced
integer iliny ! iliny=1 => assume y evenly spaced
integer ilinz ! ilinz=1 => assume z evenly spaced
C
real f(inf2,inf3,nz) ! function data
C
C f 2nd dimension inf2 must be .ge. nx; 3rd dim inf3 .ge. ny
C contents of f:
C
C f(i,j,k) = f @ x(i),y(j),z(k)
C
integer ict(8) ! code specifying output desired
C
C ict(1)=1 -- return f (0, don't)
C ict(2)=1 -- return df/dx (0, don't)
C ict(3)=1 -- return df/dy (0, don't)
C ict(4)=1 -- return df/dz (0, don't)
C ict(5)=1 -- return d2f/dxdy (0, don't)
C ict(6)=1 -- return d2f/dxdz (0, don't)
C ict(7)=1 -- return d2f/dydz (0, don't)
C ict(8)=1 -- return d3f/dxdydz (0, don't)
C
C output arguments:
C =================
C
real fval(8) ! output data
integer ier ! error code =0 ==> no error
C
C fval(1) receives the first output (depends on ict(...) spec)
C fval(2) receives the second output (depends on ict(...) spec)
C fval(3) receives the third output (depends on ict(...) spec)
C fval(4) receives the 4th output (depends on ict(...) spec)
C fval(5-8) receive 5th thru 8th outputs (if required by ict(...) spec)
C
C examples:
C on input ict = [1,1,1,1,0,0,0,0]
C on output fval= [f,df/dx,df/dy,df/dz]
C
C on input ict = [1,0,0,0,0,0,0,0]
C on output fval= [f] ... elements 2-8 never referenced
C
C on input ict = [0,1,1,0,0,0,0,0]
C on output fval= [df/dx,df/dy] ... elements 3-8 never referenced
C
C on input ict = [0,0,0,0,1,0,0,0]
C on output fval= [d2f/dxdy] ... elements 2-8 never referenced.
C
C ier -- completion code: 0 means OK
The v_spline subroutine, by Wayne Houlberg, is used throughout
the pspline library to actually solve for 1d spline coefficients.
Caution -- the convention for representation of non-compact
spline coefficients is slightly different between v_spline and
the rest of pspline. In v_spline, the convention is
f(i,j) = (i-1)'th derivative of the spline at the left hand
side of interval j
and the formulae for evaluation, for x(j).le.x.le.x(j+1), is
delx=x-x(j)
s1(x) = f(1,j)+
delx*(f(2,j)+delx*(1.0/2.0)*(f(3,j)+delx*(1.0/3.0)*f(4,j)))
s1'(x) = f(2,j)+delx*(f(3,j)+delx*(1.0/2.0)*f(4,j))
s1''(x) = f(3,j)+delx*f(4,j)
s1'''(x) = f(4,j)
But, in the explicit spline coefficients used by pspline (e.g. the 1d
routine cspline), the f(i,j) is the coefficient of the (i-1)'th power
term in the cubic polynomial, i.e. the formula is
s2(x) = f(1,j)+
delx*(f(2,j)+delx*(f(3,j)+delx*f(4,j)))
s2'(x) = f(2,j)+2.0*delx*(f(3,j)+1.5*delx*f(4,j))
s2''(x) = 2.0*f(3,j) + 6.0*delx*f(4,j)
s2'''(x) = 6.0*f(4,j)
i.e. pspline's f(3,j) = (1.0/2.0) * v_spline's f(3,j)
pspline's f(4,j) = (1.0/6.0) * v_spline's f(4,j)
Bearing this in mind, users who want access to v_spline directly
may do so. The interface is as follows:
SUBROUTINE V_SPLINE(k_bc1,k_bcn,n,x,f,wk)
!***********************************************************************
!V_SPLINE evaluates the coefficients for a 1d cubic interpolating spline
!References:
! Forsythe, Malcolm, Moler, Computer Methods for Mathematical
! Computations, Prentice-Hall, 1977, p.76
! Engeln-Muellges, Uhlig, Numerical Algorithms with Fortran, Springer,
! 1996, p.251
! W.A.Houlberg, D.McCune 3/2000
!Input:
! k_bc1-option for BC at x(1)
! =-1 periodic, ignore k_bcn
! =0 not-a-knot
! =1 s'(x1) = input value of f(2,1)
! =2 s''(x1) = input value of f(3,1)
! =3 s'(x1) = 0.0
! =4 s''(x1) = 0.0
! =5 match first derivative to first 2 points
! =6 match second derivative to first 3 points
! =7 match third derivative to first 4 points
! =else use not-a-knot
! k_bcn-option for boundary condition at x(n)
! =0 not-a-knot
! =1 s'(x1) = input value of f(2,1)
! =2 s''(x1) = input value of f(3,1)
! =3 s'(x1) = 0.0
! =4 s''(x1) = 0.0
! =5 match first derivative to first 2 points
! =6 match second derivative to first 3 points
! =7 match third derivative to first 4 points
! =else use knot-a-knot
! n-number of data points or knots-(n.ge.2)
! x(n)-abscissas of the knots in strictly increasing order
! f(1,i)-ordinates of the knots
! f(2,1)-input value of s'(x1) for k_bc1=1
! f(2,n)-input value of s'(xn) for k_bcn=1
! f(3,1)-input value of s''(x1) for k_bc1=2
! f(3,n)-input value of s''(xn) for k_bcn=2
! wk(n)-scratch work area for periodic BC
!Output:
! f(2,i)=s'(x(i))
! f(3,i)=s''(x(i))
! f(4,i)=s'''(x(i))
[Added to pspline -- DMC April 2007]
The Pspline library now allows the definition of hybrid interpolating
functions for data sets of dimensionality 2 or higher.
In a hybrid interpolant, the interpolation method is controlled independently
along each dimension of the data.
In this context there are 4 interpolation methods that can be mixed:
-1: zonal step function
--no boundary condition
--number of data points = dimension grid size - 1
0: piecewise linear
--no boundary condition
--number of data points = dimension grid size
1: C1 Akima Hermite
--boundary condition: default, periodic, or control 1st derivative.
--number of data points = dimension grid size
2: C2 Compact Spline
--boundary condition: default, periodic, control 1st or 2nd derivative.
--number of data points = dimension grid size
These interpolation methods are specified separately for each dimension
of the data to be interpolated.
=> For technical reasons in the implementation, there is one restriction:
cubic interpolation methods (1) and (2) cannot be mixed.
Note that a zonal step function can be specified along 1 or more of the
data dimensions. But, if this is used, the size of the data along zonal
step function dimensions must be ONE LESS than the size of the corresponding
grid (as used e.g. in genxpkg), which specifies the zone boundaries.
Although hybrid setup and evaluation routines could be provided for datasets
of dimensionality greater than 3, this has not been done as of April 2007.
As with all other pspline library routines, the floating point arguments
of the routines named below are of precision REAL. To get routines with
arguments of precision REAL*8, prepend the prefix "R8" to the subroutine
name.
Routines for generating Hybrid data sets:
=========================================
for 2d data sets f(i,j) = f(x(i),y(j)): mkintrp2d
(real*8 version: r8mkintrp2d)
for 3d data sets f(i,j,k) = f(x(i),y(j),z(k)): mkintrp3d
(real*8 version: r8mkintrp3d)
Vectorized Evaluation Routines (recommended):
=============================================
given a location vector, do lookup & evaluation:
vecintrp2d vecintrp3d
given a grid (different than the original grid), to
map the function from the old grid to the new grid:
gridintrp2d gridintrp3d
if zone location data is already computed, to do just
the vectorized spline evaluation:
fvintrp2d fvintrp3d
Scalar Evaluation Routines
==========================
evintrp2d evintrp3d
Derivatives can be requested along any dimension, but, the value returned
will be ZERO if the derivative order exceeds the capability of the
interpolation method along that dimension.
Thus: any derivative along a zonal step function dimension results in ZERO.
A 2nd or higher derivative along a piecewise linear dimension yields ZERO.
A 2nd or higher derivative along a Hermite dimension yields ZERO.
subroutine mkintrp2d(x,nx,y,ny,jspline,
> f,icoeff,ixdim,iydim,
> ibcxmin,bcxmin,ibcxmax,bcxmax,
> ibcymin,bcymin,ibcymax,bcymax,
> ier)
C
C setup bicubic spline, or bicubic hermite, or Hybrid linear/zonal with
C 1d cubic Hermite or Spline interpolation in 1 dimension
C
implicit NONE
C
C input:
integer nx ! length of x vector
integer ny ! length of y vector
real x(nx) ! x vector, strict ascending
real y(ny) ! y vector, strict ascending
C
integer :: jspline(2) ! interpolation method control
C (1) -- 1st dimension: -1: zone, 0: pclin, 1: Hermite, 2: spline
C (2) -- 2nd dimension: -1: zone, 0: pclin, 1: Hermite, 2: spline
C
C Standard interpolation-- all jspline values match
C e.g. jspline(1)=jspline(2)=2 for bicubic spline
C Hybrid interpolation occurs when multiple jspline values are specified.
C
C RESTRICTION: jspline = (/ 1, 2 /) or (/ 2, 1 /) not allowed. I.e.
C Hermite and Spline interpolation cannot be mixed in a hybrid interpolant.
C This restriction exists because of technical issues in the
C implementation (it could be removed in principle but the work to do
C this has not been scheduled).
C
C coefficient buffer dimensioning:
integer :: icoeff ! #coefficients per data point
integer :: ixdim ! nx; nx-1 if jspline(1)==-1
integer :: iydim ! ny; ny-1 if jspline(2)==-1
C input/output:
real f(icoeff,ixdim,iydim) ! data & spline coefficients
C
C input bdy condition data:
C meaning of arguments as described in "mkbicubw.for"
C *ignored* for dimensions with piecewise linear or zonal interpolation;
C For Hermite interpolation only the values {-1,0,1} are accepted.
C
integer ibcxmin ! bc flag for x=xmin
real bcxmin(iydim) ! bc data vs. y at x=xmin
integer ibcxmax ! bc flag for x=xmax
real bcxmax(iydim) ! bc data vs. y at x=xmax
C
integer ibcymin ! bc flag for y=ymin
real bcymin(ixdim) ! bc data vs. x at y=ymin
integer ibcymax ! bc flag for y=ymax
real bcymax(ixdim) ! bc data vs. x at y=ymax
C
integer ier ! =0 on exit if there is no error.
c ier -- completion code, 0 for normal
subroutine mkintrp3d(x,nx,y,ny,z,nz,jspline,
> f,icoeff,ixdim,iydim,izdim,
> ibcxmin,bcxmin,ibcxmax,bcxmax,
> ibcymin,bcymin,ibcymax,bcymax,
> ibczmin,bczmin,ibczmax,bczmax,
> ier)
c
c setup a tricubic spline, or tricubic Hermite, or hybrid linear/zonal 2d or
c 1d with 1d or 2d cubic or Hermite spline interpolation
c
implicit NONE
C
C input:
integer nx ! length of x vector
integer ny ! length of y vector
integer nz ! length of z vector
real x(nx) ! x vector, strict ascending
real y(ny) ! y vector, strict ascending
real z(nz) ! z vector, strict ascending
c
integer :: jspline(3) ! interpolation method control
C (1) -- 1st dimension: -1: zone, 0: pclin, 1: Hermite, 2: spline
C (2) -- 2nd dimension: -1: zone, 0: pclin, 1: Hermite, 2: spline
C (3) -- 3rd dimension: -1: zone, 0: pclin, 1: Hermite, 2: spline
C
C Standard interpolation-- all jspline values match
C e.g. jspline(1)=jspline(2)=jspline(3)=2 for tricubic spline
C
C Hybrid interpolation-- not all jspline values are the same.
C
C RESTRICTION: if any jspline(...) element has value 1, none can have
C value 2. I.e. Spline and Hermite interpolation cannot currently be mixed.
C This restriction exists because of technical issues in the
C implementation (it could be removed in principle but the work to do
C this has not been scheduled).
c
c coefficient buffer dimensions
integer :: icoeff ! #coefficients per data point
integer :: ixdim ! nx; nx-1 if jspline(1)==-1
integer :: iydim ! ny; ny-1 if jspline(2)==-1
integer :: izdim ! nz; nz-1 if jspline(3)==-1
c input/output:
real f(icoeff,ixdim,iydim,izdim) ! data and spline coefficients
c
C boundary condition data
C Contiguous storage is assumed-- 1st dimension size must match
C actual use
C
C
integer ibcxmin,ibcxmax ! BC type flag @xmin, xmax
integer ibcymin,ibcymax ! BC type flag @ymin, ymax
integer ibczmin,ibczmax ! BC type flag @zmin, zmax
C
C
real bcxmin(iydim,izdim),bcxmax(iydim,izdim) ! xmin & xmax BC data
real bcymin(ixdim,izdim),bcymax(ixdim,izdim) ! ymin & ymax BC data
real bczmin(ixdim,iydim),bczmax(ixdim,iydim) ! zmin & zmax BC data
c
c where BC data is not required, dummy scalars may be passed.
C the ibc* flags determine whether BC data isneeded.
c
c BC data not required for zonal or piecewise linear interpolation
c for Hermite interpolation ibc* values from set {-1,0,1} are accepted.
c
c BC data: bcxmin & bcxmax: BC vs. y,z @xmin,xmax
C bcymin & bcymax: BC vs. x,z @ymin,ymax
C bczmin & bczmax: BC vs. x,y @zmin,zmax
C
c ibcxmin -- indicator for boundary condition at xmin=x(1):
c bcxmin(...) -- boundary condition data
c =-1 -- use periodic boundary condition
c =0 -- use "not a knot"
c =1 -- match slope, specified at x(1),y(iy),z(iz) by bcxmin(iy,iz)
c =2 -- match 2nd derivative, specified at x(1),y(iy),z(iz)
c by bcxmin(iy,iz
c =3 -- boundary condition is slope=0 (df/dx=0) at x(1), all y(j)
c =4 -- boundary condition is d2f/dx2=0 at x(1), all y(j)
c =5 -- df/dx BC from 1st divided difference
c =6 -- d2f/dx2 BC from 2nd divided difference (parabolic fit)
c =7 -- d3f/dx3 BC from 3rd divided difference (cubic fit)
c ***NOTE bcxmin(...) referenced ONLY if ibcxmin=1 or ibcxmin=2
c
c ibcxmax -- indicator for boundary condition at x(nx):
c bcxmax(...) -- boundary condition data
c (interpretation as with ibcxmin, bcxmin)
c NOTE: if ibcxmin=-1 then the periodic BC applies on both sides
c and ibcxmax, bcxmax are ignored.
c
c interpretation of ibcymin,bcymin,ibcymax,bcymax
c is same as with ibcxmin,...
c
c interpretation of ibczmin,bczmin,ibczmax,bczmax
c is same as with ibcxmin,...
c
c the explicit bdy condition arrays are referenced only if the
c corresponding "ibc" flag values are set > 0.
c
integer ier ! exit code
c ier -- completion code, 0 for normal
subroutine vecintrp2d(ict,ivec,xvec,yvec,ivd,fval,
> nx,xpkg,ny,ypkg,jspline,fspl,icoeff,ixdim,iydim,
> iwarn,ier)
c
implicit NONE
c
c vectorized hybrid spline evaluation routine -- 2d
c 1. call vectorized zone lookup routine
c 2. call vectorized hybrid spline evaluation routine
c
c--------------------------
c input:
integer ict(6) ! selector:
c ict(1)=1 for f (don't evaluate if ict(1)=0)
c ict(2)=1 for df/dx (don't evaluate if ict(2)=0)
c ict(3)=1 for df/dy (don't evaluate if ict(3)=0)
c ict(4)=1 for d2f/dx2 (don't evaluate if ict(4)=0)
c ict(5)=1 for d2f/dy2 (don't evaluate if ict(5)=0)
c ict(6)=1 for d2f/dxdy (don't evaluate if ict(6)=0)
c
c ict(1)=3 followed by ict(2:5) = 1 or 0 allow 3rd derivatives to
c be selected: fxxx fxxy fxyy fyyy
c
c ict(1)=4 followed by ict(2:4) allow 4th derivatives to
c be selected: fxxxy fxxyy fxyyy; fxxxx=fyyyy=0
c
c ict(1)=5 followed by ict(2:4) allow 4th derivatives to
c be selected: fxxxyy fxxyyy
c
c ict(1)=6 specifies 6th derivative: fxxxyyy (step fcn)
c
c in hybrid spline evaluation, any derivative along a dimension
c with zonal interpolation (jspline(i)=-1) gives zero;
c
c piecewise linear and Hermite interpolation give zero if a 2nd or
c higher derivative is sought along any dimension.
c
integer ivec ! vector dimensioning
c
c ivec-- number of vector pts (spline values to look up)
c
c list of (x,y) pairs:
c
real xvec(ivec) ! x-locations at which to evaluate
real yvec(ivec) ! y-locations at which to evaluate
c
integer ivd ! 1st dimension of output array
c
c ivd -- 1st dimension of fval, .ge.ivec
c
c output:
real fval(ivd,*) ! output array
c
c fval(1:ivec,1) -- values as per 1st non-zero ict(...) element
c fval(1:ivec,2) -- values as per 2nd non-zero ict(...) element
c --etc--
c
c input:
integer nx,ny ! dimension of spline grids
real xpkg(nx,4) ! x grid "package" (cf genxpkg)
real ypkg(ny,4) ! y grid "package" (cf genxpkg)
c
integer :: jspline(2) ! interpolation method on each dim:
c jspline(1) for x; jspline(2) for y
c -1: zonal step fcn; 0: pcwise linear; 1: Hermite; 2: compact spline
c
integer icoeff ! #coeffs per data node
integer ixdim ! x dim for fspl
integer iydim ! y dim for fspl
real fspl(icoeff,ixdim,iydim) ! hybrid spline coefficients
c
c ixdim=nx-1 for zonal step function along x dimension; o.w. expect ixdim=nx
c iydim=ny-1 for zonal step function along y dimension; o.w. expect iydim=ny
c
c output:
c condition codes, 0 = normal return
integer iwarn ! =1 if an x value was out of range
integer ier ! =1 if argument error detected
subroutine vecintrp3d(ict,ivec,xvec,yvec,zvec,ivd,fval,
> nx,xpkg,ny,ypkg,nz,zpkg,jspline,fspl,icoeff,ixdim,iydim,izdim,
> iwarn,ier)
c
implicit NONE
c
c vectorized hybrid spline evaluation routine -- 3d
c 1. call vectorized zone lookup routine
c 2. call vectorized hybrid spline evaluation routine
c
c--------------------------
c input:
integer ict(10) ! selector:
c ict(1)=1 for f (don't evaluate if ict(1)=0)
c ict(2)=1 for df/dx (don't evaluate if ict(2)=0)
c ict(3)=1 for df/dy (don't evaluate if ict(3)=0)
c ict(4)=1 for df/dy (don't evaluate if ict(4)=0)
c ict(5)=1 for d2f/dx2 (don't evaluate if ict(5)=0)
c ict(6)=1 for d2f/dy2 (don't evaluate if ict(6)=0)
c ict(7)=1 for d2f/dz2 (don't evaluate if ict(7)=0)
c ict(8)=1 for d2f/dxdy (don't evaluate if ict(8)=0)
c ict(9)=1 for d2f/dxdz (don't evaluate if ict(9)=0)
c ict(10)=1 for d2f/dydz (don't evaluate if ict(10)=0)
c
c higher derivatives can be selected by setting ict(1)=3,-3,4,-4,...
c see fvtricub.for comments.
c
c in hybrid spline evaluation, any derivative along a dimension
c with zonal interpolation (jspline(i)=-1) gives zero;
c
c piecewise linear and Hermite interpolation give zero if a 2nd or
c higher derivative is sought along any dimension.
c
integer ivec ! vector dimensioning
c
c ivec-- number of vector pts (spline values to look up)
c
c list of (x,y,z) triples:
c
real xvec(ivec) ! x-locations at which to evaluate
real yvec(ivec) ! y-locations at which to evaluate
real zvec(ivec) ! z-locations at which to evaluate
c
integer ivd ! 1st dimension of output array
c
c ivd -- 1st dimension of fval, .ge.ivec
c
c output:
real fval(ivd,*) ! output array
c
c fval(1:ivec,1) -- values as per 1st non-zero ict(...) element
c fval(1:ivec,2) -- values as per 2nd non-zero ict(...) element
c --etc--
c
c input:
integer nx,ny,nz ! dimension of spline grids
real xpkg(nx,4) ! x grid "package" (cf genxpkg)
real ypkg(ny,4) ! y grid "package" (cf genxpkg)
real zpkg(nz,4) ! z grid "package" (cf genxpkg)
c
integer :: jspline(3) ! interpolation method on each dim:
c jspline(1) for x; jspline(2) for y; jspline(3) for z
c -1: zonal step fcn; 0: pcwise linear; 1: Hermite; 2: compact spline
c
integer icoeff ! #coeffs per data node
integer ixdim ! x dim for fspl
integer iydim ! y dim for fspl
integer izdim ! z dim for fspl
real fspl(icoeff,ixdim,iydim,izdim) ! hybrid spline coefficients
c
c ixdim=nx-1 for zonal step function along x dimension; o.w. expect ixdim=nx
c iydim=ny-1 for zonal step function along y dimension; o.w. expect iydim=ny
c izdim=nz-1 for zonal step function along z dimension; o.w. expect izdim=nz
c
c output:
c condition codes, 0 = normal return
integer iwarn ! =1 if an x value was out of range
integer ier ! =1 if argument error detected
subroutine gridintrp2d(
> x_newgrid,nx_new,y_newgrid,ny_new,f_new,if1,
> nx,xpkg,ny,ypkg,jspline,fspl,icoeff,ixdim,iydim,iwarn,ier)
c
c regrid a spline function f defined vs. x,y as in xpkg,ypkg
c to a new grid, given by x_newgrid, y_newgrid
c
c set warning flag if the range x_newgrid, y_newgrid
c exceeds the range of the original xpkg,ypkg.
c
c (xpkg,ypkg -- axis data, see genxpkg subroutine)
c
c input:
c
real x_newgrid(nx_new) ! new x grid
real y_newgrid(ny_new) ! new y grid
c
c output:
c
integer if1 ! 1st dimension of f_new
real f_new(if1,ny_new) ! f evaluated on this grid
c
c input:
c
integer nx ! size of old grid
real xpkg(nx,4) ! old grid "package"
integer ny ! size of old grid
real ypkg(ny,4) ! old grid "package"
c
integer jspline(2) ! interpolation type by dimension
integer :: icoeff ! coefficients per data point
integer :: ixdim,iydim ! fspl dimensions
! =nx,ny unless zonal step function interpolation is used
real fspl(icoeff,ixdim,iydim) ! spline coefficients of f
c
c output:
c condition codes, =0 for normal exit
c
integer iwarn ! =1 if new grid points out of range
integer ier ! =1 if there is an argument error
subroutine gridintrp3d(
> x_newgrid,nx_new,y_newgrid,ny_new,z_newgrid,nz_new,
> f_new,if1,if2,
> nx,xpkg,ny,ypkg,nz,zpkg,
> jspline,fspl,icoeff,ixdim,iydim,izdim,
> iwarn,ier)
c
c regrid a spline function f defined vs. x,y,z as in xpkg, etc.
c to a new grid, given by x_newgrid, y_newgrid, z_newgrid
c
c set warning flag if the range exceeds the range of the
c original x/y/zpkg's.
c
c (xpkg/ypkg/zpkg -- axis data, see genxpkg subroutine)
c
c input:
c
real x_newgrid(nx_new) ! new x grid
real y_newgrid(ny_new) ! new y grid
real z_newgrid(nz_new) ! new z grid
c
c output:
c
integer if1,if2 ! 1st dimensions of f_new
real f_new(if1,if2,nz_new) ! f evaluated on this grid
c
c input:
c
integer nx ! size of old grid
real xpkg(nx,4) ! old grid "package"
integer ny ! size of old grid
real ypkg(ny,4) ! old grid "package"
integer nz ! size of old grid
real zpkg(nz,4) ! old grid "package"
c
integer :: jspline(3) ! interpolation type, by dimension
! -1: zonal step fcn; 0: pclin; 1: Hermite; 2: Spline
integer :: icoeff ! coefficients per data point
integer :: ixdim,iydim,izdim ! fspl dimensions
! ixdim=nx unless zonal step fcn interpolation in x is used
! similar comment applies to iydim,izdim
real fspl(icoeff,ixdim,iydim,izdim) ! spline coefficients of f
c
c output:
c condition codes, =0 for normal exit
c
integer iwarn ! =1 if new grid points out of range
integer ier ! =1 if there is an argument error
C evaluate Hybrid function interpolation -- 2d fcn
C NO ERROR CHECKING
C
subroutine fvintrp2d(ict,ivec,ivecd,
> fval,ii,jj,xparam,yparam,hx,hxi,hy,hyi,
> jspline,fin,icoeff,ixdim,iydim)
C
implicit NONE
C
integer ict(6) ! requested output control
integer ivec ! vector length
integer ivecd ! vector dimension (1st dim of fval)
C
integer ii(ivec),jj(ivec) ! target cells (i,j)
real xparam(ivec),yparam(ivec)
! normalized displacements from (i,j) corners
C
real hx(ivec),hy(ivec) ! grid spacing, and
real hxi(ivec),hyi(ivec) ! inverse grid spacing 1/(x(i+1)-x(i))
! & 1/(y(j+1)-y(j))
C
integer :: jspline(2) ! interpolation type control
integer :: icoeff ! coefficient dimension
integer :: ixdim,iydim ! x & y dimensions
real fin(icoeff,ixdim,iydim) ! data on which to interpolate...
C
real fval(ivecd,6) ! output returned
C
C for detailed description of fin, ict and fval see subroutine
C evintrp2d comments. Note ict is not vectorized; the same output
C is expected to be returned for all input vector data points.
C
C note that the index inputs ii,jj and parameter inputs
C xparam,yparam,hx,hxi,hy,hyi are vectorized, and the
C output array fval has a vector ** 1st dimension ** whose
C size must be given as a separate argument
C
C to use this routine in scalar mode, pass in ivec=ivecd=1
C
C evaluate Hybrid function interpolation -- 3d fcn
C NO ERROR CHECKING
C
subroutine fvintrp3d(ict,ivec,ivecd,
> fval,ii,jj,kk,xparam,yparam,zparam,
> hx,hxi,hy,hyi,hz,hzi,
> jspline,fin,icoeff,ixdim,iydim,izdim)
C
implicit NONE
C
C use mkintrp3d to set up spline coefficients...
C
integer ict(10) ! requested output control
integer ivec ! vector length
integer ivecd ! vector dimension (1st dim of fval)
C
integer ii(ivec),jj(ivec),kk(ivec) ! target cells (i,j,k)
real xparam(ivec),yparam(ivec),zparam(ivec)
! normalized displacements from (i,j,k) corners
C
real hx(ivec),hy(ivec),hz(ivec) ! grid spacing, and
real hxi(ivec),hyi(ivec),hzi(ivec) ! inverse grid spacing
! 1/(x(i+1)-x(i)) & 1/(y(j+1)-y(j)) & 1/(z(k+1)-z(i))
C
integer :: jspline(3) ! interpolation type control
integer :: icoeff ! coefficient dimension
integer :: ixdim,iydim,izdim ! x & y & z dimensions
real fin(icoeff,ixdim,iydim,izdim) ! data on which to interpolate...
C
real fval(ivecd,10) ! output returned
C
C for detailed description of fin, ict and fval see subroutine evintrp3d
C comments. Note ict is not vectorized; the same output
C is expected to be returned for all input vector data points.
C
C note that the index inputs ii,jj,kk and parameter inputs
C xparam,yparam,zparam,hx,hxi,hy,hyi,hz,hzi are vectorized, and the
C output array fval has a vector ** 1st dimension ** whose
C size must be given as a separate argument
C
C to use this routine in scalar mode, pass in ivec=ivecd=1
C
subroutine evintrp2d(xget,yget,x,nx,y,ny,jspline,
> f,icoeff,ixdim,iydim,ict,fval,ier)
C
implicit NONE
C
C evaluate a 2d Hybrid interpolant on a rectilinear grid
C
C this subroutine calls three subroutines:
C vecin2d_argchk -- error check
C herm2xy -- find cell containing (xget,yget)
C fvintrp2d -- evaluate interpolant function and (optionally) derivatives
C
C input arguments:
C ================
C
real xget,yget ! target of this interpolation
integer nx,ny ! grid sizes
real x(nx) ! ordered x grid
real y(ny) ! ordered y grid
C
integer :: jspline(2) ! interpolation method for each
C dimension: jspline(1) for x; jspline(2) for y
C =-1: zonal step function; =0: pc lin; =1: Hermite; =2: Spline
C
integer :: icoeff ! no. of coefficients per data point
integer :: ixdim,iydim ! dimensioning:
! ixdim=nx-1 for zonal step in x; otherwise nx
! iydim=ny-1 for zonal step in y; otherwise ny
real f(icoeff,ixdim,iydim) ! function data w/coefficients
C
integer ict(6) ! code specifying output desired
C
C Note on derivatives: for dimensions along which zonal step function
C interpolation is done, ANY derivative request returns ZERO.
C For dimensions along which piecewise linear or Hermite interpolation
C are done, more than one differentiation returns ZERO!
C
C Derivative controls are the same as for the compact 2d spline evaluation
C routine (evbicub):
C
C ict(1)=1 -- return f (0, don't)
C ict(2)=1 -- return df/dx (0, don't)
C ict(3)=1 -- return df/dy (0, don't)
C ict(4)=1 -- return d2f/dx2 (0, don't)
C ict(5)=1 -- return d2f/dy2 (0, don't)
C ict(6)=1 -- return d2f/dxdy (0, don't)
c the number of non zero values ict(1:6)
c determines the number of outputs...
c
c new dmc December 2005 -- access to higher derivatives (even if not
c continuous-- but can only go up to 3rd derivatives on any one coordinate.
c if ict(1)=3 -- want 3rd derivatives
c ict(2)=1 for d3f/dx3
c ict(3)=1 for d3f/dx2dy
c ict(4)=1 for d3f/dxdy2
c ict(5)=1 for d3f/dy3
c number of non-zero values ict(2:5) gives no. of outputs
c if ict(1)=4 -- want 4th derivatives
c ict(2)=1 for d4f/dx3dy
c ict(3)=1 for d4f/dx2dy2
c ict(4)=1 for d4f/dxdy3
c number of non-zero values ict(2:4) gives no. of outputs
c if ict(1)=5 -- want 5th derivatives
c ict(2)=1 for d5f/dx3dy2
c ict(3)=1 for d5f/dx2dy3
c number of non-zero values ict(2:3) gives no. of outputs
c if ict(1)=6 -- want 6th derivatives
c d6f/dx3dy3 -- one value is returned.
C
C output arguments:
C =================
C
real fval(6) ! output data
integer ier ! error code =0 ==> no error
C
C fval(1) receives the first output (depends on ict(...) spec)
C fval(2) receives the second output (depends on ict(...) spec)
C fval(3) receives the third output (depends on ict(...) spec)
C fval(4) receives the fourth output (depends on ict(...) spec)
C fval(5) receives the fifth output (depends on ict(...) spec)
C fval(6) receives the sixth output (depends on ict(...) spec)
C
C examples:
C on input ict = [1,1,1,0,0,1]
C on output fval= [f,df/dx,df/dy,d2f/dxdy], elements 5 & 6 not referenced.
C
C on input ict = [1,0,0,0,0,0]
C on output fval= [f] ... elements 2 -- 6 never referenced.
C
C on input ict = [0,0,0,1,1,0]
C on output fval= [d2f/dx2,d2f/dy2] ... elements 3 -- 6 never referenced.
C
C on input ict = [0,0,1,0,0,0]
C on output fval= [df/dy] ... elements 2 -- 6 never referenced.
C
C ier -- completion code: 0 means OK
subroutine evintrp3d(xget,yget,zget,x,nx,y,ny,z,nz,jspline,
> f,icoeff,ixdim,iydim,izdim,ict,fval,ier)
C
implicit NONE
C
C use mkintrp3d to set up spline coefficients...
C
C evaluate a 3d Hybrid interpolant on a rectilinear grid
C
C this subroutine calls two subroutines:
C vecin3d_argchk -- error check
C herm3xyz -- find cell containing (xget,yget,zget)
C fvintrp3d -- evaluate the spline function (w/derivatives if req.)
C
C input arguments:
C ================
C
real xget,yget,zget ! target of this interpolation
integer nx,ny,nz ! grid sizes
real x(nx) ! ordered x grid
real y(ny) ! ordered y grid
real z(nz) ! ordered z grid
C
integer :: jspline(3) ! interpolation method for each
C dimension: jspline(1) for x; jspline(2) for y; jspline(3) for z
C =-1: zonal step function; =0: pc lin; =1: Hermite; =2: Spline
C
integer :: icoeff ! no. of coefficients per data point
integer :: ixdim,iydim,izdim ! dimensioning:
! ixdim=nx-1 for zonal step in x; otherwise nx
! iydim=ny-1 for zonal step in y; otherwise ny
! izdim=nz-1 for zonal step in z; otherwise nz
C
real f(icoeff,ixdim,iydim,izdim) ! function data
C
integer ict(10) ! code specifying output desired
C
C Note on derivatives: for dimensions along which zonal step function
C interpolation is done, ANY derivative request returns ZERO.
C For dimensions along which piecewise linear or Hermite interpolation
C are done, more than one differentiation returns ZERO!
C
C Derivative controls are the same as for the compact 3d spline evaluation
C routine (evtricub):
C
C ict(1)=1 -- return f (0, don't)
C ict(2)=1 -- return df/dx (0, don't)
C ict(3)=1 -- return df/dy (0, don't)
C ict(4)=1 -- return df/dz (0, don't)
C ict(5)=1 -- return d2f/dx2 (0, don't)
C ict(6)=1 -- return d2f/dy2 (0, don't)
C ict(7)=1 -- return d2f/dz2 (0, don't)
C ict(8)=1 -- return d2f/dxdy (0, don't)
C ict(9)=1 -- return d2f/dxdz (0, don't)
C ict(10)=1-- return d2f/dydz (0, don't)
c
c (new dmc Dec 2005 -- higher derivatives available)
c ict(1)=3 --> 3rd derivative, .le.2 diff. in any coordinate
c ict(2:8) select: fxxy fxxz fxyy fxyz fxzz fyyz fyzz
c ->note ict(1)=3, ict(5)=1 gives fxyz = d3f/dxdydz
c ict(1)=-3 --> 3rd derivative, 3 in one coordinate
c ict(2:4) select: fxxx fyyy fzzz
c ict(1)=4 --> 3rd derivative, .le.2 diff. in any coordinate
c ict(2:7) select: fxxyy fxxyz fxxzz fxyyz fxyzz fyyzz
c ict(1)=-4 --> 3rd derivative, 3 in one coordinate
c ict(2:7) select: fxxxy fxxxz fxyyy fxzzz fyyyz fyzzz
c ict(1)=5 --> 3rd derivative, .le.2 diff. in any coordinate
c ict(2:4) select: fxxyyz fxxyzz fxyyzz
c ict(1)=-5 --> 3rd derivative, 3 in one coordinate
c ict(2:10) select: fxxxyy fxxxyz fxxxzz fxxyyy fxxzzz
c fxyyyz fxyzzz fyyyzz fzzzyy
c ict(1)=6 --> 3rd derivative, .le.2 diff. in any coordinate
c fxxyyzz
c ict(1)=-6 --> 3rd derivative, 3 in one coordinate
c ict(2:10) select: fxxxyyy fxxxyyz fxxxyzz fxxxyyz
c fxxyyyz fxxyzzz fxyyyzz fxyyzzz fyyyzzz
c ict(1)=-7 --> 7th derivative
c ict(2:7) select: fxxxyyyz fxxxyyzz fxxxyzzz
c fxxyyyzz fxxyyzzz fxyyyzzz
c ict(1)=-8 --> 8th derivative
c ict(2:4) select: fxxxyyyzz fxxxyyzzz fxxyyyzzz
c ict(1)=-9 --> 9th derivative: fxxxyyyzzz
c
C
C output arguments:
C =================
C
real fval(10) ! output data
integer ier ! error code =0 ==> no error
C
C fval(1) receives the first output (depends on ict(...) spec)
C fval(2) receives the second output (depends on ict(...) spec)
C fval(3) receives the third output (depends on ict(...) spec)
C fval(4) receives the 4th output (depends on ict(...) spec)
C fval(5-10) receive 5th thru 10th outputs (if required by ict(...) spec)
C
C examples:
C on input ict = [1,1,1,1,0,0,0,0,0,0,0]
C on output fval= [f,df/dx,df/dy,df/dz]
C
C on input ict = [1,0,0,0,0,0,0,0,0,0,0]
C on output fval= [f] ... elements 2-10 never referenced
C
C on input ict = [0,1,1,0,0,0,0,0,0,0,0]
C on output fval= [df/dx,df/dy] ... elements 3-10 never referenced
C
C on input ict = [0,0,0,0,1,0,0,0,0,0,0]
C on output fval= [d2f/dx2] ... elements 2-10 never referenced.
C
C ier -- completion code: 0 means OK
The PSPLINE software is distributed with a standalone program
`pspltest' which demonstrates the spline and hermite software
by using it to fit analytic expressions in 1d, 2d, and 3d, and
comparing the fit to the original. The location of the test
program will depend on site specific installation decisions.
TRANSP sites can find a copy of `pspltest' in the TRANSP
executables area.
`pspltest' runs the single precision test. For the R8 routines
there is a double precision test program, `pspltes8'.
Another standalone program, `qk_pspline', demonstrates the evaluation
of all available derivatives of splines, spot tests that expected
continuity conditions apply across cell boundaries for all function
values and derivatives, and that the R8 precision "compact" splines
(ezspline) and non-compact splines match (they should to near machine
precision).
Another pair of test programs exist for exercising the "ezspline"
interface. These are: ezspline_test_r4, ezspline_test_r8.
All the programs come with reference outputs. Comparing run outputs
with reference outputs takes some care. Some output quantities, e.g.
cpu timings, are bound to differ from machine to machine. And, binary
results will not be bit for bit identical, so, outputs of "residuals"
will show small differences from machine to machine.
This Document was created by hlptohtml
Written By:Manish Vachharajani(mvachhar@pppl.gov)