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