PRO dustem_plot_fit_sed,st,sed, $
                         res=res,errors=errors,chi2=chi2,rchi2=rchi2, $
                         no_spec_error=no_spec_error,title=title,$;xtit=xtit,ytit=ytit,xr=xr,yr=yr $
                         help=help,win=win,ps=ps, $
                         legend_xpos=legend_xpos,legend_ypos=legend_ypos,legend_offset=legend_offset,fpol=fpol, $
                         _extra=_extra

; 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.
;-

;stop

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,set_character_size=[170,250], /encapsulated
ENDIF ELSE BEGIN

  set_plot,'X'
    IF keyword_set(win) then window,win;,xsize=600,ysize=800
ENDELSE

fact = 1.e4*(*!dustem_HCD)/(4.*!pi)/(3.e8/1.e-6/st.sed.wav)*1.e20/1.e7
;use_col_data_filt=70
;use_col_sed_spec=170
use_col_data_filt='blue'
;use_col_sed_spec='red'
use_col_sed_spec='grey'
IF not keyword_set(col_sed) THEN BEGIN
  ;use_col_sed_filt=250     ;red
  use_col_sed_filt='red'     ;red
ENDIF ELSE BEGIN
  use_col_sed_filt=col_sed
ENDELSE
IF not keyword_set(col_tot) THEN BEGIN
  ;use_col_tot=200
  use_col_tot='black'
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


spec = st.sed.em_tot * fact
if keyword_set(fpol) then specpol = st.polsed.em_tot * fact
;ADDING PLUGIN(S) TO SPECTRUM----------------
;if n_tags(!dustem_data.sed) gt 1 then begin 
scopes=tag_names((*!dustem_scope))
IF scopes[0] NE 'NONE' THEN BEGIN
;IF ptr_valid(!dustem_plugin) THEN BEGIN
  FOR i=0L,n_tags(*!dustem_scope)-1 DO BEGIN
    IF total(strsplit((*(*!dustem_scope).(i)),'+',/extract) EQ 'ADD_SED') THEN spec+=(*(*!dustem_plugin).(i))[*,0]
    if keyword_set(fpol) then begin
        IF total(strsplit((*(*!dustem_scope).(i)),'+',/extract) EQ 'REPLACE_QSED') THEN specpol=sqrt(((*(*!dustem_plugin).(i))[*,1])^2+((*(*!dustem_plugin).(i))[*,2])^2)
    endif
  ENDFOR
  FOR i=0L,n_tags(*!dustem_scope)-1 DO BEGIN
    if keyword_set(fpol) then begin    
        IF total(strsplit((*(*!dustem_scope).(i)),'+',/extract) EQ 'ADD_QSED') THEN specpol+=sqrt(((*(*!dustem_plugin).(i))[*,1])^2+((*(*!dustem_plugin).(i))[*,2])^2)
    endif  
  ENDFOR
ENDIF
;endif
;------------------------------------------

;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_cols=dustem_grains_colors(Ngrains,/cgplot)
use_lines=replicate(0,Ngrains)

norm = sed * 0. + 1

;====== PLOT THE SED

IF keyword_set(title) THEN title = title ELSE title = 'Spectral Energy Distribution (Running)'

IF keyword_set(xr) THEN xr = xr ELSE xr = [1.00E+00,6.00E+04]

IF keyword_set(yr) THEN yr = yr ELSE yr = [1.00E-7,5.00E03]

IF keyword_set(xtit) THEN xtit = xtit ELSE xtit = textoidl('\lambda (\mum)')

IF keyword_set(ytit) THEN ytit = ytit ELSE ytit = textoidl('Brightness/N_H (MJy/sr/H)')


;deffo define a title variable here. The procedures are not communicating with each other. 

;###############################
if keyword_set(fpol) then begin
    if !run_lin then begin
      IF !d.name NE 'PS' THEN cgDisplay, 600, 500
      cgplot,(*!dustem_data.sed).wav,(*!dustem_data.sed).values/norm,/nodata,xtit='',ytit=ytit,tit='',/ylog,/xlog,/ys,/xs,position=[0.12,0.25,0.96,0.76],xtickformat='(A1)',_extra=_extra,charsize=1.3,/noerase    
    endif
