PRO dustem_plot_polar,st,p_dim,ps=ps,_Extra=_extra,help=help,UV=UV,SED=SED,Pfrac=Pfrac,print_ratio=print_ratio,win=win,xr=xr,yr=yr,donotclose=donotclose,aligned=aligned,noerrbars=noerrbars,multi=multi,almabands=almabands,nsmooth=nsmooth ;+ ; NAME: ; dustem_plot_polar ; PURPOSE: ; Plots a Dustem polarization curve ; CATEGORY: ; Dustem ; CALLING SEQUENCE: ; dustem_plot_curve,st,,/help][,_extra=] ; INPUTS: ; st = Dustem model output structure ; OPTIONAL INPUT PARAMETERS: ; _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: ; model is plotted ; RESTRICTIONS: ; The dustem idl wrapper must be installed ; PROCEDURE: ; None ; EXAMPLES ; dustem_init,/wrap ; dir_in=getenv('DUSTEM_SOFT_DIR')+'/src/dustem3.8_web/' ; dir=getenv('DUSTEM_RES') ; st=dustem_read_all_res(dir_in) ; loadct,13 ; dustem_plot_curve,st,yr=[1.e-6,1.e0],/ysty,xr=[1.,400],/xsty ; MODIFICATION HISTORY: ; ; Adapted by VG from dustem_plot_ext.pro ; ; 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_polar' goto,the_end ENDIF if not keyword_set(nsmooth) then nsmooth = 1 nsmooth = 4 almabands = 1. if not keyword_set(multi) then multi=0 if multi ne 2 then begin 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 begin if ptr_valid(!dustem_data.polsed) then window,win,title='DUSTEM WRAP (POLSED)' else window,win,title='DUSTEM WRAP (POLEXT)' endif endelse endif ;stop use_col_sed_filt = 70 use_col_data_filt = 250 Ngrains=n_tags(st.sed)-2 colors=dustem_grains_colors(Ngrains,ls,ps=ps) if ptr_valid(!dustem_data.polsed) then xtit=textoidl('\lambda (\mum)') else xtit=textoidl('\lambda^{-1} (\mum^{-1})') ;ext_tot=st.ext.abs_grain(0)*0. ;polext_tot=st.polext.abs_grain(0)*0. ;FOR i=0L,Ngrains-1 DO BEGIN ;polext_tot = polext_tot + (st.polext.abs_grain(i) + st.polext.sca_grain(i)) ;ext_tot = ext_tot + (st.ext.abs_grain(i) + st.ext.sca_grain(i)) ;ENDFOR ; Are observational data defined ? defsysv,'!dustem_data',exists=data_exists IF ptr_valid(!dustem_data.polext) THEN BEGIN interpolext = interpol(st.polext.ext_tot,st.polext.wav,(*!dustem_data.polext).wav) IF keyword_set(print_ratio) THEN BEGIN ; Correction de couleur cc = 1.11 Bband = 0.44 Vband = 0.55 Rband = 0.66 Iband = 0.8 Jband = 1.25 Hband = 1.6 Kband = 2.2 ; pV and tauV are In units of 1e-21 cm2/H pV = interpol(st.polext.ext_tot,st.polext.wav,Vband) pR = interpol(st.polext.ext_tot,st.polext.wav,Rband) pI = interpol(st.polext.ext_tot,st.polext.wav,Iband) pJ = interpol(st.polext.ext_tot,st.polext.wav,Jband) pH = interpol(st.polext.ext_tot,st.polext.wav,Hband) pK = interpol(st.polext.ext_tot,st.polext.wav,Kband) j=where(st.polext.wav gt 5 and st.polext.wav lt 20) pSil=max(st.polext(j).ext_tot,jmax) Silband = st.polext(j(jmax)).wav j=where(st.polext.wav gt 0.4 and st.polext.wav lt 0.75) pmeanoptical = mean(st.polext(j).ext_tot) pV_data = interpol((*!dustem_data.polext).values,(*!dustem_data.polext).wav,Vband) tauV = interpol(st.ext.ext_tot,st.ext.wav,Vband) tauV_data = interpol((*!dustem_data.ext).values,(*!dustem_data.ext).wav,Vband) tauB = interpol(st.ext.ext_tot,st.ext.wav,Bband) tauB_data = interpol((*!dustem_data.ext).values,(*!dustem_data.ext).wav,Bband) print ; I353 and P353 are in fact nu_Pnu, in unit of erg/s/H, not erg/s/sr/H ; 1 W/m2/sr/1d20H = 1e-17 erg/s/sr/H ENDIF IF ptr_valid(!dustem_data.polsed) THEN BEGIN dustem_polsed = dustem_compute_polsed(p_dim,st);dustem_compute_polsed(p_dim,st.polsed,cont=cont) x=min(((*!dustem_data.polsed).wav-850)^2,jmin) P353 = dustem_polsed(jmin) P353_data = (*!dustem_data.polsed).values(jmin) dustem_sed = dustem_compute_sed(p_dim,st);dustem_compute_sed(p_dim,st,cont=cont) x=min(((*!dustem_data.sed).wav-850)^2,jmin) I353 = dustem_sed(jmin) I353_data = (*!dustem_data.sed).values(jmin) ; Data is Inu, not nInu, in MJy/sr/1e20 ; Convert data, which is nu_Inu SED in MJy/sr/1d20H into nu_Inu SED model in erg/s/H fact = (4*!pi) * 1e-17 / 1e20 * 353e9 print print,"Observational constraint" print,"-------------------------" print,"I353/NH data ",I353_data," MJy/sr for NH=1e20" print,"P353/NH data ",P353_data," MJy/sr for NH=1e20" print,"P353/I353 data ",P353_data/I353_data*100," %" print print,"Model results" print,"-------------" print,"I353/NH model ",I353," MJy/sr for NH=1e20" print,"P353/NH model ",P353," MJy/sr for NH=1e20" print,"P353/I353 model ",P353/I353*100," %" IF keyword_set(print_ratio) THEN BEGIN print print,"Observational constraint" print,"-------------------------" print,"tau_V data ",tauV_data print,"(p/tau)_V data ",pV_data/tauV_data*100," %" print,"p_V/A_V data ",pV_data/(tauV_data*1.086)*100," %" print,"with p_V data ",pV_data," %" print,"Rs/v data ",(P353_data/I353_data)/(pv_data/tauv_data) print,"RP/p data ",(P353_data/1e20) / (pv_data*1e-21)," MJy/sr" print,"I353/Av data ",(I353_data/1e20) / (1.086*tauv_data*1e-21)," MJy/sr" print print,"(p/tau)_V model ",pV/tauV*100," %" print,"p_V/A_V model ",pV/(tauV*1.086)*100," %" print,"p_V/E(B-V) model ",pV/((tauB-tauV)*1.086)*100," %" print,"with p_V model ",pV," %" print,"Rs/v model ",P353/I353/(pV/tauV) print,"Rs/EBV model ",P353/I353/(pV/(tauB-tauV)) print,"RP/p model ",(P353/1e20) / (pv*1e-21)," MJy/sr" print,"- R band = 0.66 um ",(P353/1e20) / (pR*1e-21)," MJy/sr" print,"- I band = 0.8 um ",(P353/1e20) / (pI*1e-21)," MJy/sr" print,"- J band = 1.25 um ",(P353/1e20) / (pJ*1e-21)," MJy/sr" print,"- H band = 1.6 um ",(P353/1e20) / (pH*1e-21)," MJy/sr" print,"- K band = 2.2 um ",(P353/1e20) / (pK*1e-21)," MJy/sr" print,"- "+string(Silband,format='(F5.2)')+" um band ",(P353/1e20) / (pSil*1e-21)," MJy/sr" print,"- mean optical ",(P353/1e20) / (pmeanoptical*1e-21)," MJy/sr" print,"I353/Av model ",(I353/1e20) / (1.086*tauv*1e-21)," MJy/sr" print ;print,"Mixt results" ;print,"-------------" ;print,"(pV_model/tauV_data) ",pV/tauV_data*100," %" ;print,"I353_model/AV_data ",(I353/1e20) / (1.086*tauv_data*1e-21)," MJy/sr" ;print,"P353model/Idata ",P353/I353_data*100," %" ;print,"Rs/v P353model/I353data/(pVmodel/tauVdata) ",P353/I353_data/(pV/tauV_data) ;print,"RP/p P353model/pVdata ",(P353/1e20)/(pV_data*1e-21)," MJy/sr" ;print ENDIF ENDIF IF keyword_set(UV) THEN BEGIN if not keyword_set(ps) then tit='Polarization Cross-section' ytit=textoidl('\sigma_{pol}/N_H (cm^2/H)') if ptr_valid(!dustem_data.polext) then begin if !dustem_show_plot EQ 2 or multi eq 2 then begin ;normpol = polext_tot normpol = st.polext.ext_tot dustem_polext = interpol(st.polext.ext_tot,st.polext.wav,(*!dustem_data.polext).wav) yr=[0.,2.] ylog = 0 endif else begin dustem_polext = (*!dustem_data.polext).wav * 0. + 1.d20 normpol = st.polext.ext_tot * 0. + 1.d20 ylog = 1 endelse if multi eq 0 then cgplot,st.polext.wav,st.polext.ext_tot/normpol,/nodata,xtit=xtit,ytit=ytit,tit=tit,xr=xr,yr=yr,ylog=ylog,xlog=1 if multi eq 1 then cgplot,st.polext.wav,st.polext.ext_tot/normpol,/nodata,_Extra=_extra,xtit='',ytit=ytit,tit=tit,xr=xr,yr=yr,ylog=ylog,position=[0.17,0.35,0.93,0.93],xtickformat='(A1)' ,xlog=1 if multi eq 2 then cgplot,st.polext.wav,st.polext.ext_tot/normpol,/nodata,_Extra=_extra,xtit=xtit,ytit='Normalized',tit='',xr=xr,yr=[0,2],ylog=ylog,position=[0.17,0.14,0.93,0.35],/noerase,yticks=3,ymino=5,xticklen=0.1,xlog=1 FOR i=0L,Ngrains-1 DO BEGIN if aligned(i) then cgoplot,1/st.polext.wav,(st.polext.abs_grain(i)+st.polext.sca_grain(i))/normpol,color=colors(i+1),psym=psym,line=ls(i+1) ENDFOR if total(aligned gt 0) then cgoplot,1/st.polext.wav,st.polext.ext_tot/normpol,psym=psym if data_exists then BEGIN if ptr_valid(!dustem_data.polext) then begin plotsym,0,/FILL cgoplot,1/(*!dustem_data.polext).wav,(*!dustem_data.polext).values/dustem_polext,psym=16,color='Olive Drab' if multi ne 2 then cgoplot,1/(*!dustem_data.polext).wav,interpolext/dustem_polext,psym=6,syms=2,color='red' if not keyword_set(noerrbars) then err_bar,(*!dustem_data.polext).wav,(*!dustem_data.polext).values/dustem_polext,yrms=(*!dustem_data.polext).sigma/dustem_polext,col=cgColor('Olive Drab') endif ENDIF ENDIF ENDIF if multi eq 0 then begin if not keyword_set(ps) then tit='Polarization fraction in extinction' ytit='Polarization fraction' cgplot,st.polext.wav,st.polext.ext_tot/st.ext.ext_tot,/nodata,_Extra=_extra,xtit=xtit,ytit=ytit,tit=tit,xr=xr,yr=yr FOR i=0L,Ngrains-1 DO BEGIN cgoplot,st.polext.wav,smooth((st.polext.abs_grain(i)+st.polext.sca_grain(i))/(st.ext.abs_grain(i)+st.ext.sca_grain(i)),nsmooth),color=colors(i+1),psym=psym,line=ls(i+1) ENDFOR if total(aligned gt 0) then cgoplot,st.polext.wav,st.polext.ext_tot/st.ext.ext_tot,psym=psym cgoplot,st.polsed.wav,st.polsed.(Ngrains+1)/st.sed.(Ngrains+1),psym=psym,linestyle=2 if data_exists and ptr_valid(!dustem_data.polext) and ptr_valid(!dustem_data.ext) then begin values_ext_for_polext = interpol((*!dustem_data.ext).values,(*!dustem_data.ext).wav,(*!dustem_data.polext).wav) sigma_ext_for_polext = interpol((*!dustem_data.ext).sigma ,(*!dustem_data.ext).wav,(*!dustem_data.polext).wav) if not keyword_set(noerrbars) then err_bar,(*!dustem_data.polext).wav,(*!dustem_data.polext).values/values_ext_for_polext, $ yrms=1./values_ext_for_polext * (*!dustem_data.polext).sigma + (*!dustem_data.polext).values/values_ext_for_polext^2 * sigma_ext_for_polext endif endif ENDIF IF ptr_valid(!dustem_data.polsed) THEN BEGIN if not keyword_set(ps) then tit='Polarization Intensity';' in emission' ytit='Polarization Intensity' ;is this line used? ytit=textoidl('\nuP_\nu/N_H (W/m^2/sr/H)') ; Convert from erg/s/H into W/m2/sr/H fact=1./(4.*!pi)*1.e17/1d20 ;dustem_polsed = dustem_compute_polsed(p_dim,st.polsed,cont=cont)*fact if not keyword_set(xr) then xr=[5e2,5e3] if !dustem_show_plot EQ 2 or multi eq 2 then begin ;this is when plotting relative to SED value normpol = st.polsed.(Ngrains+1)*fact yr=[0.,2.5] ylog = 0 endif else begin normpol = st.polsed.(Ngrains+1)*0. + 1. dustem_polsed = ((*!dustem_data.polsed).values)*0. + 1 ylog=1 if not keyword_set(yr) then yr=[1e-34,1e-30] endelse if multi eq 0 then cgplot,st.polsed.wav,/nodata,st.polsed.(Ngrains+1)*fact/normpol,xtit=xtit,ytit=ytit,tit=tit,ylog=ylog,/xlog,xr=xr,yr=yr,/ys,/xs if multi eq 1 then cgplot,st.polsed.wav,/nodata,st.polsed.(Ngrains+1)*fact/normpol,xtit='',ytit=ytit,tit=tit,ylog=ylog,/xlog,xr=xr,yr=yr,/ys,/xs,position=[0.17,0.35,0.93,0.93],xtickformat='(A1)' if multi eq 2 then cgplot,st.polsed.wav,/nodata,st.polsed.(Ngrains+1)*fact/normpol,xtit=xtit,ytit='Normalized',tit='',/xlog,xr=xr,/ys,/xs,yr=[0,2],ylog=0,position=[0.17,0.14,0.93,0.35],/noerase,yticks=2,ymino=5,xticklen=0.1 for i = 0, Ngrains-1 do begin if aligned(i) then cgoplot,st.polsed.wav,st.polsed.(i+1)*fact/normpol,color=colors(i+1),psym=psym,line=ls(i+1) ;investigate this line of the plotting of the different components in polarization endfor if total(aligned gt 0) && multi ne 2 then cgoplot,st.polsed.wav,st.polsed.(Ngrains+1)*fact/normpol,color='Teal' if total(aligned gt 0) && multi eq 2 then cgoplot,st.polsed.wav,st.polsed.(Ngrains+1)*fact/normpol,color='black' if data_exists then BEGIN if ptr_valid(!dustem_data.polsed) then begin ind_filt=where((*!dustem_data.polsed).filt_names NE 'SPECTRUM',count_filt) ;=== Plot the data IF count_filt NE 0 THEN BEGIN if !dustem_show_plot EQ 2 or multi eq 2 then begin endif else begin ; Convert data which Pnu in MJy/sr/1d20 into nu*Pnu in W/m2/sr/1d20 fact = (3.e8/(1.e-6*(*!dustem_data.polsed).wav(ind_filt)))/1.d40 endelse dustem_polsed = dustem_compute_polsed(p_dim,st)*fact;dustem_compute_polsed(p_dim,st,cont=cont)*fact if !dustem_show_plot EQ 2 or multi eq 2 then begin normpol = dustem_polsed endif else begin normpol = dustem_polsed * 0. + 1 endelse plotsym,0,/FILL cgoplot,((*!dustem_data.polsed).wav)(ind_filt),((*!dustem_data.polsed).values)(ind_filt)*fact/normpol(ind_filt),psym=8,_extra=_extra,color='Orchid',thick=2;use_col_sed_filt if multi ne 2 then cgoplot,((*!dustem_data.polsed).wav),dustem_polsed/normpol,psym=6,color='red',symsize=2 if not keyword_set(noerrbars) then err_bar,((*!dustem_data.polsed).wav)(ind_filt),((*!dustem_data.polsed).values)(ind_filt)*fact/normpol(ind_filt),yrms=3.*((*!dustem_data.polsed).sigma)(ind_filt)/2.*fact/normpol(ind_filt),color=cgColor('Orchid');use_col_sed_filt ; PLot color corrected POLSED if keyword_set(ps) then plotsym,8,thick=8 else plotsym,8 if multi eq 0 then cgoplot,((*!dustem_data.polsed).wav)(ind_filt),dustem_polsed(ind_filt)/normpol(ind_filt),psym=8,syms=2,_extra=_extra,col='red' ENDIF IF tag_exist(*!dustem_scope,'SYNCHROTRON') THEN BEGIN ;bad test but let's leave it like this for the moment ;indx=where((*(*!dustem_fit).param_descs) EQ 'dustem_plugin_synchrotron_3',counttx) ;if total(counttx) then frac_synch=p_dim(indx) else frac_synch=0.3 ;THIS LINE IS TO TAKE INTO ACCOUNT WHEN THE AMPLITUDE OF THE SYNCHROTRON IS FROZEN synchrotron_p=sqrt((*(*!dustem_plugin).synchrotron)[*,1]^2+(*(*!dustem_plugin).synchrotron)[*,2]^2) cgoplot,st.polsed.wav,synchrotron_p*(3.e8/(1.e-6*(st.polsed.wav)))/1.d40,color='Crimson',linestyle=3 cgoplot,st.polsed.wav,st.polsed.em_tot*1./(4.*!pi)*1.e17/1d20+synchrotron_p*(3.e8/(1.e-6*(st.polsed.wav)))/1.d40,color='black' ENDIF endif endif endif ;else BEGIN if keyword_set(Pfrac) then begin ; DEPRECATED ; Also show prediction for polarized SED if not keyword_set(ps) then tit='Polarization fraction in emission' ytit=textoidl('P/I') plot_oi,st.polsed.wav,st.polsed.(Ngrains+1)/st.sed.(Ngrains+1),/nodata,_extra=_extra,xtit=xtit,ytit=ytit,tit=tit,/xlog,xr=xr,yr=yr if (total(aligned) gt 0) then oplot,st.polsed.wav,st.polsed.(Ngrains+1)/st.sed.(Ngrains+1) for i = 0, Ngrains-1 do if aligned(i) then oplot,st.polsed.wav,st.polsed.(i+1)/st.sed.(Ngrains+1),color=colors(i+1),psym=psym,linestyle=ls(i+1) if keyword_Set(almabands) then begin oplot,870*[1,1],yr,linestyle=1 oplot,1300*[1,1],yr,linestyle=1 oplot,3100*[1,1],yr,linestyle=1 endif if data_exists then begin if ptr_valid(!dustem_data.polfrac) then begin ;=== Plot the data ind_filt_polfrac=where((*!dustem_data.polfrac).filt_names NE 'SPECTRUM',count_filt_polfrac) IF count_filt_polfrac NE 0 THEN BEGIN plotsym,0,/FILL oplot,((*!dustem_data.polfrac).wav)(ind_filt_polfrac),((*!dustem_data.polfrac).values)(ind_filt_polfrac),psym=8,_extra=_extra,color=use_col_sed_filt if not keyword_set(noerrbars) then err_bar,((*!dustem_data.polfrac).wav)(ind_filt_polfrac),((*!dustem_data.polfrac).values)(ind_filt_polfrac),yrms=2*((*!dustem_data.polfrac).sigma)(ind_filt_polfrac),color=use_col_sed_filt print,((*!dustem_data.polfrac).wav)(ind_filt_polfrac),((*!dustem_data.polfrac).values)(ind_filt_polfrac) endif ; PLot color corrected POLFRAC ; VG : P & I have almost identical color correction => no effect on P/I ind_filt=where((*!dustem_data.polsed).filt_names NE 'SPECTRUM',count_filt) if keyword_set(ps) then plotsym,8,thick=8 else plotsym,8 dustem_polsed = dustem_compute_polsed(p_dim,st.polsed) dustem_sed_tmp = dustem_compute_sed(p_dim,st);dustem_compute_sed(p_dim,st,cont=cont) dustem_sed = dustem_polsed for i = 0, count_filt-1 do begin j=where((*!dustem_data.polsed).filt_names(ind_filt(i)) eq (*!dustem_data.sed).filt_names,count) if count eq 1 then dustem_sed(i) = dustem_sed_tmp(j(0)) endfor endif endif endif ;else begin ; endelse ; endelse if keyword_set(ps) && not keyword_set(donotclose) then begin device,/close set_plot,'X' endif the_end: END