diff options
Diffstat (limited to 'pkg/utilities/nttools/trebin/tucspl.f')
| -rw-r--r-- | pkg/utilities/nttools/trebin/tucspl.f | 52 |
1 files changed, 52 insertions, 0 deletions
diff --git a/pkg/utilities/nttools/trebin/tucspl.f b/pkg/utilities/nttools/trebin/tucspl.f new file mode 100644 index 00000000..511211b8 --- /dev/null +++ b/pkg/utilities/nttools/trebin/tucspl.f @@ -0,0 +1,52 @@ + subroutine tucspl (xa, ya, n, work, y2) +C +C Compute the second derivative of YA at each point in the array. +C This is the initialization needed in preparation for interpolating +C using cubic splines by the subroutine TUISPL. Input and output +C are all double precision. +C +C This routine was copied with slight modifications from the SPLINE +C subroutine in Numerical Recipes by Press, Flannery, Teukolsky and +C Vetterling. +C +C N i: number of elements in each array +C XA i: array of independent-variable values +C YA i: array of dependent-variable values +C WORK io: scratch array used for work space +C Y2 o: second derivative of YA at each point +C +CH Phil Hodge, 14-Apr-1988 Subroutine copied from Numerical Recipes SPLINE. +C + integer n + double precision xa(n), ya(n), work(n), y2(n) +C-- + integer i + double precision p, sig + +C These values (and y2(n) = 0) are for a "natural" spline. + y2(1) = 0. + work(1) = 0. +C +C This is the decomposition loop of the tridiagonal algorithm. +C Y2 and WORK are used for temporary storage of the decomposed factors. +C + do 10 i = 2, n-1 + sig = (xa(i) - xa(i-1)) / (xa(i+1) - xa(i-1)) + p = sig * y2(i-1) + 2. + y2(i) = (sig - 1.) / p + work(i) = (6. * ((ya(i+1) - ya(i)) / (xa(i+1) - xa(i)) + + - (ya(i) - ya(i-1)) / (xa(i) - xa(i-1))) / + + (xa(i+1) - xa(i-1)) - sig * work(i-1)) / p + 10 continue + +C "natural" spline + y2(n) = 0. +C +C This is the backsubstitution loop of the tridiagonal algorithm. +C + do 20 i = n-1, 1, -1 + y2(i) = y2(i) * y2(i+1) + work(i) + 20 continue + + return + end |