; NAME: 
;      Spectf
; 
; PURPOSE: 
;      Calculates shifts due to telescope flexure and extinction due to 
;      atmospheric water amounts when dividing object spectrum by a 
;      standard star spectrum.
; 
; CALLING SEQUENCE: 
;      Spectf, ast_spec, star_spec, ast_ext, star_ext, sig_list, order, $ 
;              thresh, band, fspec, Swater, Awater, Aoff, badp, help=help
; 
; INPUTS: 
;      ast_spec - Object reflectance spectrum
;      star_spec - Standard star reflectance spectrum
;      ast_ext - ATRAN model array for the object
;      star_ext - ATRAN model array for the standard star
;      sig_list - Errors for resulting spectrum
;      order - Polynomial fitting order to use for flagging bad points
;      thresh - Standard deviation threshold to use in finding bad points
;      band - Band number to fit (1=vis, 2=J, 3=K)
; 
; OPTIONAL KEYWORD INPUTS: 
;      help - Set to nonzero scalar to print syntax
; 
; OUTPUTS: 
;      fspec - Resulting spectrum
;      Swater - Calculated atmospheric water for standard star
;      Awater - Calculated atmospheric water for object
;      Aoff - Offset of object with respect to standard star
;      badp - Number of points flagged as bad
; 
; EXTERNAL CALLS: 
;      water_corr
; 
; MODIFICATION HISTORY: 
;      - Written by Bobby Bus, Univ. of Hawaii IfA, June 2002
;      - Combined Spect1, Spect2, Spect3; Renamed Spectf; added help keyword,
;        Eric Volquardsen, NASA IRTF, July 2003
; 
; =========================================================================

PRO Spectf, ast_spec, star_spec, ast_ext, star_ext, sig_list, order, $ 
                 thresh, band, fspec, Swater, Awater, Aoff, badp, 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_elements(ast_spec) LT 650 THEN help=1
     IF n_elements(star_spec) LT 650 THEN help=1
     IF n_elements(ast_ext) LT 650 THEN help=1
     IF n_elements(star_ext) LT 650 THEN help=1
     IF n_elements(sig_list) LT 80 THEN help=1
     IF n_elements(order) NE 1 THEN help=1
     IF n_elements(thresh) NE 1 THEN help=1
     IF n_elements(band) NE 1 THEN band=1
     IF help EQ 1 THEN BEGIN
          print, '*********************************************************'
          print, '***                   SPECTF SYNTAX                   ***'
          print, '*** > Spectf, ast_spec, star_spec, ast_ext, star_ext, ***' 
          print, '***           sig_list, order, thresh, band, fspec,   ***' 
          print, '***           Swater, Awater, Aoff, badp              ***' 
          print, '***                                                   ***' 
          print, '*** ast_spec - Object reflectance spectrum            ***' 
          print, '*** star_spec - Standard star reflectance spectrum    ***' 
          print, '*** ast_ext - ATRAN model array for the object        ***' 
          print, '*** star_ext - ATRAN model array for standard star    ***' 
          print, '*** sig_list - Errors for resulting spectrum          ***' 
          print, '*** order - Polynomial fitting order to use for       ***' 
          print, '***         flagging bad points                       ***' 
          print, '*** thresh - Standard deviation threshold to use in   ***' 
          print, '***          finding bad points                       ***' 
          print, '*** band - Band number to fit (1=vis, 2=J, 3=K)       ***' 
          print, '*** fspec - Resulting spectrum (output)               ***' 
          print, '*** Swater - Calculated atmospheric water for         ***' 
          print, '***          standard star (output)                   ***' 
          print, '*** Awater - Calculated atmospheric water for object  ***' 
          print, '***          (output)                                 ***' 
          print, '*** Aoff - Offset of object with respect to standard  ***' 
          print, '***        star (output)                              ***' 
          print, '*** badp - Number of points flagged as bad (output)   ***' 
          print, '*********************************************************' 
          fspec = 0 
          RETURN
     ENDIF
     CASE band OF 
          2 : wl = [80, 166]
          3 : wl = [160, 334]
          ELSE : BEGIN
               wl = [0, 86]
               band = 1 
          END
     ENDCASE

; INITIALIZE VALUES
     WSet, band
     arlen = wl[1]-wl[0]+1
     flags = Replicate(1,arlen)
     index = IndGen(arlen)
     wts = 1.0 / (sig_list)^3
     x_arr = (IndGen(15) * 0.001) - 0.007
     n_arr = FltArr(2,15)

; CALCULATE ATM WATER FOR STAR
     water_corr, star_spec, star_ext, band, starout, Swater
     wave_short = ast_spec[0,wl[0]:wl[1]]  

; ROUGH SHIFT AMOUNTS
     FOR j=0, 14 DO BEGIN

; Correct for water (amount not returned)
            water_corr, ast_spec, ast_ext, band, astout, offst=x_arr[j]

