; NAME: 
;      SPLINESM
; 
; PURPOSE: 
;      To create smoothing spline for spectra or similar data and
;      return fitted values.
; 
; CALLING SEQUENCE: 
;      splinesm, x, y, dy, s, xfit, yfit, [help=help]
;
; INPUTS: 
;      x- 1 x N array of abscissa values to use in creating spline.
;      y- 1 x N array of ordinate values to use in creating spline.
;      dy- 1 x N array of error values in y to use in creating spline. 
;      s- Non-negative scalar which controls smoothing by limiting 
;         the sum of squares of diferences between calcuated and 
;         actual y values divided by error in y.  Set to zero will
;         cause the output to be equal to the intput. 
;         If s is negative SPLINESM will take absolute value.
;      xfit- 1 x M array of abscissa values to be used to get smoothed
;            ordinate values (yfit) from smoothing spline.
; 
; OPTIONAL KEYWORD INPUTS: 
;      help- Set to non-zero scalar to have SPLINESM print syntax.
; 
; OUTPUT:
;      yfit- Resulting values from application of smoothing spline
;            fit to x values.  -1 is returned if unable to calculate
; 
; EXTERNAL CALLS: 
;      NONE 
; 
; COMMON BLOCKS: 
;      NONE 
; 
; MODIFICATION HISTORY:
;      -Adapted to IDL by Eric Volquardsen, NASA IRTF, March 2003
;       from FORTRAN routine of same name by Dave Schleicher, Nov. 1993
; 
; ===========================================================================
PRO splinesm, x, y, dy, s, xfit, yfit, help=help

; CHECK INPUTS: 
     IF n_elements(help) EQ 1 THEN BEGIN
          IF help NE 0 THEN help = 1 
     ENDIF ELSE help = 0 
     IF n_params() NE 6 THEN GOTO, jumphp
     sx = size(x, /dimensions)  
     sy = size(y, /dimensions)
     sz = size(dy, /dimensions)
     sxf = size(xfit, /dimensions)
     IF n_elements(sx) NE 2 THEN help = 1 $ 
          ELSE IF sx[0] NE 1 THEN help = 1 
     IF n_elements(sy) NE 2 THEN help = 1 $ 
          ELSE IF sy[0] NE 1 THEN help = 1 
     IF n_elements(sz) NE 2 THEN help = 1 $ 
          ELSE IF sz[0] NE 1 THEN help = 1 
     IF n_elements(s) NE 1 THEN help = 1 ELSE IF s LT 0 THEN s = abs(s)
     IF help EQ 1 THEN BEGIN
          jumphp: 
          print, '****************************************************'
          print, '***             SPLINESM SYNTAX                  ***'
          print, '*** > splinesm, x, y, dy, s, xfit, yfit          ***'
          print, '***                                              ***'
          print, '*** x- 1 x N array of abscissa values            ***'
          print, '*** y- 1 x N array of ordinate values            ***'
          print, '*** dy- 1 x N array of error values in y         ***'
          print, '*** s- Non-negative scalar which controls        ***' 
          print, '***    smoothing by limiting the sum of squares  ***' 
          print, '***    of diferences between calcuated and       ***'
          print, '***    actual y values divided by error in y.    ***'
          print, '*** xfit- 1 x M array of output abscissa values  ***'
          print, '***       to use in calculating output ordinates ***' 
          print, '***       (yfit).                                ***'  
          print, '*** yfit- 1 x M array of fitted ordinate values  ***
          print, '***       (output).                              ***' 
          print, '****************************************************'
          RETURN 
     ENDIF
     sx = sx[1]
     IF sxf[0] NE 1 THEN sxf = sxf[0] ELSE sxf = sxf[1]

; INITIALIZE VALUES
     quitloop = 1
     m2 = sx - 2  ; was 3 greater (inconsistent w/size of vectors) 
     n = sx + 1 
     a = fltarr(n)
     b = fltarr(n)
     c = fltarr(n)
     d = fltarr(n)
     h = x[1] - x[0]
     f = (y[1] - y[0])/h
     t = fltarr(n)
     t1 = fltarr(n)
     r = fltarr(n)
     r1 = fltarr(n)
     r2 = fltarr(n)
     p = s   ; was set to zero (s didnt work)
     u = fltarr(n)
     v = fltarr(n) 
     FOR i = 2, m2 DO BEGIN  
          g = h 
          h = x[(i+1)] - x[i]
          e = f 
          f = (y[(i+1)] - y[i])/h
          a[i] = f - e 
          t[i] = 2*(g+h)/3
          t1[i] = h/3
          r2[i] = dy[(i-1)]/g
          r[i] = dy[(i+1)]/h
          r1[i] = (-1)*(dy[i]/g + dy[i]/h)
     ENDFOR