endif ELSE cgplot,(*!dustem_data.sed).wav,(*!dustem_data.sed).values/norm,/nodata,xtit='',ytit=ytit,tit=title,/ylog,/xlog,/ys,/xs,position=[0.12,0.35,0.96,0.90],xtickformat='(A1)',_extra=_extra,charsize=1.3
;###############################


;cgplot,(*!dustem_data.sed).wav,(*!dustem_data.sed).values/norm,/nodata,xtit='',ytit=textoidl('Brightness/N_H (MJy/sr/H)'),tit=title,ylog=1,/xlog,xr=xr,yr=yr,/ys,/xs,position=[0.12,0.35,0.96,0.90],xtickformat='(A1)'


;stop
;cgplot,st.polsed.wav,Q_sed*fact,xtit='',ytit=ytitstq,tit=titstq,ylog=ylog,/xlog,xr=xr,yr=yr,/ys,/xs,position=[0.17,0.35,0.95,0.95],xtickformat='(A1)'

;this is where you will add the normalized sed


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
; following lines for fawlty compatibility
defsysv,'!psym',exists=pexist
if pexist eq 0 then defsysv,'!psym',0
IF count_spec NE 0 THEN BEGIN
  plotsym,0,/fill
  xx=((*!dustem_data.sed).wav)[ind_spec]
  yy=((*!dustem_data.sed).values)[ind_spec]/norm[ind_spec]
  rms=3.*((*!dustem_data.sed).sigma)[ind_spec]/2./norm[ind_spec]
  cgoplot,xx,yy,psym=8,syms=0.5,color=use_col_sed_spec
  IF not keyword_set(no_spec_error) THEN BEGIN
    cgerrplot,xx,yy-rms,yy+rms,color=use_col_sed_spec
    ;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
  ;stop
  xx=((*!dustem_data.sed).wav)[ind_filt]
  yy=((*!dustem_data.sed).values)[ind_filt]/norm[ind_filt]
  rms=3.*((*!dustem_data.sed).sigma)[ind_filt]/2./norm[ind_filt]
  plotsym,0,/fill
  cgoplot,xx,yy,psym=8,color='Dodger Blue';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
  cgerrplot,xx,yy-rms,yy+rms,color='Dodger Blue';use_col_data_filt
ENDIF
;=== Plot the computed SED
IF count_filt NE 0 THEN BEGIN
  plotsym,8
  xx=((*!dustem_data.sed).wav)[ind_filt]
  yy=sed[ind_filt]/norm[ind_filt]
  cgoplot,xx,yy,color=use_col_sed_filt,psym=6,syms=2
ENDIF
IF count_spec NE 0 THEN BEGIN
  plotsym,0
  xx=((*!dustem_data.sed).wav)[ind_spec]
  yy=sed[ind_spec]/norm[ind_spec]
  cgoplot,xx,yy,color=use_col_sed_filt,psym=6,syms=2
ENDIF


IF !dustem_show_plot EQ 2 THEN BEGIN
	norm = spec
ENDIF ELSE BEGIN
	norm = spec * 0. + 1
ENDELSE
use_cols[1]='Cornflower'
FOR i=0L,Ngrains-1 DO BEGIN
  cgoplot,st.sed.wav,st.sed.(i+1)*fact/norm,color=use_cols[i],linestyle=use_lines[i]
ENDFOR

;stop

;PLOTTING OF THE PLUGIN(S)--------------- AUTOMATE THIS. QUITE FEASIBLE
IF tag_exist(*!dustem_scope,'CONTINUUM') THEN BEGIN
  cgoplot,st.sed.wav,(*(*!dustem_plugin).continuum)[*,0],color='Teal',linestyle=3
ENDIF
IF tag_exist(*!dustem_scope,'FREEFREE') THEN BEGIN
  cgoplot,st.sed.wav,(*(*!dustem_plugin).freefree)[*,0],color='Dark Red',linestyle=3
ENDIF
IF tag_exist(*!dustem_scope,'SYNCHROTRON') THEN BEGIN
  cgoplot,st.sed.wav,(*(*!dustem_plugin).synchrotron)[*,0],color='Crimson',linestyle=3
