; file: splinapp_nv ; init: Sep 18 2013 Rob Rutten, Deil ; last: Apr 11 2017 Rob Rutten Deil ; note: copy of splinesmooth.pro in otherlibs/vitas ;+ ; NAME: ; splinapp_nv ; ; PURPOSE: ; This procedure computes coefficients of cubic splines ; for a given observational set and smoothing parameter ; lambda. The method is coded according to Pollock D.S.G. ; (1999), "A Handbook of Time-Series Analysis, Signal ; Processing and Dynamics, Academic Press", San Diego ; ; CATEGORY: ; OAP ; ; CALLING SEQUENCE: ; SMTH = splinapp_nv([X,] Y, [SIGMA], LAMBDA=LAMBDA) ; ; INPUTS: ; X = 1D Array (independent variable) ; Y = 1D Array (function) ; SIGMA = 1D Array (weight of each measurement) By default ; all the measurements are of the same weight. ; ; KEYWORDS: ; LAMBDA = Smoothing parameter (It can be determined ; empirically, by the LS method or by cross- ; validation, eg. see book of Pollock.) LAMBDA ; equals 0 results in a cubic spline interpolation. ; In the other extreme, for a very large LAMBDA ; the result is smoothing by a linear function. ; ; EXAMPLE: ; Smth = SPLINESMOOTH(X, Y, LAMBDA = 1.d5) ; ; OUTPUT: ; SMTH: 1D array with smoothed values of a given function ; in given nodes. ; ; MODIFICATION HISTORY: ; Jan 2000 Nikola Vitas: start ; Mar 2007 Nikola Vitas: turned into function ; Sep 18 2013 RR: (..) to [..], result shifted 1 place in x ; Sep 6 2014 RR: corrected 1st value ;- FUNCTION splinapp_nv, XX, YY, SS, LAMBDA = LAMBDA CASE N_PARAMS() OF 1: BEGIN Y = XX X = INDGEN(N_ELEMENTS(Y)) SIGM = FLTARR(N_ELEMENTS(Y))+1 END 2: BEGIN Y = YY X = XX SIGM = FLTARR(N_ELEMENTS(Y))+1 END 3: BEGIN Y = YY X = XX SIGM = SS END ELSE: MESSAGE, 'Wrong number of arguments' ENDCASE NUM = SIZE(X, /N_ELEMENTS) N = NUM-1 IF NOT(KEYWORD_SET(LAMBDA)) THEN MESSAGE, 'Parameter lambda is not defined.' ; Definition of the help variables H = DBLARR(NUM) & R = DBLARR(NUM) & F = DBLARR(NUM) & P = DBLARR(NUM) Q = DBLARR(NUM) & U = DBLARR(NUM) & V = DBLARR(NUM) & W = DBLARR(NUM) ; Definition of the unknown coefficients A = DBLARR(NUM) & B = DBLARR(NUM) & C = DBLARR(NUM) & D = DBLARR(NUM) ; Computation of the starting values H[0] = X[1] - X[0] R[0] = 3.D/H[0] ; Computation of all H, R, F, P & Q FOR I = 1, N - 1 DO BEGIN H[I] = X[I+1] - X[I] R[I] = 3.D/H[I] F[I] = -(R[I-1] + R[I]) P[I] = 2.D * (X[I+1] - X[I-1]) Q[I] = 3.D * (Y[I+1] - Y[I])/H[I] - 3.D * (Y[I] - Y[I-1])/H[I-1] ENDFOR ; Compute diagonals of the matrix: W + LAMBDA T' SIGMA T FOR I = 1, N - 1 DO BEGIN U[I] = R[I-1]^2 * SIGM[I-1] + F[I]^2 * SIGM[I] + R[I]^2 * SIGM[I+1] U[I] = LAMBDA * U[I] + P[I] V[I] = F[I] * R[I] * SIGM[I] + R[I] * F[I+1] * SIGM[I+1] V[I] = LAMBDA * V[I] + H[I] W[I] = LAMBDA * R[I] * R[I+1] * SIGM[I+1] ENDFOR ; Decomposition in the form L' D L V[1] = V[1]/U[1] W[1] = W[1]/U[1] FOR J = 2, N-1 DO BEGIN U[J] = U[J] - U[J-2] * W[J-2]^2 - U[J-1] * V[J-1]^2 V[J] = (V[J] - U[J-1] * V[J-1] * W[J-1])/U[J] W[J] = W[J]/U[J] ENDFOR ; Gaussian eliminations to solve Lx = T'y Q[0] = 0.D FOR J = 2, N-1 DO Q[J] = Q[J] - V[J-1] * Q[J-1] - W[J-2] * Q[J-2] FOR J = 1, N-1 DO Q[J] = Q[J]/U[J] ; Gaussian eliminations to solve L'c = D^{-1}x Q[N-2] = Q[N-2] - V[N-2]*Q[N-1] FOR J = N-3, 1, -1 DO Q[J] = Q[J] - V[J] * Q[J+1] - W[J] * Q[J+2] ; Coefficients in the first segment D[0] = Y[0] - LAMBDA * R[0] * Q[1] * SIGM[0] ;RR corrected next line = email Nikola Apr 11 2017 D[1] = Y[1] - LAMBDA * (F[1] * Q[1] + R[1] * Q[2]) * SIGM[1] A[0] = Q[1]/(3.D * H[0]) B[0] = 0.D C[0] = (D[1] - D[0])/H[0] - Q[1] * H[0]/3.D ; Other coefficients FOR J = 1, N-1 DO BEGIN A[J] = (Q[J+1]-Q[J])/(3.D * H[J]) B[J] = Q[J] C[J] = (Q[J] + Q[J-1]) * H[J-1] + C[J-1] D[J] = R[J-1] * Q[J-1] + F[J] * Q[J] + R[J] * Q[J+1] D[J] = Y[J] - LAMBDA * D[J] * SIGM[J] ENDFOR D[N] = Y[N] - LAMBDA * R[N-1] * Q[N-1] * SIGM[N] Y1 = DBLARR(NUM) FOR I = 0, N_ELEMENTS(Y)-2 DO $ Y1[I+1] = $ ;RR changed [I] to [I+1] D[I] + C[I]*(X[I+1]-X[I]) + B[I]*(X[I+1]-X[I])^2 $ + A[I]*(X[I+1]-X[I])^3 ; ;RR added the second next line ;;; Y1[0] = D[0] + C[0]*(X[1]-X[0]) + B[0]*(X[1]-X[0])^2 + A[0]*(X[1]-X[0])^3 Y1[0] = D[0] Y1[N] = D[N] + C[N]*X[N] + B[N]*X[N]^2 + A[N]*X[N]^3 RETURN, Y1 END