OV-GSO-DC / dustem-wrapper_idl

Blame view

src/idl/dustem_plugin_ff_hii.pro 4.25 KB
Edit Raw Normal View History
93b2bb2f   Ilyes Choubani   other freefree pl... 
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
 
FUNCTION dustem_plugin_ff_hii,key=key,val=val,scope=scope,paramtag=paramtag,paramdefault=paramdefault,help=help

;+
; NAME:
;    dustem_plugin_ff_hii
; 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
;    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_ff_hi'
  output=0.
  goto,the_end
ENDIF

;default values of input parameters
Te=1.0E4             ;default electron temperature (K)
;Adding an amplitude just for the sake of it (still not done)
;Amplitude=       ;Amplitude for Te=1.E4 and nu=  GHz (arb) and lambda=1E4 microns
;Not so sure about this emission measure
EM = 100.             ;default emission measure  (cm^-6 pc)
smallp=0.0            ;default polarization fraction
psi=0.                ;default polarization angle
scope='ADD_SED'
paramtag=[textoidl('T_{e}')+' [K]','EM','p','Psi [deg]']
paramdefault=[Te,EM,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 Te=(val(a))(0)
  IF count2 NE 0 then EM=(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
k_B = 1.380649E-23 ;J.K^-1 
g_ff = alog(4.955*1D-2/nu*1D9)+1.5*alog(Te) ;gaunt factor
tau_ff = 3.014*1D-2*Te^(-1.5)*(nu*1D-9)^(-2)*EM*g_ff ;ff optical depth 
T_ff = Te*(1-exp(-tau_ff))

out=fltarr(Nwavs,3)
out[*,0] = 2*k_B*T_ff*nu^2/(cmic*1E-6)^2*1E20 ;convert to MJy by multiplying by 1E20

polar_ippsi2iqu,out[*,0],Q,U,replicate(smallp,Nwavs),replicate(psi,Nwavs)

out[*,1]=Q
out[*,2]=U

; ;Not so sure about these two: 
; mjy=1           ;output is in MJy/sr
; lambir_ref=10000.
; 
; 
; output=fltarr(Nwavs,3) ;newly added
; 
; 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

; polar_ippsi2iqu,output[*,0],Q,U,replicate(smallp,Nwavs),replicate(psi,Nwavs)
; 
; output[*,1]=Q
; output[*,2]=U

the_end:
RETURN,out
  
END