; Calculate goodness-of-fit
	    divout = astout[1,*] / starout[1,*]
	    b2 = svdfit(wave_short[*], divout[*], order, /DOUBLE, $ 
                        weights=wts, yfit=y_out)
            chi2 = total(wts * (divout - y_out)^2)

; Output for rough shifts
	    Plot, wave_short, divout, PSYM=4, /YNOZERO
            OPlot, wave_short, y_out
            print, x_arr[j], chi2
	    n_arr[0,j] = x_arr[j]
	    n_arr[1,j] = chi2
     ENDFOR

; Select best rough shift amount
     m_out = min(n_arr[1,*])

; FIRST ITERATION
; Generate new shifts based on rough results
     x_arr = (IndGen(11) * 0.0002) + (n_arr[0,!C] - 0.001)

; Reset goodness-of-fit array
     n_arr = FltArr(2,11)
     FOR j=0, 10 DO BEGIN

; Correct for water (amount not returned)
            water_corr, ast_spec, ast_ext, band, astout, offst=x_arr[j]

; Calculate goodness-of-fit
	    divout = astout[1,*] / starout[1,*]
	    b2 = svdfit(wave_short[*], divout[*], order, /DOUBLE, $ 
                        weights=wts, yfit=y_out)
            chi2 = total(wts * (divout - y_out)^2)

; Output for first iteration medium shifts
            Plot, wave_short, divout, PSYM=4, /YNOZERO
            OPlot, wave_short, y_out
            print, x_arr[j], chi2
	    n_arr[0,j] = x_arr[j]
	    n_arr[1,j] = chi2
     ENDFOR

; Select best shift amount
     m_out = min(n_arr[1,*])

; Generate fine shifts
     x_arr = (IndGen(5) * 0.0001) + (n_arr[0,!C] - 0.0002)

; Reset goodness-of-fit array
     n_arr = FltArr(2,5)
     FOR j=0, 4 DO BEGIN

; Correct for water (amount not returned)
            water_corr, ast_spec, ast_ext, band, astout, offst=x_arr[j]

; Calculate goodness-of-fit
	    divout = astout[1,*] / starout[1,*]
	    b2 = svdfit(wave_short[*], divout[*], order, /DOUBLE, $ 
                        weights=wts, yfit=y_out)
            chi2 = total(wts * (divout - y_out)^2)

; Output for first iteration fine shifts
            Plot, wave_short, divout, PSYM=4, /YNOZERO
            OPlot, wave_short, y_out
            print, x_arr[j], chi2
	    n_arr[0,j] = x_arr[j]
	    n_arr[1,j] = chi2
     ENDFOR

; Select best shift amount
     m_out = min(n_arr[1,*])
     Aoff = n_arr[0,!C]

; Get atmospheric water amount for best first iteration shift
     water_corr, ast_spec, ast_ext, band, astout, Awater, offst=Aoff

; Flag bad points and print number found in first iteration
     divout = astout[1,*] / starout[1,*]
     b2 = svdfit(wave_short[*], divout[*], order, /DOUBLE, weights=wts, $ 
                 yfit=y_out)
     diff1 = (divout - y_out) + (total(y_out) / arlen)
     sigma1 = stdev(diff1,mean)
     bad_list = where(abs(diff1 - mean) GT thresh*sigma1, bad_num)
     print, 'Number of bad points found in first iteration:', bad_num

; SECOND ITERATION
     IF bad_num GT 0 THEN BEGIN

; Remove bad points from iteration
          flags[bad_list] = -1
          good_list = where(flags GT 0)
          index1 = index[good_list]
          xx = wave_short[good_list]
          ee = wts[good_list]

; Perform medium shifts for second iteration
          x_arr = (IndGen(11) * 0.0002) + (Aoff - 0.001)
          n_arr = FltArr(2,11)
          FOR j=0, 10 DO BEGIN

; Correct for water (amount not returned)
               water_corr, ast_spec, ast_ext, band, astout, offst=x_arr[j]

; Calculate goodness-of-fit
               divout = astout[1,*] / starout[1,*]
               divshort = divout[good_list]
               b2 = svdfit(xx[*], divshort[*], order, /DOUBLE, weights=ee, $ 
                           yfit=y_out)
               chi2 = total(ee * (divshort - y_out)^2)

; Output of medium shifts in second iteration
               Plot, xx, divshort, PSYM=4, /YNOZERO
               OPlot, xx, y_out
               print, x_arr[j], chi2
               n_arr[0,j] = x_arr[j]
               n_arr[1,j] = chi2
          ENDFOR

; Select best shift amount from medium shifts second iteration
          m_out = min(n_arr[1,*])

; Perform fine shift for second iteration
          x_arr = (IndGen(5) * 0.0001) + (n_arr[0,!C] - 0.0002)
          n_arr = FltArr(2,5)
          FOR j=0, 4 DO BEGIN

; Correct for water (amount not returned)
               water_corr, ast_spec, ast_ext, band, astout, offst=x_arr[j]

