dustem_plot_fit_sed.pro 15.9 KB
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239 240 241 242 243 244 245 246 247 248 249 250 251 252 253 254 255 256 257 258 259 260 261 262 263 264 265 266 267 268 269 270 271 272 273 274 275 276 277 278 279 280 281 282 283 284 285 286 287 288 289 290 291 292 293 294 295 296 297 298 299 300 301 302 303 304 305 306 307 308 309 310 311 312 313 314 315 316 317 318 319 320 321 322 323 324 325 326 327 328 329 330 331 332 333 334 335 336 337 338 339 340 341 342 343 344 345 346 347 348 349 350 351 352 353 354 355 356 357 358 359 360 361 362 363 364 365 366 367 368 369 370 371 372 373 374 375 376 377 378 379 380 381 382 383 384 385 386 387 388 389 390 391 392 393 394 395 396 397 398 399 400 401 402 403 404 405 406 407 408 409 410 411 412 413 414 415 416 417 418 419 420 421 422 423 424 425 426 427
 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.
;-


;NOta bene: only delete the cgoplot commands

;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='Powder Blue'

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

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)')


;###############################
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,ytit=ytit,tit=title,/ylog,/xlog,/ys,/xs,position=[0.12,0.35,0.96,0.90],charsize=1.3,_extra=_extra,xtickformat='(A1)'
;stop  
;###############################

;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,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='Indian Red',psym=7,syms=1
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