FUNCTION dustem_plugin_phangs_stellar_continuum,key=key $

                                			   			 ,val=val $
                                			   			 ,scope=scope $
                                			   			 ,paramtag=paramtag $
                                			   			 ,age_values=age_values $
                                			   			 ,metalicity_values=metalicity_values $
                                			   			 ,paramdefault=paramdefault $

;                                			   			 ,object_distance=object_distance $    ;not used

                                			   			 ,Voronoi_Npix=Voronoi_Npix $
                                			   			 ,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_continuum=dustem_plugin_phangs_stellar_continuum([,key=][,val=][,scope=][,paramtag=][paramdefault=][,/help])

; INPUTS:
;    None
; OPTIONAL INPUT PARAMETERS:
;    key  = input parameter number

;					  First parameter: Amplitude factor for the emission
;           Second parameter: E(B-V) for extinction to stars 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_continuum = stellar continuum spectrum resampled on dustem wavelengths [MJy/sr]

; OPTIONAL OUTPUT PARAMETERS:

;    scope = if set, return the scope of the plugin

;    paramdefault = default values of parameters

;    paramtag = if set, return the return plugin parameter names as strings

; ACCEPTED KEY-WORDS:
;    help                  = if set, print this help

;    method                = SSP used. can be EMILES or CB19. default='EMILES'

; 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
    scope='ADD_SED' 

    out=0
    goto, the_scope
ENDIF 

;use_method='EMILES'
use_method='CB19'

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]

;==== 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 E(B-V)

;==== define parameter tags

IF keyword_set(paramtag) THEN BEGIN
	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
	GOTO,the_paramtag
ENDIF

paramdefault=dblarr(Nparam)
paramdefault[*]=0.D0         ;default parameter values is 0.

;==== 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(),mstars_o_Lv:ptr_new(),Mstarprems: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',/info
	CASE use_method of
		'EMILES': BEGIN

			template_flux=dustem_read_emiles_stellar_templates(age_values,metalicity_values,template_wav=template_wav,Mstars=Mstars,Ms_o_Lv=Ms_o_Lv,Mstarprems=Mstarprems)

		END
		'CB19': BEGIN
			;== this is just to get the Mstars

			bidon=dustem_read_emiles_stellar_templates(age_values,metalicity_values,template_wav=template_wav,Mstars=Mstars,Ms_o_Lv=Ms_o_Lv,Mstarprems=Mstarprems)

			;== This is to read in the CB19 templates

			;== Note, templates are already interpolated on dustem wavelenghts (!)

			template_flux=dustem_read_cb19_stellar_templates(age_values=age_values,metalicity_values=metalicity_values,template_wav=template_wav);,Mstars=Mstars)

		END
	ENDCASE

	!dustem_plugin_phangs_stellar_continuum.ssps=ptr_new(template_flux)   ; in Lsun Msun^-1 A^-1

	!dustem_plugin_phangs_stellar_continuum.wavs=ptr_new(template_wav)
	!dustem_plugin_phangs_stellar_continuum.mstars=ptr_new(Mstars)

	!dustem_plugin_phangs_stellar_continuum.mstars_o_Lv=ptr_new(Ms_o_Lv)
	!dustem_plugin_phangs_stellar_continuum.mstarprems=ptr_new(Mstarprems)

ENDIF

;stop

use_Voronoi_Npix=1.   ;deg^2
IF keyword_set(Voronoi_Npix) THEN use_Voronoi_Npix=Voronoi_Npix
;use_object_distance=1.   ;MPc
;IF keyword_set(object_distance) THEN use_object_distance=object_distance

;bin_surface_pc2=use_bin_surface*!pi^2/180.^2*(use_object_distance*1.e6)^2   ;surface in pc^2
;Npixfact=(1./(0.2/60./60.))^2        ;how much Muse pixels of 0.2 arcsec in 1 deg^2
;bin_surface_Npixmuse=use_bin_surface*Npixfact

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

Ms_o_Lv=*!dustem_plugin_phangs_stellar_continuum.mstars_o_Lv
Mstarprems=*!dustem_plugin_phangs_stellar_continuum.Mstarprems

;==== 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]

		Mstarprem=Mstarprems[i-2]
		M_o_Lv=Ms_o_Lv[i-2]

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

		;output[*,0]=output[*,0]+paramvalues[i]*Mstar*ssp
		;==== below is trying to remove Mstar

		;output[*,0]=output[*,0]+paramvalues[i]*ssp
		;==== below is trying to multiply weights by L/M
		;output[*,0]=output[*,0]+paramvalues[i]/M_o_Lv*ssp*use_bin_surface
		;output[*,0]=output[*,0]+paramvalues[i]/M_o_Lv*ssp/use_bin_surface
		;output[*,0]=output[*,0]+paramvalues[i]*ssp*Mstar/(bin_surface_pc2)
		;output[*,0]=output[*,0]+paramvalues[i]*ssp*Mstar/(bin_surface_Npixmuse)
;		output[*,0]=output[*,0]+paramvalues[i]*ssp*Mstar/(use_Voronoi_Npix)
		output[*,0]=output[*,0]+paramvalues[i]*ssp*Mstarprem/(use_Voronoi_Npix)

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

	ENDIF
ENDFOR

output[*,0]=output[*,0]*paramvalues[0]

;cgplot,wavs*1.e4,output[*,0],xrange=[4750,7000]
;stop

;flux=output[*,0]

;no polarization for now

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

;==== Caution:

;output is Flambda in Lsun/pc2/AA

;at this point, output is unredenned 
;==== must be transformed into MJy/sr

Lsun2ergs=3.828e33
pc2_to_cm2=(3.0857e18)^2
fact1=Lsun2ergs/pc2_to_cm2

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
fact=fact1*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

;This uses the prescription of Calzetti et al 2000. It is applied from 912AA to 2.2 mic, no correction applied before and beyound that.
;Note that the supplied color excess should be that derived for the stellar  continuum, EBV(stars), which (in Calzetti 2000) is related to the 
;reddening derived from the gas, EBV(gas), via the Balmer decrement by EBV(stars) = 0.44*EBV(gas)

calz_unred, lambir*1.e4, output[*,0], -1.*paramvalues[1], flux_red, R_V = R_V

;Note: we could actually use our own dustemwrap extinction curve to do that.

output[*,0]=flux_red

no_unred:

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:

the_paramtag:

the_end:
RETURN,output
  
END