dustem_plugin_freefree.pro
3.63 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
FUNCTION dustem_plugin_freefree,key=key,val=val,scope=scope,paramtag=paramtag,help=help
;+
; NAME:
; dustem_plugin_freefree
; PURPOSE:
; DUSTEM plugin to compute free-free emission
; CATEGORY:
; DUSTEM Wrapper
; CALLING SEQUENCE:
; freefree=dustem_plugin_freefree([,key=][,val=])
; INPUTS:
; None
; OPTIONAL INPUT PARAMETERS:
; key = input parameter number
; val = input parameter value
; OUTPUTS:
; freefree = free-free spectrum (on dustem wavelengths)
; OPTIONAL OUTPUT PARAMETERS:
; scope = scope of the plugin
; ACCEPTED KEY-WORDS:
; help = if set, print this help
; COMMON BLOCKS:
; None
; SIDE EFFECTS:
; None
; RESTRICTIONS:
; None
; PROCEDURE:
; This is a dustem plugin
; EXAMPLES
; dustem_init
; vec=dustem_plugin_freefree(scope=scope)
; MODIFICATION HISTORY:
; Written by JPB
;-
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
paramtag=['T_gas (K)','Amp','p','Psi (deg)']
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]] ;-newly added
IF count4 NE 0 then psi=val[d[0]] ;-newly added
ENDIF
;IF !dustem_which EQ 'DESERT' THEN BEGIN
; lambir=((*!dustem_params).gemissiv.lambir)
;ENDIF ELSE BEGIN
; lambir=((*!dustem_params).lambda.lambda)
;ENDELSE
;replaced by the following to ease the life of plugins writters
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.
output=fltarr(Nwavs,5) ;newly added
;stop
;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
convert_mk_mjy, lambir, Tb, I_Mjy, /RJ
norm=interpol(I_Mjy,lambir,lambir_ref)
output[*,0]=I_Mjy/norm*Amplitude
END
ENDCASE
scope='ADD_SED'
;polarization this actually need to
polar_ippsi2iqu,output[*,0],Q,U,replicate(smallp,Nwavs),replicate(psi,Nwavs)
output[*,1]=Q
output[*,2]=U
output[*,3]=replicate(smallp,Nwavs)
output[*,4]=replicate(psi,Nwavs)
the_end:
RETURN,output
END