PRO fitplanck,w,bl,dbl,params,rmse,noplot=noplot,noprint=noprint
; INPUTS:
; w     = wavelengths of measured brightnesses (Angstroms)
; bl    = measured brightnesses (kR/A)
; dbl   = probable errors in measured brightnesses (kR/A)
; params = 2-element array:
; params(0) = t     = trial blackbody temperature
; params(1) = meanf = input user guess at best-fit factor
; OPTIONAL INPUT:
; noplot - if nonzero, turns off plotting
; noprint - if nonzero, inhibits PRINT statements
; OUTPUT:
; params (as above under INPUTS, but params(1) will have been modified
; OPTIONAL OUTPUTS
; rmse  = RMS error of fit (in units of probable errors in measurements
planck,w,params,pl & z=bl/pl & dz=dbl/pl	; Normalize observations to blackbody
gd=WHERE(dbl lt 5)		; Treat only well-defined brightnesses
; Get minimum, maximum, and average brightness ratios of PLBS to blackbody
minf=MIN(z(gd)+dz(gd)) & maxf=MAX(z(gd)-dz(gd)) & params(1)=(minf+maxf)/2.
IF NOT KEYWORD_SET(noprint) THEN $
    PRINT,params(1),FORMAT='("Best fit factor = ",E8.2)'
dif=bl-pl*params(1)		; Errors in blackbody fit
; Get RMS error of blackbody fit relative to measurement errors
rmse=SQRT(TOTAL((dif/dbl)^2)/N_ELEMENTS(dbl))	; get MEAN squared error
IF NOT KEYWORD_SET(noprint) THEN $
    PRINT,rmse,FORMAT='("RMS error = ",F6.3)'	; display MEAN squared error
; plot observed brightnesses and error bars
IF NOT KEYWORD_SET(noplot) THEN BEGIN
    PLOT,w,bl,PSYM=7,XTITLE='Wavelength ('+angstrom()+')', $
              YTITLE='Brightness (kR/'+angstrom()+')', $
              TITLE='Blackbody Fit to PLBS Brightness', $
              XRANGE=[5000,8500],XSTYLE=1
    FOR i=0,N_ELEMENTS(bl)-1 DO $
        PLOTS,[w(i),w(i)],[bl(i)+dbl(i),bl(i)-dbl(i)],NOCLIP=0
    wp=FINDGEN(21)*200.+5000.	; Create regular wavelength grid
    planck,w,params,plp		; Planck function on regular grid
    OPLOT,wp,plp*params(1)	; Overplot normalized blackbody fit
    xl=0.95*!X.CRANGE(0)+0.05*!X.CRANGE(1)
    XYOUTS,xl,0.9*!Y.CRANGE(1),STRING(FIX(t),FORMAT='(I4," K blackbody")')
    XYOUTS,xl,0.85*!Y.CRANGE(1),STRING(params(1),FORMAT='("(*",E8.2,")")')
    XYOUTS,xl,0.8*!Y.CRANGE(1),STRING(rmse,FORMAT='("RMS error = ",F5.2," !4r!x")')
    PRINT,'Enter identification label for this plot'
    id='' & READ,id
    XYOUTS,TOTAL(!X.CRANGE)/2.,0.95*!Y.CRANGE(1),id,ALIGNMENT=0.5
ENDIF
RETURN
END