pro fitspec, wl, flux, eflux,  $
             chisqarr, sptarr, extarr, dof,  $
             Avmax = Avmax, $
             wlrange = wlrange, $
             display = display, $
             nodisp = nodisp,  $
             outfile = outfile


;+
; given a spectrum, loop through calibrator sample and fit spectrum
; for the spectral type and the reddending
;
; uses "dwarfspectra-fit.lst" for list of spectra
;
; 07/26/01 MCL
;
; added 'Avmax' parm
; for fits w/o extinction, will plot chisq vs spt and find the minima
; added /display: examine best fit for each calibrator
; can write to output file
; 09/12/01 MCL
;-

if n_params() lt 2 then begin
    print, 'pro fitspec, wl, flux, eflux,'
    print, '             chisqarr, sptarr, extarr, dof, '
    print, '             [Avmax=], '
    print, '             [wlrange=], '
    print, '             [display]'
    print, '             [nodisp]'
    print, '             [outfile=]'
    return
endif


;; range of wl's to fit
;wlrange = [1.95, 2.4]
if not(keyword_set(wlrange)) then wlrange = limits(wl)

; list of reddenings to explore
if (n_elements(Avmax) eq 0) then Avmax = 0.0
Av_list = makex(0, Avmax, 1)


; get list of calibrator spectra
; sort by spectral type
readcol2, 'dwarfspectra-fit.lst', speclist, sptname, ref, units, specdir,  $
  form = 'a,a,a,a,a', /sil
spt = sptypenum(sptname)
ss = sort(spt)
speclist = speclist(ss)
sptname = sptname(ss)
spt = spt(ss)
ref = ref(ss)
units = units(ss)
specdir = specdir(ss)


; initialize
nspec = n_elements(speclist)
next = n_elements(Av_list)
chisqarr = fltarr(nspec, next)
sptarr = fltarr(nspec, next)
extarr = fltarr(nspec, next)


; normalize object spectrum 
w = where(wl ge wlrange(0) and wl le wlrange(1), npix)
wlobj = wl(w)
fobj = flux(w)/median(flux(w))
efobj = eflux(w)/median(flux(w))
dof = npix-2  ; degrees of freedom


; loop over calibrator stars
print, form = '($,A)', 'looping over calibrators: '
;for i = 0, nspec-1, 6 do begin
if keyword_set(display) then begin
    win, 0, 500
    lincolr
endif
for i = 0, nspec-1 do begin

    ; load calibrator spectrum and shift to object wl scale
    readcol2, /silent, specdir(i)+speclist(i), wlcal0, fcal0, form = 'f,f'
    print, form = '($,A," ")', sptname(i)
    fcal = interpol(fcal0, wlcal0, wlobj)
    fcal = fcal/median(fcal)

    ;plot, wlobj, fobj, chars = 2, ps = 10, /xs
    ;oplot, wlobj, fcal, col = 2, ps = 10
    ;if not(getyn('ok?')) then stop

    ; loop over extinctions
    for j = 0, next-1 do begin
        ; store the SpT/extinction combination
        sptarr(i, j) = spt(i)
        extarr(i, j) = Av_list(j)
        ; extinct
        fcalext = fcal * mag(0.412 * wlobj^(-1.75) * Av_list(j), /rev)
        fcalext = fcalext/median(fcalext)
        ; compute g.o.f.
        chisqarr(i, j) = total( ((fcalext-fobj)/efobj)^2. )
    endfor

    ; display best fit for each calibrator
    if keyword_set(display) then begin
        whereismin, chisqarr(i, *), wm, y, /silent
        plot, wlobj, fobj, /xs, $
          tit = 'SpT = '+strc(spt(i))+', Av = '+strc(Av_list(y)), chars = 2
        oplot, wlobj, fcal * mag(0.412 * wlobj^(-1.75) * Av_list(wm), /rev), col = 2
        wait, 0.5
    endif

endfor
print

if not(keyword_set(nodisp)) or keyword_set(outfile) then begin

    ; open output device
    if keyword_set(outfile) then begin
        fdecomp, outfile, disk, dir, name, ext
        ps_open, name, thick = 5, /color
    endif else $
      win, 0, 500

    ; 2-d fit: contour plot
    if n_elements(Av_list) gt 1 then begin
        loadct, 15
        display2, chisqarr, spt, Av_list, $
          tit = 'reduced chi-square', chars = 1.5,  $
          xtit = 'sp type', ytit = 'A(V)' 
    ; 1-d fit
    endif else begin
        ; plot chi-sq vs SpT
        plot, spt, chisqarr, ps = 8, $
          xtit = 'sp type', ytit = 'reduced chi square', chars = 2
        cc = poly_fit(spt, chisqarr, 2, yfit)
        oplot, spt, yfit

        sptbest = -cc(1)/cc(2)/2.0
        plots, [0, 0]+sptbest, !y.crange, line = 1

        whereismin, chisqarr, wmn, /silent
        oplot, [spt(wmn)], [chisqarr(wmn)], ps = 6, sym = 2
        reg2dev, 0.1, 0.9, xd, yd
        xyouts, xd, yd, /dev, $
          'best SpT (fitted) = '+string(sptbest, '(F4.1)')+$
          textoidl('!Cbest SpT (min \chi2) = '+strc(spt(wmn(0)))),  $
          chars = 1.5
        datestamp, 'FITSPEC.PRO'

        ; over plot best fitting spectrum (based on chi-sq minima)
        if (!d.name eq 'X') then win, 1, 500 & lincolr, /silent
        w = (nearest(spt, sptbest))(0)
        readcol2, /silent, specdir(w)+speclist(w), wlcal0, fcal0, form = 'f,f'
        fcal = interpol(fcal0, wlcal0, wlobj)
        fcal = fcal/median(fcal)
        plot, wlobj, fobj, /xs,  $
          tit = 'best fitting spectrum (fitted chi-sq minima)', chars = 2
        oplot, wlobj, fcal, col = 2, thick = 1.5*!p.thick

    endelse
    message, 'min reduced chi-square = '+strc(min(chisqarr)), /info

    datestamp, 'FITSPEC.PRO'
    if (!d.name eq 'PS') then begin
        ps_close
        gv, outfile
    endif

endif


end