PRO dustem_plot_fit_sed,st,sed,cont, $
                        res=res,errors=errors,chi2=chi2,rchi2=rchi2, $
                        no_spec_error=no_spec_error, $
                        help=help,_extra=extra,win=win,ps=ps

;+
; NAME:
;    dustem_plot_fit_sed
; PURPOSE:
;    Plots a Dustem model and SED. Parameter values and error are
;    printed on plot. Used for plotting results during fit.
; CATEGORY:
;    Dustem
; CALLING SEQUENCE:
;    dustem_plot_fit_sed,st,sed,cont[,/no_spec_error][,res=][,errors=][,chi2=][,rchi2=][/help][_extra=]
; INPUTS:
;    st        = Dustem model output structure
;    sed       = Dustem model SED
;    cont      = Dustem model NIR continuum
; OPTIONAL INPUT PARAMETERS:
;    res       = fit result values. If set values are written on plot.
;    errors    = fit result errors. If set values are written on plot.
;    chi2      = fit chi^2. if set value is written on plot.
;    rchi2     = Reduced fit chi^2. if set value is written on plot.
;    _extra    = extra parameters for the plot routine
; OUTPUTS:
;    None
; OPTIONAL OUTPUT PARAMETERS:
;    None
; ACCEPTED KEY-WORDS:
;    help      = If set, print this help
; COMMON BLOCKS:
;    None
; SIDE EFFECTS:
;    SED and model are plotted
; RESTRICTIONS:
;    The dustem idl wrapper must be installed
; PROCEDURE:
;    None
; EXAMPLES
;    dustem_plot_fit_sed,st,sed,cont
; MODIFICATION HISTORY:
;    Written by J.-Ph. Bernard
;    see evolution details on the dustem cvs maintained at CESR
;    Contact J.-Ph. Bernard (Jean-Philippe.Bernard@cesr.fr) in case of problems.
;-

IF keyword_set(help) THEN BEGIN
  doc_library,'dustem_plot_fit_sed'
  goto,the_end
ENDIF

if keyword_set(ps) then begin
        set_plot, 'PS'
        device, filename=ps, /color, /encapsulated
endif else begin
        set_plot,'X'
        if keyword_set(win) then window,win
endelse

fact = 1.e4*1.E20/(4.*!pi)/(3.e8/1.e-6/st.sed.wav)*1.e20/1.e7

use_col_data_filt=70
use_col_sed_spec=170
IF not keyword_set(col_sed) THEN BEGIN
  use_col_sed_filt=250     ;red
ENDIF ELSE BEGIN
  use_col_sed_filt=col_sed
ENDELSE
IF not keyword_set(col_tot) THEN BEGIN
  use_col_tot=200
ENDIF ELSE BEGIN
  use_col_tot=col_tot
ENDELSE
IF not keyword_set(line_tot) THEN BEGIN
  use_line_tot=0
ENDIF ELSE BEGIN
  use_line_tot=line_tot
ENDELSE
IF not keyword_set(col_pah) THEN BEGIN
  use_col_pah=50    ;blue
ENDIF ELSE BEGIN
  use_col_pah=col_pah
ENDELSE
IF not keyword_set(line_pah) THEN BEGIN
  use_line_pah=2      ;dash, like in Desert et al.
ENDIF ELSE BEGIN
  use_line_pah=line_pah
ENDELSE
IF not keyword_set(col_vsg) THEN BEGIN
  use_col_vsg=150    ;green
ENDIF ELSE BEGIN
  use_col_vsg=col_vsg
ENDELSE
IF not keyword_set(line_vsg) THEN BEGIN
  use_line_vsg=1      ;dots, like in Desert et al.
ENDIF ELSE BEGIN
  use_line_vsg=line_vsg
ENDELSE
IF not keyword_set(col_bg) THEN BEGIN
  use_col_bg=250    ;red
ENDIF ELSE BEGIN
  use_col_bg=col_bg
ENDELSE
IF not keyword_set(line_bg) THEN BEGIN
  use_line_bg=3      ;dot-dash, like in Desert et al.
ENDIF ELSE BEGIN
  use_line_bg=line_bg
ENDELSE
IF not keyword_set(col_cont) THEN BEGIN
  use_col_cont=30
ENDIF ELSE BEGIN
  use_col_count=col_cont
ENDELSE
IF not keyword_set(line_cont) THEN BEGIN
  use_line_cont=5    ;long dash
ENDIF ELSE BEGIN
  use_line_cont=line_cont
ENDELSE


;spec = st.sed.em_tot * fact + cont
spec = st.sed.em_tot * fact
IF !dustem_idl_continuum THEN spec=spec+cont

;stop

;use_cols=[use_col_pah,use_col_vsg,use_col_bg,use_col_cont]
;use_lines=[use_line_pah,use_line_vsg,use_line_bg,use_line_cont]
col_off=30
;Ngrains=(*!dustem_params).grain.Ngrains
Ngrains=(*!dustem_params).Ngrains
use_cols=long(findgen(Ngrains)/(Ngrains-1)*(255-col_off)+col_off)
use_lines=replicate(1,Ngrains)

