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