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

                                			   ,age_values=age_values $
                                			   ,metalicity_values=metalicity_values $

                                			   ,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

;This does not work either
;st_templates=read_muse_templates_info(age_values=age_values,metalicity_values=metalicity_values,Nbins=Nbins,Nage=Nage,NZ=Nmetalicity,bins=bins)
;print,age_values
;     0.030000000     0.050000000     0.080000000      0.15000000      0.25000000      0.40000000      0.60000000       1.0000000       1.7500000       3.0000000       5.0000000       8.5000000       13.500000
;print,metalicity_values
;     -1.49000    -0.960000    -0.350000    0.0600000     0.260000     0.400000

;stop

Nage=n_elements(age_values)
Nmetalicity=n_elements(metalicity_values)

Nparam=Nage*Nmetalicity+2    ;+2 is for Amplitude and opacity

paramdefault=dblarr(Nparam)

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

;==== define parameter tags

paramtag=strarr(Nparam)
paramtag[0]='Amplitude'
paramtag[1]='E(B-V)'
FOR i=2L,Nparam-1 DO BEGIN
	ij=index2ij([i-2],[Nage,Nmetalicity])
	paramtag[i]='Age'+strtrim(ij[0,0],2)+'Met'+strtrim(ij[0,1],2)
ENDFOR

;==== 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=dblarr(Nwavs,3)

wavs=*!dustem_plugin_phangs_stellar_continuum.wavs
Mstars=*!dustem_plugin_phangs_stellar_continuum.mstars
;==== sum up the stellar spectrum
;stop

FOR i=2L,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-2]

		IF count NE 0 THEN BEGIN

			ssp=*(*!dustem_plugin_phangs_stellar_continuum.ssps)[i-2]

			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

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


;==== Caution:
;output is Flambda in ergs/s/cm2/AA
;at this point, output is unredenned 
;==== must be transformed into MJy/sr
NHref=*!DUSTEM_HCD   ;NH dustemwrap reference [H/cm^2]
fact2=lambir*1.e4      ;[AA] to go from F_lambda in ergs/s/cm2/AA to Fnu ergs/s/cm2
fact3=1./NHref         ;[cm2/H] to go from ergs/s/cm2 to ergs/s/H
fact4=1.e4*(*!dustem_HCD)/(4.*!pi)/(3.e8/1.e-6/lambir)*1.e20/1.e7   ;This goes from ergs/s/H to MJy/sr (see dustem_plot_dataset)
fact=fact2*fact3*fact4

print,'flux(1mic)=',interpol(flux,lambir,1.),' ergs/s/cm2/AA'
print,'fact2=',interpol(fact2,lambir,1.)
print,'fact3=',fact3
print,'fact4=',interpol(fact4,lambir,1.)
print,'fact=',interpol(fact,lambir,1.)
print,'flux(1mic)=',interpol(flux,lambir,1.)*interpol(fact,lambir,1.),' MJy/sr'

;goto,no_unred

;==== reden spectrum using Calzetti's extinction law (which is what Phangs does)

;cgplot,lambir,flux,/xlog,/ylog

;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[1], flux_red, R_V = R_V

output[*,0]=flux_red

no_unred:

;fact=Lsun*watts2ergs/Mpc2cm*1.e-4*4.*!pi/Mpc2cm/NHref   ;to go from F_lambda in Lsun/A/pc^2 to Watts/m2/Hz/sr to MJy/sr

;fact=Lsun/pc*(lambir/1.e6)^2/c/4./!pi/pc*1.e20   ;to go from F_lambda in Lsun/A/pc^2 to Watts/m2/Hz/sr to MJy/sr
;print,'flux=',interpol(flux,lambir,1.),' Lsun/Mpc^2/AA'
;print,'fact1=',fact1
;print,'fact2=',interpol(fact2,lambir,1.)
;print,'fact3=',fact3
;print,'fact=',interpol(fact,lambir,1.)
;print,'flux=',interpol(flux,lambir,1.)*interpol(fact,lambir,1.),' ergs/s/H'
;print,'fact4=',interpol(fact4,lambir,1.)
;print,'fact5=',fact5
;print,'flux=',interpol(flux,lambir,1.)*interpol(fact,lambir,1.)*interpol(fact4,lambir,1.),' MJy/sr'

;output=output*fact
output[*,0]=output[*,0]*fact
output[*,1]=output[*,1]*fact
output[*,2]=output[*,2]*fact

print,'output=',interpol(output[*,0],lambir,1.),' MJy/sr'
;stop

the_scope:
scope='ADD_SED'

the_paramtag:

the_end:
RETURN,output
  
END