Previous: fitsf Up: ../plot79_f.html Next: fittap
SUBROUTINE FITSM (X,Y, DYINV2, N, S, A,B,C,D, R,R1,R2, T,T1, U,V) C$ (Cubic Spline with Smoothing) C$ This subroutine implements the evaluation of a cubic spline C$ with smoothing through a set of data points (X(K),Y(K)) in C$ which the ordinate values Y(K) have uncertainties C$ associated with them. The resulting spline will generally C$ not go through the input points, but will provide a C$ smoothed curve close to them. C$ C$ The input arguments are: C$ C$ X(*)................Array of abscissa values, arranged in C$ strictly increasing order. C$ Y(*)................Array of ordinate values. C$ DYINV2(*)...........Array of the inverse square of the C$ uncertainties DY(*) in the values Y(*). C$ If available, the values DY(K) should C$ be standard deviations of the points C$ Y(K). C$ N...................Number of data points. C$ S...................A non-negative parameter which controls C$ the extent of smoothing. the spline C$ function is determined such that C$ C$ N C$ SUM ((F(X(K))-Y(K))/DY(K))**2 .LE. S. C$ K=1 C$ C$ where equality holds unless F describes C$ a straight line. If S = 0.0, a normal C$ cubic spline passing exactly through C$ the input points will be produced. A C$ natural choice for S will often be (N - C$ SQRT(2N)) .LE. S .LE. (N + SQRT(2N)). C$ C$ The output arguments are: C$ C$ A(*),B(*),C(*),D(*)......Arrays of dimension N collecting C$ the coefficients of the cubic C$ spline F(T) such that with X(K) C$ .LE. T .LE. X(K+1) and H = T - C$ X(K), we have F(T) = ((D(K)*H + C$ C(K))*H + B(K))*H + A(K). C$ Furthermore, A(N) = F(X(N)) and C$ C(N)=0, while B(N) and D(N) are C$ left undefined. C$ C$ The remaining arguments R(*), R1(*), R2(*), T(*), T1(*), C$ U(*) and V(*) are used for scratch space, and all have C$ dimension N+2. C$ C$ The method has been published by C$ C$ C. Reinsch, "Smoothing by Spline Functions", Numerische C$ Mathematik 10, 177-183 (1967). C$ C$ This routine is a FORTRAN version of Reinsch's Algol C$ procedure. Reinsch's published Algol procedure did not C$ initialize A(N) and C(N) to the stated values. This has C$ been corrected in this version. C$ (03-APR-82)