if !dustem_show_plot EQ 2 then begin
	norm = sed
	yr=[0.,2.]
	ylog = 0
endif else begin
	norm = sed * 0. + 1
	yr=[1e-4,10]
	ylog = 1
endelse
plot,(*!dustem_data.sed).wav,(*!dustem_data.sed).values/norm,/nodata,_extra=extra,yr=yr,/ys,ylog=ylog,/xlog
ind_filt=where((*!dustem_data.sed).filt_names NE 'SPECTRUM',count_filt)
ind_spec=where((*!dustem_data.sed).filt_names EQ 'SPECTRUM',count_spec)
;=== Plot the data
IF count_spec NE 0 THEN BEGIN
  plotsym,0,/fill
  oplot,((*!dustem_data.sed).wav)(ind_spec),((*!dustem_data.sed).values)(ind_spec)/norm(ind_spec),psym=8,_extra=extra,syms=0.5
  IF not keyword_set(no_spec_error) THEN BEGIN
    err_bar,((*!dustem_data.sed).wav)(ind_spec),((*!dustem_data.sed).values)(ind_spec)/norm(ind_spec),yrms=3.*((*!dustem_data.sed).sigma)(ind_spec)/2./norm(ind_spec)
  ENDIF
ENDIF
IF count_filt NE 0 THEN BEGIN
  plotsym,0,/fill
  oplot,((*!dustem_data.sed).wav)(ind_filt),((*!dustem_data.sed).values)(ind_filt)/norm(ind_filt),psym=8,_extra=extra,color=use_col_data_filt
  err_bar,((*!dustem_data.sed).wav)(ind_filt),((*!dustem_data.sed).values)(ind_filt)/norm(ind_filt),yrms=3.*((*!dustem_data.sed).sigma)(ind_filt)/2./norm(ind_filt),color=use_col_data_filt
ENDIF
;=== Plot the computed SED
IF count_filt NE 0 THEN BEGIN
  plotsym,8
  oplot,((*!dustem_data.sed).wav)(ind_filt),sed(ind_filt)/norm(ind_filt),color=use_col_sed_filt,psym=8,syms=2
ENDIF
IF count_spec NE 0 THEN BEGIN
  plotsym,0
  oplot,((*!dustem_data.sed).wav)(ind_spec),sed(ind_spec)/norm(ind_spec),color=use_col_sed_spec,psym=8,syms=0.5
ENDIF


if !dustem_show_plot EQ 2 then begin
	norm = spec
endif else begin
	norm = spec * 0. + 1
endelse
FOR i=0L,Ngrains-1 DO BEGIN
  oplot,st.sed.wav,st.sed.(i+1)*fact/norm,color=use_cols(i),linestyle=use_lines(i)
ENDFOR
;NIR continuum
IF !dustem_idl_continuum THEN BEGIN
  oplot,st.sed.wav,cont,color=use_col_cont,linestyle=use_line_cont
ENDIF
;;PAH emission
;oplot,st.sed.wav,st.sed.em_pah*fact,color=use_col_pah,linestyle=use_line_pah
;;VSG emission
;oplot,st.sed.wav,st.sed.em_vsg*fact,color=use_col_vsg,linestyle=use_line_vsg
;;BG emission
;oplot,st.sed.wav,st.sed.em_bg*fact/norm,color=use_col_bg,linestyle=use_line_bg
;Total emission
;oplot,st.sed.wav,(st.sed.(4)+st.sed.(5)+st.sed.(1)+st.sed.(2))*fact/norm,color=use_col_tot,linestyle=use_line_bg
oplot,st.sed.wav,spec/norm,color=use_col_tot,linestyle=use_line_tot

;if !dustem_show_plot EQ 2 then return

frmt0='(A20)'
frmt1='(1E10.2)'
frmt2='(F7.2)'
IF keyword_set(res) THEN BEGIN
  offset=0.2
  xy_xpos=2.
  Npar=n_elements(res)
  FOR i=0L,Npar-1 DO BEGIN
;    sstr=(str_sep((*!PDO)(i),'(*!dustem_params).'))
    sstr=(str_sep((*(*!dustem_fit).param_descs)(i),'(*!dustem_params).'))
    ssstr=sstr(n_elements(sstr)-1)
    str=string(ssstr+'=',format=frmt0)+string(res(i),format=frmt1)
    IF keyword_set(errors) THEN BEGIN
      str=str+textoidl('\pm')+string(errors(i),format=frmt1)
    ENDIF
    xyouts,xy_xpos,10^(1.-(i+1)*offset),str
  ENDFOR
ENDIF
IF keyword_set(chi2) THEN BEGIN
  xyouts,xy_xpos,10^(1.-(Npar+1)*offset),string('chi2=',format=frmt0)+string(chi2,format=frmt2)
ENDIF
IF keyword_set(rchi2) THEN BEGIN
  xyouts,xy_xpos,10^(1.-(Npar+2)*offset),string('red. chi2=',format=frmt0)+string(rchi2,format=frmt2)
ENDIF

if keyword_set(ps) then begin
        device,/close
        set_plot,'X'
endif

the_end:

END