dustem_plot_polsed.pro
3.58 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
PRO dustem_plot_polsed, st, p_dim, dustem_polsed, aligned=aligned, win=win
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,title='DUSTEM WRAP (POLSED)'
ENDELSE
;I will try to use the refreshed version of the polsed model spectrum
;This is because the default procedures for dustem
titstq='Polarization Intensity'
ytitstq=textoidl('P_\nu (W/m^2/Hz/sr) for N_H=10^{20} H/cm^2')
xtit=textoidl('\lambda (\mum)')
xr=[10,2e4]
yr=[1e-30,1e-19]
Ngrains=(*!dustem_params).Ngrains
Nwaves=(size(st.polsed.em_tot))[1]
if not keyword_set(st) then begin
;activate the plugins here
st=dustem_run(p_dim)
dustem_polsed = dustem_compute_polsed(p_dim,st)
ENDIF
fact=1/(4.*!pi)*(*!dustem_HCD)*1.E-20/(3.e8/1.e-6/(st.polsed).wav)*1.00e17 ;applied correction here to take into account the HCD
fact1=1.E-20 ;This will convert data (real data) that is in MJy/sr to W
P = st.polsed.em_tot*fact ;This is Polarized intensity P
I = st.sed.em_tot*fact;This is Total intensity I
out=0.
frac=P/I
tes=where(finite(frac) eq 0)
frac(tes)=0.
colors=dustem_grains_colors(Ngrains,ls,ps=ps)
scopes=tag_names((*!dustem_scope))
IF scopes[0] NE 'NONE' THEN BEGIN
FOR i=0L,n_tags(*!dustem_scope)-1 DO BEGIN
IF total(strsplit((*(*!dustem_scope).(i)),'+',/extract) EQ 'REPLACE_QSED') THEN BEGIN
Q_spec=(*(*!dustem_plugin).(i))[*,1]
ENDIF
IF total(strsplit((*(*!dustem_scope).(i)),'+',/extract) EQ 'REPLACE_USED') THEN BEGIN
U_spec=(*(*!dustem_plugin).(i))[*,2]
ENDIF
ENDFOR
IF ~isa(Q_spec) && ~isa(U_spec) THEN polar_ippsi2iqu,I,Q_spec,U_spec,frac,replicate(0.,Nwaves)
FOR i=0L,n_tags(*!dustem_scope)-1 DO BEGIN
IF total(strsplit((*(*!dustem_scope).(i)),'+',/extract) EQ 'ADD_QSED') THEN BEGIN
Q_spec+=(*(*!dustem_plugin).(i))[*,1]
ENDIF
IF total(strsplit((*(*!dustem_scope).(i)),'+',/extract) EQ 'ADD_USED') THEN BEGIN
U_spec+=(*(*!dustem_plugin).(i))[*,2]
ENDIF
ENDFOR
P=sqrt(Q_spec^2+U_spec^2)
ENDIF
normpol1=dustem_polsed*fact1
normpol=P
;normpol=P*0.+1. ;and then use this to plot the normalized data
; 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
ylog=1.
cgplot,st.polsed.wav,P,xtit='',ytit=ytitstq,tit=titstq,ylog=ylog,/xlog,xr=xr,yr=yr,/ys,/xs,position=[0.17,0.35,0.93,0.93],xtickformat='(A1)',charsize=1.3
cgoplot,st.polsed.wav,P,color='black'
cgoplot,(*!dustem_data.polsed).wav,(*!dustem_data.polsed).values*fact1,color='Tomato',psym=16,symsize=1,thick=2
err_bar,(*!dustem_data.polsed).wav,(*!dustem_data.polsed).values*fact1,yrms=3.*((*!dustem_data.polsed).sigma)/2.*fact1,color=cgColor('Tomato')
cgoplot,(*!dustem_data.polsed).wav,dustem_polsed*fact1,psym=6,color='red',symsize=2
cgplot,st.polsed.wav,P/P,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,ytickformat='(F6.2)',charsize=1.3
cgoplot,st.polsed.wav,P/P,color='black'
cgoplot,(*!dustem_data.polsed).wav,(*!dustem_data.polsed).values*fact1/normpol1,color='Tomato',psym=16,symsize=1,thick=2
err_bar,(*!dustem_data.polsed).wav,(*!dustem_data.polsed).values*fact1/normpol1,yrms=3.*((*!dustem_data.polsed).sigma)/2./normpol1*1E-20,color=cgColor('Tomato')
;stop
IF keyword_set(ps) THEN BEGIN
device,/close
set_plot,'X'
ENDIF
END