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

;+
; NAME:
;    dustem_plugin_phangs_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_phangs_stellar_continuum([,key=][,val=][,scope=][,paramtag=][paramdefault=][,/help])
; INPUTS:
;    None
; OPTIONAL INPUT PARAMETERS:
;    key  = input parameter number
;           First parameter: E(B-V) for extinction applied to the template, as used in Phangs
;           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 for phangs
;    It differs from dustem_plugin_emiles_stellar_continuum.pro only by the way the extinction is applied to the output spectrum
; EXAMPLES
;    dustem_init
;    vec=dustem_plugin_phangs_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_phangs_stellar_continuum'
  output=0.
  goto,the_end
ENDIF

IF keyword_set(scope) THEN BEGIN    
    out=0
    goto, the_scope
ENDIF 

;IF keyword_set(paramtag) THEN BEGIN   
;    out=0
;    goto, the_paramtag
;ENDIF

;stop

;default values of input parameters

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)

;stop

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]

;==== This does not work because the info files has many more entries than there are available ssps.
;age_values=age_values_emiles
;metalicity_values=met_values_emiles

Nage=n_elements(age_values)
Nmetalicity=n_elements(metalicity_values)
Nparam=Nage*Nmetalicity+1    ;+1 is for opacity

paramdefault=dblarr(Nparam)

paramdefault[*]=0.D0         ;default parameter values is 0.

;==== define parameter tags

;If the test below is activated (ie paramtags computed only when requested) the pluggin does not work.
;IF keyword_set(paramtag) THEN BEGIN
	paramtag=strarr(Nparam)
	paramtag[0]='E(B-V)'
	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
    ;out=0.
    ;GOTO,the_paramtag
;ENDIF

;==== decode key keyword. 
;stop

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 already present

defsysv,'!dustem_plugin_phangs_stellar_continuum',exist=exist
IF exist EQ 0 THEN BEGIN
	st={ssps:ptr_new(),wavs:ptr_new(),Mstars:ptr_new()}
	defsysv,'!dustem_plugin_phangs_stellar_continuum',st
ENDIF
IF not ptr_valid(!dustem_plugin_phangs_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_phangs_stellar_continuum.ssps=ptr_new(template_flux)
	!dustem_plugin_phangs_stellar_continuum.wavs=ptr_new(template_wav)
	!dustem_plugin_phangs_stellar_continuum.mstars=ptr_new(Mstars)
ENDIF

;stop

lambir=dustem_get_wavelengths()
Nwavs=n_elements(lambir)
output=fltarr(Nwavs,3)
wavs=*!dustem_plugin_phangs_stellar_continuum.wavs
Mstars=*!dustem_plugin_phangs_stellar_continuum.mstars
;==== sum up the stellar spectrum
;stop
FOR i=1L,Nparam-1 DO BEGIN
	IF paramvalues[i] NE 0. THEN BEGIN
		;stop
		;only interpolate within the template wavelengths
		ind=where(lambir LE max(wavs) AND lambir GE min(wavs),count)
		Mstar=Mstars[i-1]
		IF count NE 0 THEN BEGIN
			ssp=*(*!dustem_plugin_phangs_stellar_continuum.ssps)[i-1]
			output[ind,0]=output[ind,0]+paramvalues[i]*Mstar*interpol(ssp,wavs,lambir[ind])

			message,'Mstar='+strtrim(Mstar,2)+' paramvalue='+strtrim(paramvalues[i],2),/continue

		ENDIF
	ENDIF
ENDFOR
;==== Caution: at this point, output is unredenned Flambda in Lsun/pc^2/Angstroem

stop

;==== reden spectrum using Calzetti's extinction law (which is what Phangs does)
flux=output[*,0]
;factor -1 is to reden, instead of unreden. factor 1e4 is to go from microns to Angstroem
calz_unred, lambir*1.e4, flux, -1.*paramvalues[0], flux_red, R_V = R_V
output[*,0]=flux_red

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

;at that stage, output is F_lambda in Lsun/A/pc^2

;going to W/m2/sr/hz

pc=3.085677581e16 ;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   ;to go from F_lambda in Lsun/A/pc^2 to Watts/m2/Hz/sr

output=output*fact

;print,fact
;      9.59803

the_scope:
scope='ADD_SED'

the_paramtag:

the_end:
RETURN,output
  
END