; CALCULATE FIT COEFFICIENTS
     FOR i = 0, m2 DO BEGIN   
          c[i] = r[i]*r1[(i+1)] + r1[i]*r2[(i+1)]
          d[i] = r[i]*r2[(i+2)]
     ENDFOR
     b = r^2 + r1^2 + r2^2
     f2 = (-1)*s
     gstart: 
     FOR i = 2, m2 DO BEGIN  
          r1[(i-1)] = f*r[(i-1)]
          r2[(i-2)] = g*r[(i-2)]
          r[i] = 1.0/(p*b[i] + t[i] - f*r1[(i-1)] - g*r2[(i-2)]) 
          u[i] = a[i] - r1[(i-1)]*u[(i-1)] - r2[(i-2)]*u[(i-2)]
          f = p*c[i] + t1[i] - h*r1[(i-1)]  
          g = h 
          h = d[i]*p 
     ENDFOR 

; CREATE U MATRIX
     FOR i = m2, 2, (-1) DO u[i] = r[i]*u[i] - r1[i]*u[(i+1)] - $ 
          r2[i]*u[(i+2)]
     e = 0.0
     h = 0.0 

; SMOOTHING PARAMETER CALCULATIONS 
     FOR i = 2, m2 DO BEGIN  
          g = h
          h = (u[(i+1)] - u[i]) / (x[(i+1)] - x[i])
          v[i] = (h-g)*dy[i]^2
          e = e + v[i]*(h-g)
     ENDFOR
     v[sx-1] = (-1)*h*dy[sx-1]^2
     g = v[sx]
     e = e - g*h
     g = f2
     f2 = e*p^2
     go = 0 
     IF f2 GE s THEN go = 1
     IF f2 LE g THEN go = 1 
     IF go EQ 1 THEN goto, gfinish
     f = 0.0 
     h = (v[2] - v[1]) / (x[2] - x[1])

     FOR i = 2, m2 DO BEGIN  
          g = h 
          h = (v[(i+1)] - v[i]) / (x[(i+1)] - x[i])
          g = h - g - r1[(i-1)]*r[(i-1)] - r2[(i-2)]*r[(i-2)]
          f = f + g*r[i]*g
          r[i] = g 
     ENDFOR
     h = e - p*f  
     IF h LE 0 THEN GOTO, gfinish
     p = p + (s-f2)/((sqrt(s/e) + p)*h)
     quitloop = quitloop + 1 
     IF quitloop EQ 100 THEN BEGIN
          yfit = -1
          print, 'UNABLE TO CALCULATE SPLINE'
          print, string(7b), format='(a,$)'
          wait, 0.25   
          print, string(7b), format='(a,$)'
          wait, 0.25   
          print, string(7b), format='(a,$)'
          RETURN 
     ENDIF
     GOTO, gstart

; FINAL COEFFICIENT ITERATIONS 
     gfinish:  
     a = y - p*v 
     c = u
     FOR i = 0, m2 DO BEGIN    ; ****
          h = x[(i+1)] - x[i]
          d[i] = (c[(i+1)] - c[i]) / 3*h
          b[i] = (a[(i+1)] - a[i]) / h - (h*d[i] + c[i]) * h 
     ENDFOR

; EVALUATE SPLINE AT NEEDED LOCATIONS
      ii = 0
      yfit = fltarr(sxf)
      FOR i = 0, (sxf-1) DO BEGIN
           IF ii GT sx THEN ii = 1
           IF xfit[i] LT x[ii] THEN GOTO, evg1
           IF xfit[i] LE x[(ii+1)] THEN GOTO, evg3
           evg1: 
           ii = 1 
           jj = sx
           evg2: 
           kk = (ii + jj)/2
           IF xfit[i] LT x[kk] THEN jj=kk
           IF xfit[i] GE x[kk] THEN ii=kk
           IF jj GT (ii+1) THEN GOTO, evg2 
           evg3: 
           dx = xfit[i] - x[ii]
           yfit[i] = a[ii] + dx*(b[ii] + dx*(c[ii] + dx*d[ii]))
      ENDFOR

END