ARJO SEGERS
KNMI
Date: 3 March 2009
subroutine Interval( x, a, ileft, iflag ) real, intent(in) :: x(:) real, intent(in) :: a integer, intent(inout) :: ileft integer, intent(out) :: iflag end subroutine Interval Computes: ileft = max( i, 1 <= i <= n, and x(i) < a ) Input: x: a real sequence of length n, assumed nondecreasing. a: the point whose location with tespect to the sequence x is to be determined. ileft: a 'first guess' estimate of the index to be found. this variable should be set to 1 for the first call to INTERVAL. Interpretate 'ileft' as the interval number {1,..,n-1} Output: ileft iflag meaning ------ -------- ------------------------------------- 1 -1 a < x(1) i 0 x(i) <= a < x(i+1) n-1 0 a = x( n ) n 1 x( n ) < a Method: Same as de Boor's INTERV sr, but not using a local index variable to store the last index found. Instead the 'first guess' index LEFT is used to initiate a search for the 'true' LEFT.
subroutine Swap( lx, xx ) integer, intent(in) :: lx real, intent(inout) :: xx(lx) end subroutine Swap
subroutine Interp_Lin( x,y, x0, y0, last, status ) real, intent(in) :: x(n), y(n) real, intent(in) :: x0 real, intent(out) :: y0 integer, intent(inout) :: last integer, intent(out) :: status end subroutine Interp_LinInterpolate 'y' defined at positions 'x' linear to 'x0', return result in 'y0'. Index 'last' is a help variable used to speedup search algorithm in repated calls to this routine.
subroutine Interp_Lin( x,y, x0, y0, last, status ) real, intent(in) :: x(n2), y(n1,n2) real, intent(in) :: x0 real, intent(out) :: y0(n1) integer, intent(inout) :: last integer, intent(out) :: status end subroutine Interp_Lin
subroutine Interp_Lin( x,y, x0, y0, last, status ) real, intent(in) :: x(n3), y(n1,n2,n3) real, intent(in) :: x0 real, intent(out) :: y0(n1,n2) integer, intent(inout) :: last integer, intent(out) :: status end subroutine Interp_Lin
subroutine Interp_Lin_Weight( x, x0, w, status ) real, intent(in) :: x(n) real, intent(in) :: x0 real, intent(out) :: w(n) integer, intent(out) :: status end subroutine InterpolWeightLet array 'y' be defined on axis points 'x', and let 'x0' be a point on the axis. The weights 'w' are then computed such that:
y interpolated to x0 = sum( y(:) * w(:) )
subroutine CircInterp_Lin( x, p, y, x0, y0, last, status ) real, intent(in) :: x(:), y(:) real, intent(in) :: p real, intent(in) :: x0 real, intent(out) :: y0 integer, intent(inout) :: last integer, intent(out) :: status end subroutine Interp_Lin
subroutine Interp_Lin( x,y, x0, y0, status ) real, intent(in) :: x(n), y(n) real, intent(in) :: x0(q) real, intent(out) :: y0(q) integer, intent(out) :: status end subroutine Interp_Lin
subroutine Interp_MuHerm( x,y, x0, y0, status ) real, intent(in) :: x(n), y(n) real, intent(in) :: x0(q) real, intent(out) :: y0(q) integer, intent(out) :: status end subroutine Interp_MuHerm
subroutine IntervalQuad_Lin( x,y, a,b, c, ilast, status ) real, intent(in) :: x(n), y(n) real, intent(in) :: a, b real, intent(out) :: c integer, intent(inout) :: ilast integer, intent(out) :: status end subroutine IntervalQuad_LinCompute an approximation to the integral:
b int y(x) dx x=aInput is formed by two vectors 'x' and 'y' . Vector 'y' contains the function values to the abscissa values contained in 'x'. The result, stored in 'c', is computed using the trapezoid formula throughout the domain of integration. If 'a < x(1)' or 'b > x(n)' then the function is extended to 'a', resp. 'b', assuming a constant value of 'y(1)', resp. 'y(n)'.
The index 'ilast' should be set to 1 in the first call. If subsequent integrals on the same vector are to be computed, it should be kept unaltered between calls. This will speed up the lookup of values if the various integrals are ascending in 'x'.
subroutine IntervalQuad_Cos_Lin( x, y, a, b, c, ilast, status ) real, intent(in) :: x(n), y(n) real, intent(in) :: a, b real, intent(out) :: c integer, intent(inout) :: ilast integer, intent(out) :: statusCompute an approximation to the integral:
b int y(x) cos(x) dx x=aSee routine 'IntervalQuad_Lin' for method.
subroutine IntervalQuad_Const( x, y, a,b, c, ilast, status ) real, intent(in) :: x(n+1), y(n) real, intent(in) :: a, b real, intent(out) :: c integer, intent(inout) :: ilast integer, intent(out) :: status end subroutine IntervalQuad_ConstCompute integral over '[a,b]' for function specified by:
y(1) y(2) y(3) +-----o-----+ |/ / / / / /+-----o----+ +----o----+ / / / / / / / / / /| | / / / / / / / / / / / | |/ / / / / / / / / / / /| ----+------|--+-----------+--------|-+--- x(1) x(2) x(3) x(4) a bSee routine 'IntervalQuad_Lin' for method.
subroutine IntervalSum( x, y, a_in, b_in, c, ilast, status, fac ) real, intent(in) :: x(n+1), y(n) real, intent(in) :: a_in, b_in real, intent(out) :: c integer, intent(inout) :: ilast integer, intent(inout) :: status real, intent(out), optional :: fac(n) end subroutine IntervalSumCompute sum of values 'y' over interval '[a,b]',
y(1) y(2) y(3) 100% +-----o-----+ 60% 30% | ^ +-----o----+ +----o----+ | ^ ^ | | | | | | | | ----+------|--+-----------+--------|-+--- x(1) x(2) x(3) x(4) a bEach 'y(i)' contributes to the sum for the fraction of '[x(i),x(i+1)]' covered by '[a,b]'. If optional output array 'fac' is provided it is filled with the contributions for each element of 'y'. The result is stored in 'c'. If 'a<x(1)' or 'b>x(n+1)' an error is issued.
subroutine Sort( x, xsort, isort, status ) real, intent(in) :: x (n) real, intent(out) :: xsort(n) integer, intent(out) :: isort(n) integer, intent(out) :: status end subroutine
Sort the values of an array 'x' in increasing order. Sorted 'x' is stored in 'xsort', indices in 'isort' such that:
xsort(i) = x(isort(i))
subroutine Sort( x, xsort, isort, status ) real, intent(in) :: x (n1,n2) real, intent(out) :: xsort(n1,n2) integer, intent(out) :: isort(n1,n2,2) integer, intent(out) :: status end subroutine Sort
Sort the values of an array 'x' in increasing order. Sorted 'x' is stored in 'xsort', indices in 'isort' such that:
xsort(i,j) = x(isort(i,j,1),isort(i,j,2))