ENDIF

IF tag_exist(*!dustem_scope,'MBBDY_ISRF') THEN BEGIN
  cgoplot,st.sed.wav,(*(*!dustem_plugin).mbbdy_isrf)[*,0],color='Gold',linestyle=3
ENDIF

IF tag_exist(*!dustem_scope,'MBBDY') THEN BEGIN
  cgoplot,st.sed.wav,(*(*!dustem_plugin).mbbdy)[*,0],color='Cornflower',linestyle=3
ENDIF


cgoplot,st.sed.wav,spec/norm,color=use_col_tot,linestyle=use_line_tot
if keyword_set(fpol) then cgoplot,st.polsed.wav,specpol,color='Teal',linestyle=4
;plot the normalized data as well. 
;----------------------------------------

;==== print the legend
frmt0='(A36)'
frmt1='(1E10.2)'
frmt2='(F7.2)'
use_legend_xpos=0.50
if keyword_set(fpol) then use_legend_ypos=0.70 ELSE use_legend_ypos=0.84
use_legend_offset=0.03             ;This is the offset between lines of the legend in normalized units
legend_charsize=0.86

k=0.
iscond=0
n_plgns = n_tags(*!dustem_scope)
inn=n_plgns/2+1
IF n_plgns mod 2 ne 0 THEN inn=n_plgns/2+2

IF keyword_set(legend_xpos) THEN use_legend_xpos=legend_xpos
IF keyword_set(legend_ypos) THEN use_legend_ypos=legend_ypos
IF keyword_set(legend_offset) THEN use_legend_offset=legend_offset

IF !d.name NE 'PS' THEN cleanplot

tg=0
prv_str=''
IF keyword_set(res) THEN BEGIN
  Npar=n_elements(res)
  ;print,'=============='
  
  FOR i=0L,Npar-1 DO BEGIN
    parameter_description=(*(*!dustem_fit).param_descs)[i]
    parameter_type=dustem_parameter_description2type(parameter_description,string_name=string_name)
    str=string(string_name+' = ',format=frmt0)+string(res[i],format=frmt1)
    IF keyword_set(errors) THEN BEGIN
      str=str+textoidl(' \pm ')+string(errors(i),format=frmt1)
    ENDIF
   
    xxpos=use_legend_xpos*0.07
    yypos=use_legend_ypos*0.03
    xyouts,xxpos,yypos,'Model: '+!dustem_model,color=0,/normal,charsize=legend_charsize
   
    IF STRUPCASE(strmid(strtrim(parameter_description,2),0,6)) eq 'DUSTEM' THEN BEGIN
       
;     xxpos=use_legend_xpos*0.07
;     yypos=use_legend_ypos*0.03
;     xyouts,xxpos,yypos,'Model: '+!dustem_model,color=0,/normal,charsize=legend_charsize
     
    xxpos=use_legend_xpos*0.75
    yypos=use_legend_ypos*1.03
    xyouts,xxpos,yypos,'--Plugins--',color=0,/normal,charsize=legend_charsize
    
    ii = strsplit(string_name,'_',count=countx) & ii = ii(countx-1)-1 ; Locating the last underscore to automate the extraction of the plugin's keyword
    
    mm=where(tag_names(*!dustem_scope) eq strupcase(strmid(string_name,7,ii-7)),coun) ; Selecting a plugin through matching the string name of the plugin form the scope system variable with the one read from the parameter description vector
 
    ;tg+=1
    ;if strmid(string_name,7,ii-7) eq prv_str then tg-=1
    prmtg=(*(*!dustem_paramtag).(mm))
    indtg=(strmid(string_name,ii+1)) & indtg=strmid(indtg,0,/reverse_offset)   
    ;stop
    indtg=fix(indtg)
    prmtg=prmtg[indtg-1] 
    
    str=string(strmid(string_name,7,ii-7)+' ['+strmid(string_name,ii+1)+']: '+prmtg+' = ',format=frmt0)+string(res[i],format=frmt1)
   
    prv_str=strmid(string_name,7,ii-7)
    
    
    
     IF keyword_set(errors) THEN BEGIN
       str=str+textoidl(' \pm ')+string(errors(i),format=frmt1)
     ENDIF
    
    xxpos=use_legend_xpos*0.46
    
    IF n_plgns gt 3 THEN BEGIN ;adapt the display of the plugins if there are more than three (display two columns)
        IF i eq inn+k THEN BEGIN
           k+=inn     
           ;iscond=1
           inn+=1
           xxpos=use_legend_xpos*0.06
        ENDIF
    ENDIF
              
