intensity_free_free.pro
1.45 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
FUNCTION intensity_free_free,nu,T,em,I_halpha_R=I_halpha_R,nHe=nHe,ergs=ergs,mjy=mjy
;nu=frequency in Hz
;em is emission measure in cm^-6*pc
;T is the gas temperature
;example:
;nu=range_gen(1000,[1.,1000.])*1.e9 & T=5.e4 & em=1.
;I_ff=intensity_free_free(nu,T,em,I_halpha=I_halpha_R)
;cgplot,nu/1.e9,I_ff/I_halpha_R,/xlog,/ylog,yr=[1.e-2,1.e4],/ysty,xr=[1,1e3],/xsty
;print,I_halpha_R
; 0.0750229
nu_GHz=nu/1.e9
nu_10=nu_GHz/10. ;in units of 10 GHz
c=2.998e8 ;m/s
k=1.38e-23 ;J/K
;compute I Halpha for caseB recombination in Rayleigh
I_halpha=intensity_h_alpha(em,T) ;this is in erg/cm2/s/sr
R=5.661e-18 ;1R in ergs/s/cm^2/arcsec^2
sr=4.254517e10 ;1sr in arcsec^2
RR=R*sr
I_halpha_R=I_halpha/RR ;in Rayleigh
Z=1.
;stop
T4=T/1.e4
;ionization potentials of He
chi_he=24.5874 ;eV
chi_Heii=54.41778 ;eV
use_nHe=0.
IF keyword_set(nhe) THEN use_nHe=nHe
;not sure what the double 0.08 really means ....
T_muK=I_halpha_R*14.*T4^0.317*10^(0.029/T4)*(1.+0.08*chi_Heii/chi_he*use_nHe/0.08)*Z^2*gaunt_ff(nu,T)/nu_10^2
;T_muK=I_halpha_R*14.*T4^0.317*10.^(0.029/T4)*Z^2*gaunt_ff(nu,T)/nu_10^2
;T_muK is the free-free brightness temperature (I=2*h*nu^2*k*Tb/c^2) in muK
Ifreefree=T_muK
IF keyword_set(ergs) THEN BEGIN
Ifreefree=2.*nu^2*k*T_muK*1.e6/c^2 ;This is ergs/s/m^2/Hz/sr
ENDIF
;convert to MJy/sr if needed
IF keyword_set(mjy) THEN BEGIN
lambda=c/nu*1.e6 ;mic
T_mK=T_muK/1.e3
convert_mk_mjy, c/nu, T_mK, I_Mjy, /RJ
Ifreefree=I_Mjy
ENDIF
RETURN,Ifreefree
END