PRO test_intensity_free_free

;This should reproduce the plot of Fig 2 in Walls-Gabaud PASA 1998, 15 111
nu=range_gen(1000,[1.,1000.])*1.e9 & T=5.e3 & em=1.
I_ff=intensity_free_free(nu,T,em,I_halpha=I_halpha_R)
print,I_halpha_R
;    0.0750229

cgplot,alog10(nu/1.e9),alog10(I_ff/I_halpha_R),yr=[-2.,4],/ysty,xr=[0,3],/xsty,ytit='log(Tb [muK]/IHalpha [R])',xtit='log(nu [GHz])'

T=1.e4
I_ff=intensity_free_free(nu,T,em,I_halpha=I_halpha_R)
cgoplot,alog10(nu/1.e9),alog10(I_ff/I_halpha_R),linestyle=2

T=2.e4
I_ff=intensity_free_free(nu,T,em,I_halpha=I_halpha_R)
cgoplot,alog10(nu/1.e9),alog10(I_ff/I_halpha_R),linestyle=3

;==== test converion to MJy/sr or ergs
nu=range_gen(1000,[1.,1000.])*1.e9 & T=5.e3 & em=1.
I_ff_erg=intensity_free_free(nu,T,em,I_halpha=I_halpha_R,/erg)
print,I_halpha_R

cgplot,alog10(nu/1.e9),alog10(I_ff_erg/I_halpha_R),/ysty,xr=[0,3],/xsty,ytit='log(Tb [erg/s/m2/sr/Hz]/IHalpha [R])',xtit='log(nu [GHz])'
I_ff_mmjjyy=I_ff_erg*1.e20  ;in MJy/sr

I_ff_mjy=intensity_free_free(nu,T,em,I_halpha=I_halpha_R,/mjy)
cgplot,alog10(nu/1.e9),alog10(I_ff_mjy/I_halpha_R),/ysty,xr=[0,3],/xsty,ytit='log(Tb [MJy/sr]/IHalpha [R])',xtit='log(nu [GHz])'
cgoplot,alog10(nu/1.e9),alog10(I_ff_mmjjyy/I_halpha_R),color='red'

END