;     IF iscond THEN BEGIN
;     inn+=1
;     xxpos=use_legend_xpos*0.06
;     ENDIF

    yypos=use_legend_ypos-(i-k)*use_legend_offset
    xyouts,xxpos,yypos,str,color=0,/normal,charsize=legend_charsize
    endif else begin
    IF STRUPCASE(strmid(strtrim(parameter_description,2),0,24)) eq '(*!DUSTEM_PARAMS).GRAINS' and STRUPCASE(strmid(strtrim(parameter_description,2),3,1,/reverse_offset)) EQ 'O' then begin
     indpop=fix(STRUPCASE(strmid(strtrim(parameter_description,2),12,1,/reverse_offset)))
     ;oo+=1
     ;string_name=(((*!dustem_params).GRAINS).grain_type)[oo-1]
     string_name=(((*!dustem_params).GRAINS).grain_type)[indpop]
     ;stop
    endif
    k+=1
    xxpos=use_legend_xpos*1.5
    yypos=use_legend_ypos*1.03
    xyouts,xxpos,yypos,'--Parameters--',color=0,/normal,charsize=legend_charsize
    yypos=use_legend_ypos-(i)*use_legend_offset 
    xxpos=use_legend_xpos*1.14
    str=string(string_name+' = ',format=frmt0)+string(res[i],format=frmt1)
    xyouts,xxpos,yypos,str,color=0,/normal,charsize=legend_charsize    
    endelse
  ENDFOR
  ;stop
ENDIF
IF keyword_set(chi2) THEN BEGIN
  xxpos=1.29*use_legend_xpos
  yypos=1.13*0.84;use_legend_ypos
  xyouts,xxpos,yypos,string('chi2=',format=frmt0)+string(chi2,format=frmt2),color=0,/normal,charsize=legend_charsize
ENDIF
IF keyword_set(rchi2) THEN BEGIN
  xxpos=1.29*use_legend_xpos
  yypos=1.13*0.84-use_legend_offset;use_legend_ypos-use_legend_offset
  xyouts,xxpos,yypos,string('red. chi2=',format=frmt0)+string(rchi2,format=frmt2),color=0,/normal,charsize=legend_charsize
ENDIF

xtit=textoidl('\lambda (\mum)')
;cgplot,(*!dustem_data.sed).wav,(*!dustem_data.sed).values/sed,/nodata,xtit=xtit,ytit='Normalized',tit='',/xlog,xr=xr,/ys,/xs,yr=[0,2],ylog=0,position=[0.12,0.14,0.96,0.35],/noerase,yticks=2,ymino=2,xticklen=0.1
;stop

IF keyword_set(_extra) THEN BEGIN
  extra_kept={XR:[0.,0.]}
  extra_tags=tag_names(_extra)
  ind=where(extra_tags EQ 'XR',count)
  IF count NE 0 THEN extra_kept.XR=_extra.(ind[0]) ;ELSE extra_kept=0
ENDIF

