FUNCTION dustem_plugin_emiles_stellar_continuum,key=key $
                                			   ,val=val $
                                			   ,scope=scope $
                                			   ,paramtag=paramtag $
                                			   ,paramdefault=paramdefault $
                                			   ,help=help

;+
; NAME:
;    dustem_plugin_emiles_stellar_continuum
; PURPOSE:
;    DustEMWrap plugin to add EMILES stellar template emission to Dustemwrap SEDs
; CATEGORY:
;    DustEM, Distributed, Mid-Level, Plugin
; CALLING SEQUENCE:
;    stellar_brightness=dustem_plugin_emiles_stellar_continuum([,key=][,val=][,scope=][,paramtag=][paramdefault=][,/help])
; INPUTS:
;    None
; OPTIONAL INPUT PARAMETERS:
;    key  = input parameter number
;           First parameter: factor to apply to tau for redenning (1 means value derived from emission SED)
;           Following 13*6 parameters are weights by which to multiply the MILES templates [Msun/pc^2]
;    val  = corresponding input parameter value
; OUTPUTS:
;    stellar_brightness = stellar spectrum spectrum resampled on dustem wavelengths [MJy/sr]
; 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_emiles_stellar_continuum(scope=scope)
; MODIFICATION HISTORY:
;    Written by JPB June 2023
;    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_emiles_stellar_continuum'
  output=0.
  goto,the_end
ENDIF

;stop

;default values of input parameters
scope='ADD_SED'

age_values=[0.03, 0.05, 0.08, 0.15, 0.25, 0.40, 0.60, 1.0, 1.75, 3.0, 5.0, 8.5, 13.5]
metalicity_values=[-1.48630, -0.961400, -0.351200, +0.0600000, +0.255900, +0.397100]

;read template info
info_st=dustem_read_emiles_stellar_templates(age_values,metalicity_values,/info_only)
age_values_emiles=info_st.age
order=sort(age_values_emiles)
age_values_emiles=age_values_emiles[order]
un=uniq(age_values_emiles)
age_values_emiles=age_values_emiles[un]
met_values_emiles=info_st.z
order=sort(met_values_emiles)
met_values_emiles=met_values_emiles[order]
un=uniq(met_values_emiles)
met_values_emiles=met_values_emiles[un]

Nage=n_elements(age_values)
Nmetalicity=n_elements(metalicity_values)
Nparam=Nage*Nmetalicity+1    ;+1 is for opacity
paramtag=strarr(Nparam)
paramdefault=dblarr(Nparam)
paramdefault[*]=0.D0
i=0L
paramtag[0]='tau_f'
ii=1L
FOR i=1L,Nparam-1 DO BEGIN
	ij=index2ij([i-1],[Nage,Nmetalicity])
	paramtag[ii]='Age'+strtrim(ij[0,0],2)+'Met'+strtrim(ij[0,1],2)
	ii=ii+1
ENDFOR
paramvalues=paramdefault
IF keyword_set(key) THEN BEGIN
  ind=where(key EQ 1,count)
  IF count NE 0 THEN paramvalues[0]=val[ind[0]]
  FOR i=1L,Nparam-1 DO BEGIN
  	ind=where(key EQ i+1,count)
  	IF count NE 0 THEN paramvalues[i]=val[ind[0]]
  ENDFOR
ENDIF

;==== initialize templates if not present
defsysv,'!dustem_plugin_emiles_stellar_continuum',exist=exist
IF exist EQ 0 THEN BEGIN
	st={ssps:ptr_new(),wavs:ptr_new(),Mstars:ptr_new()}
	defsysv,'!dustem_plugin_emiles_stellar_continuum',st
ENDIF
IF not ptr_valid(!dustem_plugin_emiles_stellar_continuum.ssps) THEN BEGIN
	message,'Initializing SSP templates',/continue
	template_flux=dustem_read_emiles_stellar_templates(age_values,metalicity_values,template_wav=template_wav,Mstars=Mstars)
	!dustem_plugin_emiles_stellar_continuum.ssps=ptr_new(template_flux)
	!dustem_plugin_emiles_stellar_continuum.wavs=ptr_new(template_wav)
	!dustem_plugin_emiles_stellar_continuum.mstars=ptr_new(Mstars)
ENDIF


lambir=dustem_get_wavelengths()
Nwavs=n_elements(lambir)
output=fltarr(Nwavs,3)
wavs=*!dustem_plugin_emiles_stellar_continuum.wavs
Mstars=*!dustem_plugin_emiles_stellar_continuum.mstars
;==== sum up the stellar spectrum
;stop
FOR i=1L,Nparam-1 DO BEGIN
	;IF paramvalues[i] NE 0. AND ptr_valid((*!phangs_ssps)[i-1]) THEN BEGIN
	IF paramvalues[i] NE 0. THEN BEGIN
		;stop
		;only interpolate within the template wavelengths
		;ind=where(lambir LE max(*!phangs_wavs) AND lambir GE min(*!phangs_wavs),count)
		ind=where(lambir LE max(wavs) AND lambir GE min(wavs),count)
		;Mstar=(*!phangs_Mstars)[i]
		Mstar=Mstars[i-1]
		IF count NE 0 THEN BEGIN
			;output[ind,0]=output[ind,0]+paramvalues[i]*interpol(*(*!phangs_ssps)[i-1],*!phangs_wavs,lambir[ind])
			;spectrum=*(*!phangs_ssps)[i-1]
			;ssp=*(*!phangs_ssps)[i-1]
			ssp=*(*!dustem_plugin_emiles_stellar_continuum.ssps)[i-1]
			;stop
			output[ind,0]=output[ind,0]+paramvalues[i]*Mstar*interpol(ssp,wavs,lambir[ind])
		ENDIF
	ENDIF
ENDFOR
;==== Caution: at this point, output is in Lsun/pc^2/Angstroem
;Caution: I am not sure about the nomalisation of the extinction curve that is in 
IF ptr_valid(!dustem_current) THEN BEGIN
	;stop
	tau_fact=paramvalues[0]
	extinction_curve=((*!dustem_current).ext).ext_tot*(*!dustem_HCD)/1.0e21    ;This is tau
	;*(*!dustem_HCD)/1.0e21
	;wav_ref=551.e-9*1.e6 ;visible 551 nm in microns
	;tau_ref=interpol(extinction_curve,lambir,wav_ref)
	;normalize extinction curve by given tau fact
	extinction_curve=extinction_curve*tau_fact
	;stop
	output[*,0]=output[*,0]*exp(-1.*extinction_curve)
ENDIF

;no polarization for now
output[*,1]=0.
output[*,2]=0. 

;at that stage, output is F_lambda in Lsun/A/pc^2
;goint to MJy/sr (1e-20 W/m2/sr/hz)
pc=3.085677581e16 ;m
;kpc=3.085677581e19  ;m
Lsun=3.828e26  ;Watts
c=3.e5*1e3     ;m/s
;Fact=Lsun       ;Watts/A
;fact=fact*c     ;Watts/hz
;fact=fact/4./!pi  ;Watts/Hz/sr
;fact=fact/pc^2   ;Watts/m2/Hz/sr

fact=Lsun/pc*c/4./!pi/pc
output=output*fact

the_end:
RETURN,output
  
END