dustem_plot_fit_sed.pro 7.93 KB
PRO dustem_plot_fit_sed,st,sed,cont,freefree,synchrotron, $
                        res=res,errors=errors,chi2=chi2,rchi2=rchi2, $
                        no_spec_error=no_spec_error, $
                        help=help,win=win,ps=ps, $
                        legend_xpos=legend_xpos,legend_ypos=legend_ypos, $
                        _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, /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
use_col_data_filt='blue'
use_col_sed_spec='red'
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
;color and linestyle for continuum
IF not keyword_set(col_cont) THEN BEGIN
  ;use_col_cont=30
  use_col_cont='Navy'
ENDIF ELSE BEGIN
  use_col_count=col_cont
ENDELSE
IF not keyword_set(line_cont) THEN BEGIN
  use_line_cont=0
ENDIF ELSE BEGIN
  use_line_cont=line_cont
ENDELSE
;color and linestyle for freefree
IF not keyword_set(col_freefree) THEN BEGIN
  ;use_col_cont=30
  use_col_freefree='Dark Red'
ENDIF ELSE BEGIN
  use_col_freefree=col_freefree
ENDELSE
IF not keyword_set(line_freefree) THEN BEGIN
  use_line_freefree=0
ENDIF ELSE BEGIN
  use_line_freefree=line_freefree
ENDELSE
IF not keyword_set(col_synchrotron) THEN BEGIN
  ;use_col_cont=30
  use_col_synchrotron='Dark Red'
ENDIF ELSE BEGIN
  use_col_synchrotron=col_synchrotron
ENDELSE
IF not keyword_set(line_synchrotron) THEN BEGIN
  use_line_synchrotron=0
ENDIF ELSE BEGIN
  use_line_synchrotron=line_synchrotron
ENDELSE

;spec = st.sed.em_tot * fact + cont
spec = st.sed.em_tot * fact
IF !dustem_idl_continuum THEN spec=spec+cont
IF !dustem_idl_freefree THEN spec=spec+freefree
IF !dustem_idl_synchrotron THEN spec=spec+synchrotron

;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

;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

;stop

;plot,(*!dustem_data.sed).wav,(*!dustem_data.sed).values/norm,/nodata,_extra=extra,yr=yr,/ys,ylog=ylog,/xlog
;stop
;===== waring: _extra not well taken by cgplot for unknown keywords ...
;cgplot,(*!dustem_data.sed).wav,(*!dustem_data.sed).values/norm,/nodata,_extra=extra,yr=yr,/ys,ylog=ylog,/xlog
cgplot,(*!dustem_data.sed).wav,(*!dustem_data.sed).values/norm,/nodata,yr=yr,/ys,ylog=ylog,/xlog,_extra=_extra
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
  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,_extra=extra,syms=0.5
  IF not keyword_set(no_spec_error) THEN BEGIN
    cgerrplot,xx,yy-rms,yy+rms
    ;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,_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
  cgerrplot,xx,yy-rms,yy+rms
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=8,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_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
  cgoplot,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
  cgoplot,st.sed.wav,cont,color=use_col_cont,linestyle=use_line_cont
ENDIF
IF !dustem_idl_freefree THEN BEGIN
  cgoplot,st.sed.wav,freefree,color=use_col_freefree,linestyle=use_line_freefree
ENDIF
IF !dustem_idl_synchrotron THEN BEGIN
  cgoplot,st.sed.wav,synchrotron,color=use_col_synchrotron,linestyle=use_line_synchrotron
ENDIF
cgoplot,st.sed.wav,spec/norm,color=use_col_tot,linestyle=use_line_tot

;==== print the legend
frmt0='(A20)'
frmt1='(1E10.2)'
frmt2='(F7.2)'
;WE have to find a way to pass this through keywords ...
;use_legend_xpos=30.
;use_legend_ypos=-1
use_legend_xpos=2.e3
use_legend_ypos=2.5
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(res) THEN BEGIN
  offset=0.2
  ;xy_xpos=30.
  ;ystart=-1
  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)
    ;stop
    IF keyword_set(errors) THEN BEGIN
      str=str+textoidl('\pm')+string(errors(i),format=frmt1)
    ENDIF
    xyouts,use_legend_xpos,10^(use_legend_ypos-(i+1)*offset),str,color=0
  ENDFOR
ENDIF
IF keyword_set(chi2) THEN BEGIN
  xyouts,use_legend_xpos,10^(use_legend_ypos-(Npar+1)*offset),string('chi2=',format=frmt0)+string(chi2,format=frmt2),color=0
ENDIF
IF keyword_set(rchi2) THEN BEGIN
  xyouts,use_legend_xpos,10^(use_legend_ypos-(Npar+2)*offset),string('red. chi2=',format=frmt0)+string(rchi2,format=frmt2),color=0
ENDIF

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

;stop

the_end:

END