PRO plotfit,w,bl,dbl,params,rmse,ymax=ymax
; 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) = factor = multiplier to adjust blackbody luminosity to PLBS
; OPTIONAL OUTPUTS
; rmse   = RMS error of fit (in units of probable errors in measurements)
IF N_Elements(ymax) EQ 0 THEN ymax=100  ; Safe upper bound
planck,w,params,pl 		; Calculate Planck function
dif=bl-pl			; 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
; plot observed brightnesses and error bars
PLOT,w,bl,PSYM=7 $
    ,XTITLE='Wavelength ('+angstrom()+')',XRANGE=[5000,8500],XSTYLE=1 $
    ,YTITLE='Brightness (kR/'+angstrom()+')',YRANGE=[0,MAX(bl+dbl)<ymax] $
    ,TITLE='Blackbody Fit to PLBS Brightness'
          
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,wp,params,plp		; Planck function on regular grid
OPLOT,wp,plp			; Overplot normalized blackbody fit
xl=0.95*!X.CRANGE(0)+0.05*!X.CRANGE(1)
XYOUTS,xl,0.925*!Y.CRANGE(1),STRING(params(0),FORMAT='(F5.0," K blackbody")')
XYOUTS,xl,0.875*!Y.CRANGE(1),STRING(params(1),FORMAT='("(*",E9.2,")")')
XYOUTS,xl,0.825*!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
timestmp
RETURN
END