First commit

Jean-Philippe Bernard
1 parent ab47fb7c
Showing 4 changed files with 817 additions and 0 deletions Show diff stats
src/idl/dustem_fit_intensity_emiles_example.pro
src/idl/dustem_plugin_emiles_stellar_continuum.pro
src/idl/dustem_plugin_phangs_stellar_continuum.pro
src/idl/dustem_read_emiles_stellar_templates.pro
@@ -0,0 +1,319 @@
+PRO dustem_fit_intensity_emiles_example,model=model $
+                                ,sed_file=sed_file $
+                                ,Nitermax=Nitermax $
+                                ,postscript=postscript $
+                                ,fits_save=fits_save $
+                                ,help=help $
+                                ,wait=wait $
+                                ,noobj=noobj $
+                                ,verbose=verbose
+
+;+
+; NAME:
+;    dustem_fit_intensity_emiles_example
+;
+; PURPOSE:
+; This routine is an example of how to fit an observational SED
+; (StokesI only) with DustEM and DustEMWrap, and using the EMiles Strellar template pluggin
+;  
+; For this example, the code uses the SED in the file example_SED_1.xcat,
+; which is distributed in the Data/EXAMPLE_OBSDATA/ directory 
+;
+; The example SED has Stokes I photometric data points from
+; IRAC, MIPS and IRAS. Examples illustrating running DustEMWrap to
+; fit spectral data, polarisation data and extinction data
+; are provided in other _example routines in the src/idl/
+; directory. See the DustEMWrap User Guide for more information.
+;
+; CATEGORY:
+;    DustEMWrap, Distributed, High-Level, User Example
+;
+; CALLING SEQUENCE:
+;    dustem_fit_intensity_emiles_example[,model=][sed_file=][,postscript=][,Nitermax=][,fits_save=][,/help,/wait,/verbose]
+;
+; INPUTS:
+;    None
+;
+; OPTIONAL INPUT PARAMETERS:
+;    None
+;
+; OUTPUTS:
+;    None
+;
+; OPTIONAL OUTPUT PARAMETERS:
+;    Plots, results structure in binary FITS table format
+;
+; ACCEPTED KEY-WORDS:
+;    model = specifies the interstellar dust mixture used by
+;            DustEM. See userguide or dustem_test_model_exists.pro
+;            for more details about available models in current release.
+;    sed_file = string naming the path to text file in .xcat format that
+;          describes the observational SED. If not set, the file
+;          'Data/EXAMPLE_OBSDATA/example_SED_1.xcat' is used.  
+;    postscript = if set, final plot is saved as postscript file
+;    Nitermax = maximum number of fit iterations. Default is 5.
+;    fits_save = if set, save the fit results in a binary
+;               FITS file. 
+;    help      = if set, print this help
+;    wait      = if set, wait this many seconds between each step of
+;                the code (for illustration purposes)
+;    verbose      = if set, subroutines will run in verbose mode
+;    noobj     = if set, runs with no object graphics
+;
+; COMMON BLOCKS:
+;    None
+;
+; SIDE EFFECTS:
+;    None
+;
+; RESTRICTIONS:
+;    The DustEM fortran code must be installed
+;    The DustEMWrap IDL code must be installed
+;
+; PROCEDURES AND SUBROUTINES USED:
+;    
+; EXAMPLES
+;    dustem_fit_intensity_example
+;    dustem_fit_intensity_example,Nitermax=1,fits_save='/tmp/mysavefile.fits'
+;    dustem_fit_intensity_example,model='DBP90'
+;
+; MODIFICATION HISTORY:
+;    Written by JPB Apr-2011
+;    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_fit_intensity_emiles_example'
+  goto,the_end
+END
+
+IF keyword_set(model) THEN BEGIN
+  use_model=strupcase(model)
+ENDIF ELSE BEGIN
+  use_model='DBP90'    ;Example with default keywords uses the DBP90 model
+ENDELSE
+
+exists=dustem_test_model_exists(use_model)
+if exists ne 1 then $
+   message,'Unknown dust model'
+
+use_polarization=0   ; initialize Dustemwrap in no polarization mode 
+use_window=2          ; default graphics window number to use for plotting the results
+use_verbose=0
+if keyword_set(verbose) then use_verbose=1
+use_Nitermax=5        ; maximum number of iterations for the fit
+IF keyword_set(Nitermax) THEN use_Nitermax=Nitermax
+
+dustem_define_la_common
+
+;=== Set the (model-dependent) parameters that you want to fit (pd),
+;=== their initial values (iv)
+;=== and whether they are bounded (ulimed,llimed,llims,ulims).
+;=== Fixed parameters (fpd) and their values (fiv) are also set here.
+;=== Refer to the DustEM and DustEMWrap User guides for an explanation
+;=== of the physical meaning of dust model and plug-in parameters, and
+;=== how to specify them.
+
+;=== Examples are provided for some of the dust models. 
+;=== To try them, uncomment the model that you want to try and re-run
+
+;=== AN EXAMPLE FOR DBP90
+;=== Here we fit the dust abundances of the DBP90 model, the
+;=== intensity of the dust-heating radiation field, as well as a plug-in
+;=== for synchrotron emission
+;=== The free parameters are all lower-bounded at zero.
+;=== use_model='DBP90' ; you should specify this above!
+
+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.49, -0.96, -0.35, +0.06, +0.26, +0.4]
+Nage=n_elements(age_values)
+Nmetalicity=n_elements(metalicity_values)
+ia=where(age_values EQ 0.2500,counta)
+iz=where(metalicity_values EQ -0.35,countb)
+ij=lonarr(1,2)
+ij[0,0]=ia & ij[0,1]=iz
+ind1=ij2index(ij,[Nage,Nmetalicity])+1
+ia=where(age_values EQ 0.05,counta)
+iz=where(metalicity_values EQ +0.06,countb)
+ij=lonarr(1,2)
+ij[0,0]=ia & ij[0,1]=iz
+ind2=ij2index(ij,[Nage,Nmetalicity])+1
+
+; pd = [ $
+;      '(*!dustem_params).G0', $                           ;G0
+;      '(*!dustem_params).grains(0).mdust_o_mh',$          ;PAH0 mass fraction
+;      '(*!dustem_params).grains(1).mdust_o_mh',$          ;VSG mass fraction
+;      '(*!dustem_params).grains(2).mdust_o_mh',$          ;BG mass fraction
+;      'dustem_plugin_emiles_stellar_continuum_'+strtrim(ind1[0]+1,2), $   ;Emiles coefficient for 0.25Gyr,Z/H=-0.35
+;      'dustem_plugin_emiles_stellar_continuum_'+strtrim(ind2[0]+1,2)]   ;Emiles coefficient for 0.05Gyr,Z/H=+0.06
+
+; ; initial values vector for run without synchrotron plugin
+; iv =   [1.6, 2.2e-4, 5.7e-4, 3.4e-3,100.,100.]
+
+pd = [ $
+     '(*!dustem_params).G0', $                           ;G0
+     '(*!dustem_params).grains(0).mdust_o_mh',$          ;PAH0 mass fraction
+     '(*!dustem_params).grains(1).mdust_o_mh',$          ;VSG mass fraction
+     '(*!dustem_params).grains(2).mdust_o_mh',$          ;BG mass fraction
+     'dustem_plugin_emiles_stellar_continuum_1', $       ;tau factor to apply for reddening [-]
+     'dustem_plugin_emiles_stellar_continuum_'+strtrim(ind1[0]+1,2), $   ;Emiles coefficient for 0.25Gyr,Z/H=-0.35
+     'dustem_plugin_emiles_stellar_continuum_'+strtrim(ind2[0]+1,2)]   ;Emiles coefficient for 0.05Gyr,Z/H=+0.06
+
+iv =   [1.6, 2.2e-4, 5.7e-4, 3.4e-3,0.5,33.,31.]
+Npar=n_elements(pd)
+ulimed=replicate(0,Npar)
+llimed=replicate(1,Npar) 
+llims=replicate(1.e-15,Npar)
+ulims=replicate(1.e-15,Npar)
+ulimed[4]=1. & ulims[4]=1.
+fpd=[] & fiv=[]
+
+; pd = [ $
+;      '(*!dustem_params).G0', $                           ;G0
+;      '(*!dustem_params).grains(0).mdust_o_mh',$          ;PAH0 mass fraction
+;      '(*!dustem_params).grains(1).mdust_o_mh',$          ;VSG mass fraction
+;      '(*!dustem_params).grains(2).mdust_o_mh',$          ;BG mass fraction
+;      'dustem_plugin_phangs_stellar_continuum_1', $       ;tau factor to apply for reddening [-]
+;      'dustem_plugin_phangs_stellar_continuum_'+strtrim(ind1[0]+1,2), $   ;Emiles coefficient for 0.25Gyr,Z/H=-0.35
+;      'dustem_plugin_phangs_stellar_continuum_'+strtrim(ind2[0]+1,2)]   ;Emiles coefficient for 0.25Gyr,Z/H=-0.35
+; ; 
+; ; initial values vector for run without synchrotron plugin
+; iv =   [1.6, 2.2e-4, 5.7e-4, 3.4e-3,0.1,100.,100.]
+; Npar=n_elements(pd)
+; ulimed=replicate(0,Npar)
+; llimed=replicate(1,Npar) 
+; llims=replicate(1.e-15,Npar)
+; ulims=replicate(1.e-15,Npar)
+; fpd=[] & fiv=[]
+
+
+if keyword_set(wait) then begin
+   message,'Finished setting dust model and plug-in parameters: '+use_model,/info
+   wait,wait
+end
+
+;== INITIALISE DUSTEM
+dustem_init,model=use_model,polarization=use_polarization
+!dustem_nocatch=1
+!dustem_verbose=use_verbose
+IF keyword_set(noobj) THEN !dustem_noobj=1
+!EXCEPT=2 ; for debugging
+
+;=== READ EXAMPLE SED DATA
+dir=!dustem_wrap_soft_dir+'/Data/EXAMPLE_OBSDATA/'
+file=dir+'example_SED_5.xcat'
+IF keyword_set(sed_file) THEN file=sed_file
+sed=read_xcat(file,/silent)
+
+if keyword_set(wait) then begin
+   message,'Finished reading SED data: '+file,/info
+   wait,wait
+end
+
+;;=== ADJUST THE UNCERTAINTIES FROM WITHIN THE CODE (FOR ILLUSTRATION)
+ind=where(sed.sigmaII LT (0.2*sed.StokesI)^2,count)
+IF count NE 0 THEN sed[ind].sigmaII=(0.2*sed[ind].StokesI)^2
+
+;== SET THE OBSERVATIONAL STRUCTURE
+;== sed is passed twice in the call -- the first occurrence is the SED that you
+;== wish to fit, the second occurrence is the SED that you wish to visualise. 
+dustem_set_data,m_fit=sed,m_show=sed
+
+;== SET INITIAL VALUES AND LIMITS OF THE PARAMETERS THAT WILL BE
+;== ADJUSTED DURING THE FIT
+dustem_init_params,use_model,pd,iv,fpd=fpd,fiv=fiv,ulimed=ulimed,llimed=llimed,ulims=ulims,llims=llims
+
+
+if keyword_set(wait) then begin
+   message,'Finished initializing DustEMWrap, including plugins and fixed parameters',/info
+   wait,wait
+end
+
+;=== INFORMATION TO RUN THE FIT
+tol=1.e-16     ;fit tolerence
+
+;=== INFORMATION TO MAKE THE PLOT
+yr=[1.00e-4,1.00E2] ; y-axis limits
+xr=[1.00E-1,6.00e4] ; x-axis limits
+tit='FIT INTENSITY EXAMPLE' ; plot title
+ytit=textoidl('I_\nu (MJy/sr) for N_H=10^{20} H/cm^2') ; y-axis title
+xtit=textoidl('\lambda (\mum)') ; x-axis title
+
+
+;===  RUN THE FIT
+t1=systime(0,/sec)
+res=dustem_mpfit_data(tol=tol,Nitermax=use_Nitermax,gtol=gtol $
+                      ,/xlog,/ylog,xr=xr,yr=yr,xtit=xtit,ytit=ytit,title=tit $
+                      ,legend_xpos=legend_xpos,legend_ypos=legend_ypos $
+                      ,errors=errors,chi2=chi2,rchi2=rchi2,show_plot=show_plot)
+t2=systime(0,/sec)
+
+if keyword_set(wait) then begin
+   message,'Finished running DustEMWrap, using Niters: '+strtrim(string(use_Nitermax),2),/info
+   message,'Time taken [sec]: '+sigfig(t2-t1,2,/sci),/info
+   wait,wait
+end
+
+;=== MAKE THE FINAL PLOT
+IF keyword_set(postscript) THEN BEGIN
+;    dir_ps='./'
+    mydevice=!d.name
+    set_plot,'PS'
+;    ps_file=dir_ps+postscript
+    ps_file=postscript
+    device,filename=ps_file,/color
+ENDIF
+
+;stop
+
+IF !dustem_noobj THEN BEGIN
+  dustemwrap_plot_noobj,*(*!dustem_fit).CURRENT_PARAM_VALUES,st=dummy,xr=xr,/xstyle,yr=yr,/ysty,/ylog,/xlog,title=tit+' (Final fit)'
+ ENDIF ELSE BEGIN
+  dustemwrap_plot,*(*!dustem_fit).CURRENT_PARAM_VALUES,st=dummy,xr=xr,/xstyle,yr=yr,/ysty,/ylog,/xlog,title=tit+' (Final fit)'
+ENDELSE
+
+IF keyword_set(postscript) THEN BEGIN
+  device,/close
+  set_plot,mydevice
+  message,'Wrote '+ps_file,/info
+ENDIF
+
+if keyword_set(wait) then begin
+   message,'Made the plot of the final results',/info
+   wait,wait
+end
+
+IF keyword_set(fits_save) THEN BEGIN
+   message,'Writing out results structure: '+fits_save,/info
+   dustem_write_fits_table,filename=fits_save,help=help
+;=== At this point, you could erase all dustem system variables, or exit idl... all the
+;=== information needed to recover the results and remake the plots has been saved in the FITS table
+
+;=== Moved the following to dustem_fits_io_example
+  ;; dustem_read_fits_table,filename=fits_save,dustem_st=dustem_spectra_st
+  ;; ;stop
+  ;; ;!dustem_show=ptr_new(*!dustem_data)   ;test
+  ;; !dustem_noobj=1
+  ;; ;==== plot result taken from the saved fits table
+  ;; res=*(*!dustem_fit).CURRENT_PARAM_VALUES
+  ;; IF !dustem_noobj THEN BEGIN
+  ;;   ;dustemwrap_plot_noobj,res,st=dustem_spectra_st,xr=xr,/xstyle,yr=yr,/ysty,/ylog,/xlog,title=tit+' (From Saved FITS file)'
+  ;;   dustemwrap_plot_noobj,res,xr=xr,/xstyle,yr=yr,/ysty,/ylog,/xlog,title=tit+' (From Saved FITS file)'
+  ;; ENDIF ELSE BEGIN
+  ;;   ;dustemwrap_plot,res,st=dustem_spectra_st,xr=xr,/xstyle,yr=yr,/ysty,/ylog,/xlog,title=tit+' (From Saved FITS file)'
+  ;;   dustemwrap_plot,res,xr=xr,/xstyle,yr=yr,/ysty,/ylog,/xlog,title=tit+' (From Saved FITS file)'
+  ;; ENDELSE
+  ;; ;stop
+
+  IF keyword_set(wait) THEN BEGIN
+     message,'Saved the results as FITS in the file: '+fits_save,/info
+     wait,wait
+  ENDIF
+ENDIF
+
+
+the_end:
+message,'Finished dustem_fit_intensity_example',/info
+
+END
@@ -0,0 +1,177 @@
+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
@@ -0,0 +1,165 @@
+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
+
+;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]='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
+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_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])
+		ENDIF
+	ENDIF
+ENDFOR
+;==== Caution: at this point, output is unredenned Flambda in Lsun/pc^2/Angstroem
+
+;==== 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
+;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
@@ -0,0 +1,156 @@
+FUNCTION dustem_read_emiles_stellar_templates,age_values, $
+											  metalicity_values, $
+											  template_wav=template_wav, $
+											  Mstars=Mstars, $
+											  info_only=info_only, $
+											  info_st=info_st, $
+											  help=help
+
+;+
+; NAME:
+;    dustem_read_emiles_stellar_templates
+; PURPOSE:
+;    reads EMILES stellar population templates
+; CATEGORY:
+;    DustEMWrap, Distributed, High-Level, User Example
+; CALLING SEQUENCE:
+;    templates=dustem_read_emiles_stellar_templates(age_values,metalicity_values[,template_waves=][,Mstars=][,/info_only][,info_st=][,/help])
+; INPUTS:
+;    age_values : array of age values
+;    metalicity_values : array of metalicity values
+; OPTIONAL INPUT PARAMETERS:
+;    None
+; OUTPUTS:
+;    templates   : SSP templates (F_lambda) [Lsun Msun^-1 A^-1]
+; OPTIONAL OUTPUT PARAMETERS:
+;    template_wav  = common wavelength array for templates [Angstroem]
+;    Mstars        = Mass of stars [Msun]
+;    info_st       = info structure for the templates
+; ACCEPTED KEY-WORDS:
+;    info_only   : if set, read only the information structure
+;    help        : if set, prints this help
+; COMMON BLOCKS:
+;    None
+; SIDE EFFECTS:
+;    None
+; RESTRICTIONS:
+;    The DustEM fortran code must be installed
+;    The DustEMWrap IDL code must be installed
+; PROCEDURES AND SUBROUTINES USED:
+;    
+; EXAMPLES
+;    
+; 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_read_emiles_stellar_templates'
+  templates=-1
+  goto,the_end
+END
+
+dir_templates='/Users/jpb/Soft_Libraries/dustem-wrapper_idl/Data/SSPs/ssp_emiles-ch/'
+
+;read templates info file
+;IMF slope   Z      Age    Mtotal  M(*+remn)  M*   Mremn   Mgas  M(*+remn)/Lv M*/Lv    Mv          unknown
+file=dir_templates+'templates_info.txt'
+format='A,F,F,F,F,F,F,F,F,F,F,F,F'
+readcol,file,imf,slope,z,age,Mtot,Mstar_rem,Mstar,Mremn,Mgas,Mstar_rem_o_Lv,Mstar_o_Lv,Mv,unknown,format=format
+one_st={imf:'',slope:0.,z:0.,age:0.,Mtot:0.,Mstar_rem:0.,Mstar:0.,Mremn:0.,Mgas:0.,Mstar_rem_o_Lv:0.,Mstar_o_Lv:0.,Mv:0.,unknown:0.}
+info_st=replicate(one_st,n_elements(imf))
+info_st.imf=imf
+info_st.slope=slope
+info_st.z=z
+info_st.age=age
+info_st.mtot=Mtot
+info_st.Mstar_rem=Mstar_rem
+info_st.Mstar=Mstar
+info_st.Mremn=Mremn
+info_st.Mgas=Mgas
+info_st.Mstar_rem_o_Lv=Mstar_rem_o_Lv
+info_st.Mstar_o_Lv=Mstar_o_Lv
+info_st.Mv=Mv
+info_st.unknown=unknown
+IF keyword_set(info_only) THEN BEGIN
+	templates=info_st
+	goto,the_end
+ENDIF
+
+;stop
+
+;file name format is ssp_emiles-ch/Ech1.30Zp0.40T13.5000_iTp0.00_baseFe_linear_FWHM_variable.fits
+;=== producing age and metalicity strings for file names 
+age_values_4str=string(age_values,format='(F7.4)')
+FOR i=0L,n_elements(age_values_4str)-1 DO age_values_4str[i]=textoidl_str_replace(age_values_4str[i],' ','0')
+metalicity_values_4str=round(metalicity_values*100.)/100.   ;This is the array for metalicity values in fits files names
+metalicity_values_4str=string(metalicity_values_4str,format='(F5.2)')
+FOR i=0L,n_elements(metalicity_values_4str)-1 DO metalicity_values_4str[i]=textoidl_str_replace(metalicity_values_4str[i],' ','p')
+FOR i=0L,n_elements(metalicity_values_4str)-1 DO metalicity_values_4str[i]=textoidl_str_replace(metalicity_values_4str[i],'-','m')
+
+;Ntemplates=n_elements(paramvalues)-1
+Ntemplates=n_elements(age_values)*n_elements(metalicity_values)
+Nage=n_elements(age_values)
+Nmetalicity=n_elements(metalicity_values)
+
+templates=ptrarr(Ntemplates)
+
+;file names are like: Ech1.30Zm0.35T00.2500_iTp0.00_baseFe_linear_FWHM_variable.fits
+file_start='Ech1.30'
+file_end='_iTp0.00_baseFe_linear_FWHM_variable.fits'
+ii=0L
+Nage_str=6
+Nmet_str=4
+Mstars=fltarr(Ntemplates)
+files_not_found=['']
+FOR i=1L,Ntemplates DO BEGIN
+	zsign_str='p'
+	ij=index2ij([i-1],[Nage,Nmetalicity])
+	iage=ij[0,0]
+	imet=ij[0,1]
+	age=age_values[iage]
+	met=metalicity_values[imet]
+	met_str=metalicity_values_4str[imet]
+	age_str=age_values_4str[iage]
+	message,'reading template for Age '+age_str+' Met '+met_str,/continue
+	file=dir_templates+file_start+'Z'+met_str+'T'+age_str+file_end
+	fi=file_info(file)
+	IF fi.exists THEN BEGIN
+		message,'Reading SSP '+file,/continue
+		spec=readfits(file,h,2)
+		templates[ii]=ptr_new(spec)
+		make_axis,h,wavs
+		ind=where(info_st.z EQ met AND info_st.age EQ age,count)
+        IF count EQ 0 THEN BEGIN
+        	message,'age,met not found in template info',/continue
+        	stop
+        ENDIF ELSE BEGIN
+        	Mstars[i-1]=info_st[ind[0]].Mstar
+        ENDELSE
+	ENDIF ELSE BEGIN
+		message,'file not found '+file,/continue
+		files_not_found=[files_not_found,file]
+	ENDELSE
+	ii=ii+1
+ENDFOR
+
+IF n_elements(files_not_found) NE 1 THEN BEGIN
+	files_not_found=files_not_found[1:*]
+	message,strmid(n_elements(files_not_found),2)+' SSP files were not found :',/continue
+	hprint,files_not_found
+	stop
+ENDIF
+
+;Nspec=n_elements(spec)
+;print,min(wavs),max(wavs)
+;       3500.2000       12000.250
+angstroem2mic=1.e6/1.e10
+template_wav=wavs*angstroem2mic   ;mic
+
+the_end:
+
+RETURN,templates
+
+END
 \ No newline at end of file