;stop
;,position=[0.12,0.14,0.96,0.35] old position before your polfrac display modification
if keyword_set(fpol) then begin
cgplot,(*!dustem_data.sed).wav,(*!dustem_data.sed).values/sed,_extra=extra_kept,/nodata,xtit=xtit,ytit='Normalized',tit='',/xlog,/ys,/xs,yr=[0,2],ylog=0,position=[0.12,0.1,0.96,0.25],/noerase,yticks=2,ymino=2,xticklen=0.1,ytickformat='(F6.2)',charsize=1.3
endif else cgplot,(*!dustem_data.sed).wav,(*!dustem_data.sed).values/sed,_extra=extra_kept,/nodata,xtit=xtit,ytit='Normalized',tit='',/xlog,/ys,/xs,yr=[0,2],ylog=0,position=[0.12,0.14,0.96,0.35],/noerase,yticks=2,ymino=2,xticklen=0.1,ytickformat='(F6.2)',charsize=1.3
;plot the normalized data as well. 
IF count_spec NE 0 THEN BEGIN
  xx=((*!dustem_data.sed).wav)[ind_spec]
  yy=((*!dustem_data.sed).values/sed)[ind_spec]
  cgoplot,xx,yy,psym=16,symsize=1,thick=2,color=use_col_sed_spec
  trois_sigma=(3.*((*!dustem_data.sed).sigma)/2./sed)[ind_spec]
  cgerrplot,xx,yy-trois_sigma,yy+trois_sigma,color=use_col_sed_spec
ENDIF
IF count_filt NE 0 THEN BEGIN
  xx=((*!dustem_data.sed).wav)[ind_filt]
  yy=((*!dustem_data.sed).values/sed)[ind_filt]
  cgoplot,xx,yy,psym=16,symsize=1,thick=2,color='Dodger Blue'
  trois_sigma=(3.*((*!dustem_data.sed).sigma)/2./sed)[ind_filt]
  cgerrplot,xx,yy-trois_sigma,yy+trois_sigma,color='Dodger Blue'
ENDIF

;stop
;cgoplot,(*!dustem_data.sed).wav,(*!dustem_data.sed).values/sed,psym=16,symsize=1,thick=2,color='Dodger Blue'
cgoplot,10^!x.crange,replicate(1.,2),color='black',linestyle=0
;cgoplot,st.sed.wav,spec/spec,color='black'
;cgerrplot,((*!dustem_data.sed).wav),((*!dustem_data.sed).values)/sed-3.*((*!dustem_data.sed).sigma)/2./sed,((*!dustem_data.sed).values)/sed+3.*((*!dustem_data.sed).sigma)/2./sed,color='Dodger Blue'


if keyword_set(fpol) then begin
    plotsym,0,/fill
    cgplot, (*!dustem_data.polfrac).wav,(*!dustem_data.polfrac).values*100,/nodata,xtit='',tit=title,ytit='',/xlog,/ylog,/ys,/xs,position=[0.12,0.76,0.96,0.91],charsize=1.3,yr=[1.0E-1,30.0],/noerase,xticklen=0.1,xtickformat='(A1)',_extra=extra_kept;,xr=_extra.xr
    xxpos=use_legend_xpos*0.35
    yypos=use_legend_ypos*1.25
    xyouts,xxpos,yypos,'Polarization fraction (%)',color=cgcolor('purple'),charsize=legend_charsize,/normal
    yypos=use_legend_ypos*1.20
    stringg='Polarization SED '+'('+textoidl('P_{\nu}')+')'
    xyouts,xxpos,yypos,stringg,color=cgcolor('teal'),charsize=legend_charsize,/normal
    pilotfx=fpol(0);((*!dustem_data.polfrac).values)(0);fpol(0);
    cgoplot, (*!dustem_data.polfrac).wav,(*!dustem_data.polfrac).values*100,charsize=1.3,psym=8,syms=1,thick=2,color='Dodger Blue'   
    cgoplot, st.polsed.wav,(specpol/spec)*100/pilotfx*((*!dustem_data.polfrac).values)(0),color='purple'
    cgoplot, (*!dustem_data.polfrac).wav,fpol*100/pilotfx*((*!dustem_data.polfrac).values)(0),color='Red',psym=6,syms=1
    ;fpol=abs(fpol)
    ;delp=(*!dustem_data.polfrac).values*100-fpol*100/pilotfx*((*!dustem_data.polfrac).values)(0)
    ;cgoplot, ((*!dustem_data.polfrac).wav),((*!dustem_data.polfrac).values*100+delp),charsize=1.3,color='black',linestyle=2
    ;stop
endif

IF keyword_set(ps) THEN BEGIN
  device,/close
  set_plot,'X'
ENDIF

;stop

the_end:

END