; Calculate goodness-of-fit
               divout = astout[1,*] / starout[1,*]
               divshort = divout[good_list]
               b2 = svdfit(xx[*], divshort[*], order, /DOUBLE, weights=ee, $ 
                           yfit=y_out)
               chi2 = total(ee * (divshort - y_out)^2)

; Output of fine shifts for second iteration
               Plot, xx, divshort, PSYM=4, /YNOZERO
               OPlot, xx, y_out
               print, x_arr[j], chi2
               n_arr[0,j] = x_arr[j]
               n_arr[1,j] = chi2
          ENDFOR

; Select best shift amount for second iteration
          m_out = min(n_arr[1,*])
          Aoff = n_arr[0,!C]

; Get atmospheric water amount for best second iteration shift
          water_corr, ast_spec, ast_ext, band, astout, Awater, offst=Aoff

; Flag bad points and print number found in second iteration
          divout = astout[1,*] / starout[1,*]
          divshort = divout[good_list]
          b2 = svdfit(xx[*], divshort[*], order, /DOUBLE, weights=ee, $ 
                      yfit=y_out)
          diff = (divshort - y_out) + (total(y_out) / (arlen - bad_num))
          sigma = stdev(diff,mean)
          bad_list = where(abs(diff - mean) GT thresh*sigma, bad_num)
          print, 'Number of bad points found in second iteration:', bad_num

; THIRD ITERATION
          IF bad_num GT 0 THEN BEGIN

; Remove bad points from iteration
               bad = index1[bad_list]
               flags[bad] = -1
               good_list = where(flags GT 0)
               index2 = index[good_list]
               xx = wave_short[good_list]
               ee = wts[good_list]

; Perform medium shifts for third iteration
               x_arr = (IndGen(11) * 0.0002) + (Aoff - 0.001)
               n_arr = FltArr(2,11)
               FOR j=0, 10 DO BEGIN

; Correct for water (amount not returned)
                    water_corr, ast_spec, ast_ext, band, astout, offst=x_arr[j]

; Calculate goodness-of-fit
                    divout = astout[1,*] / starout[1,*]
                    divshort = divout[good_list]
                    b2 = svdfit(xx[*], divshort[*], order, /DOUBLE, $ 
                                weights=ee, yfit=y_out)
                    chi2 = total(ee * (divshort - y_out)^2)

; Output of medium shifts in third iteration
                    Plot, xx, divshort, PSYM=4, /YNOZERO
                    OPlot, xx, y_out
                    print, x_arr[j], chi2
                    n_arr[0,j] = x_arr[j]
                    n_arr[1,j] = chi2
               ENDFOR

; Select best shift amount from medium shifts third iteration
               m_out = min(n_arr[1,*])

; Perform fine shift for third iteration
               x_arr = (IndGen(5) * 0.0001) + (n_arr[0,!C] - 0.0002)
               n_arr = FltArr(2,5)
               FOR j=0, 4 DO BEGIN

; Correct for water (amount not returned)
                    water_corr, ast_spec, ast_ext, band, astout, offst=x_arr[j]

; Calculate goodness-of-fit
                    divout = astout[1,*] / starout[1,*]
                    divshort = divout[good_list]
                    b2 = svdfit(xx[*], divshort[*], order, /DOUBLE, $ 
                                weights=ee, yfit=y_out)
                    chi2 = total(ee * (divshort - y_out)^2)

; Output of fine shifts in third iteration
                    Plot, xx, divshort, PSYM=4, /YNOZERO
                    OPlot, xx, y_out
                    print, x_arr[j], chi2
                    n_arr[0,j] = x_arr[j]
                    n_arr[1,j] = chi2
               ENDFOR

; Select best shift amount from third iteration
               m_out = min(n_arr[1,*])
               Aoff = n_arr[0,!C]

; Get atmospheric water amount for best third iteration shift
               water_corr, ast_spec, ast_ext, band, astout, Awater, offst=Aoff

; Calculate goodness-of-fit
               divout = astout[1,*] / starout[1,*]
               divshort = divout[good_list]
               b2 = svdfit(xx[*], divshort[*], order, /DOUBLE, weights=ee, $ 
                           yfit=y_out)
          ENDIF

; Flag bad points and print number found in third iteration
          fit_out = Interpol(y_out, xx, wave_short, /SPLINE)
          y_out = fit_out
          diff = (divout - y_out) + (total(y_out) / arlen)
          fmean = mean(diff)
          bad_list = where(abs(diff - fmean) GT thresh*sigma, bad_num)
          print, 'Number of bad points found in final pass:', bad_num
          badp = bad_num
          flags = Replicate(1,arlen)
          flags[bad_list] = -1
     ENDIF

; FINAL OUTPUT
     print, n_arr[0,!C], n_arr[1,!C]
     fspec = FltArr(4,arlen)
     fspec[0,*] = wave_short
     fspec[1,*] = divout
     fspec[2,*] = sig_list
     fspec[3,*] = flags
     good_pts = where(flags GT 0)
     Plot, fspec[0,good_pts], fspec[1,good_pts], PSYM=4, /YNOZERO
     OPlot, fspec[0,*], y_out

END