dustem_plugin_freefree.pro
4.37 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
FUNCTION dustem_plugin_freefree,key=key $
,val=val $
,scope=scope $
,paramtag=paramtag $
,paramdefault=paramdefault $
,help=help
;+
; NAME:
; dustem_plugin_freefree
;
; PURPOSE:
; DustEMWrap plugin to compute free-free emission
;
; CATEGORY:
; DustEM, Distributed, Mid-Level, Plugin
;
; CALLING SEQUENCE:
; freefree=dustem_plugin_freefree([,key=][,val=])
;
; INPUTS:
; None
;
; OPTIONAL INPUT PARAMETERS:
; key = input parameter number
; 1 : Tgas [K]
; 2 : free-free Amplitude []
; 3 : polarization fraction [%] (default=0)
; 4 : polarization angle [deg]
; val = input parameter value
;
; OUTPUTS:
; freefree = free-free spectrum (on dustem wavelengths)
;
; OPTIONAL OUTPUT PARAMETERS:
; scope = scope of the plugin
; paramdefault = default values of parameters
; paramtag = plugin parameter names as strings
;
; ACCEPTED KEY-WORDS:
; help = if set, print this help
;
; COMMON BLOCKS:
; None
;
; SIDE EFFECTS:
; None
;
; RESTRICTIONS:
; The DustEMWrap IDL code must be installed
;
; PROCEDURE:
; This is a DustEMWrap plugin
;
; EXAMPLES
; dustem_init
; vec=dustem_plugin_freefree(scope=scope)
;
; MODIFICATION HISTORY:
; Written by JPB 2022
; Evolution details on the DustEMWrap gitlab.
; See http://dustemwrap.irap.omp.eu/ for FAQ and help.
;-
IF keyword_set(help) THEN BEGIN
doc_library,'dustem_plugin_freefree'
output=0.
goto,the_end
ENDIF
;default values of input parameters
Tgas=10000. ;default gas temperature
Amplitude=1. ;Amplitude
smallp=0.0 ;default polarization fraction
psi=0. ;default polarization angle
scope='ADD_SED'
paramtag=[textoidl('T_{gas}')+' [K]','Amp','p','Psi [deg]']
paramdefault=[Tgas,Amplitude,smallp,psi]
IF keyword_set(key) THEN BEGIN
a=where(key EQ 1,count1)
b=where(key EQ 2,count2)
c=where(key EQ 3,count3) ;default polarization fraction -newly added
d=where(key EQ 4,count4) ;default polarization angle -newly added
IF count1 NE 0 then Tgas=val[a[0]]
IF count2 NE 0 then Amplitude=val[b[0]]
IF count3 NE 0 then smallp=val[c[0]]
IF count4 NE 0 then psi=val[d[0]]
ENDIF
lambir=dustem_get_wavelengths()
Nwavs=n_elements(lambir)
cmic=3.e14
nu=cmic/lambir ;Hz
mjy=1 ;output is in MJy/sr
lambir_ref=10000. ;in mic (this is 1cm)
output=fltarr(Nwavs,3)
;use_method='Deschenes2008'
;use_method='WallsGabaud1998'
use_method='Dickinson2003_norm'
CASE use_method OF
'WallsGabaud1998':BEGIN
em=1. ;This is a stupid value for the Emission measure. Result will be rescaled based on I_halpha_R anyway
output[*,0]=intensity_free_free(nu,Tgas,em,I_halpha_R=I_halpha_R,nHe=nHe,ergs=ergs,mjy=mjy)
fact=Amplitude/I_halpha_R
;=== normalize to the requested Halpha value.
output[*,0]=output[*,0]*fact
END
'WallsGabaud1998_scaled':BEGIN
em=1. ;This is a stupid value for the Emission measure. Result will be rescaled based on I_halpha_R anyway
output[*,0]=intensity_free_free(nu,Tgas,em,I_halpha_R=I_halpha_R,nHe=nHe,ergs=ergs,mjy=mjy)
norm=interpol(output[*,0],lambir,lambir_ref)
output[*,0]=output[*,0]/norm*Amplitude
END
'Deschenes2008':BEGIN
beta_freefree=2.+1./(10.48+1.5*alog(Tgas/8.e3)-alog(nu/1.e9))
output[*,0]=nu^(-1.*beta_freefree)
norm=interpol(output[*,0],lambir,lambir_ref)
output[*,0]=output[*,0]/norm*Amplitude
END
'Dickinson2003_norm':BEGIN
;stop
em=1.
T4=Tgas/1.e4
nu_ghz=nu/1.e9
a=0.366*(nu_ghz/1.e9)^(-0.15)*(alog(4.995*1e-2*(nu_ghz)^(-1.))+1.5*alog(Tgas))
Tb=8.369e3*a*(nu_ghz)^(-2.)*T4^0.667*10^(0.029*T4)*(1.+0.08) ;in mK for 1 Rayleigh
psm_convert_mk2mjy, lambir, Tb, I_Mjy, /RJ
norm=interpol(I_Mjy,lambir,lambir_ref)
output[*,0]=I_Mjy/norm*Amplitude
END
ENDCASE
polar_ippsi2iqu,output[*,0],Q,U,replicate(smallp,Nwavs),replicate(psi,Nwavs)
output[*,1]=Q
output[*,2]=U
the_end:
RETURN,output
END