FUNCTION dustem_read_cb19_stellar_templates,age_values=age_values $
											,metalicity_values=metalicity_values $
											,template_wav=template_wav $
											,Mstars=Mstars $
											,info_only=info_only $
											,info_st=info_st $
											,help=help

;+
; NAME:
;    dustem_read_cb19_stellar_templates
; PURPOSE:
;    reads EMILES stellar population templates
; CATEGORY:
;    DustEMWrap, Distributed, High-Level, User Example
; CALLING SEQUENCE:
;    templates=dustem_read_cb19_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
;    st=dustem_read_cb19_stellar_templates(age_values=age_values,metalicity_values=metalicity_values,template_wav=template_wav)
; 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_cb19_stellar_templates'
  templates=-1
  goto,the_end
END

;=== This is what we want
age_values_wanted=[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_wanted=[-1.48630, -0.961400, -0.351200, +0.0600000, +0.255900, +0.397100]
NagesWanted=n_elements(age_values_wanted)
NZWanted=n_elements(metalicity_values_wanted)

;=== This is what we have
Z=[0.0001, 0.0002, 0.0005, 0.001, 0.002, 0.004, 0.006, 0.008, 0.010, 0.014, 0.017, 0.020, 0.030, 0.040, 0.060]
;la métallicité solaire pour ce set étant: Z=0.01524, et donc ces modèles vont en gros de [Z/Zo]= -2.18 à +0.59
Z0=0.01524
Z=alog10(Z/Z0)
;print,Z
;     -2.18299     -1.88196     -1.48401     -1.18298    -0.881955    -0.580925    -0.404834    -0.279895    -0.182985   -0.0368569    0.0474640     0.118045     0.294136     0.419075     0.595166
NZ=n_elements(Z)

dir_templates=!phangs_data_dir+'/ISRF/SSPs/CB2019/CB19_Chabrier/Mu100/'

;dir_templates=!dustem_wrap_soft_dir+'/Data/SSPs/CB2019/CB19_Chabrier/Mu100/'

files=file_search(dir_templates+'*.fits')
Nfiles=n_elements(files)
;metalicities=fltarr(Nfiles)
;ffile=strarr(Nfiles)
;FOR i=0L,Nfiles-1 DO BEGIN
;	file=files[i]
;	pos=strposmulti(file,'/')
;	ffile[i]=strmid(file[i],pos[n_elements(pos)-1]+1,1000)
;	pos=strposmulti(ffile[i],'_')
;	metalicities[i]=float(strmid(ffile[i],pos[0]+2,pos[1]-pos[0]-2))
;	stop
;	;metalicities[i]=
;ENDFOR

;read one fits file for header infos
file=dir_templates+'cb2019_z001_chab_hr_xmilesi_ssp.fits'
d=mrdfits(file,1,h)
wavs=d.wavelength
tagnames=tag_names(d)
ages=tagnames[1:*]
Nages=n_elements(ages)
ages=strmid(ages,1,1000)
aages=fltarr(Nages)
FOR i=0L,Nages-1 DO BEGIN
	pos=strposmulti(ages[i],'P',count)
	IF count NE 0 THEN BEGIN
		aages[i]=float(textoidl_str_replace(ages[i],'P','.'))/1.e9   ;Gyr
	ENDIF
ENDFOR
ages=aages
Nages=n_elements(ages)

all_ssps=ptrarr(NZ,Nages)
FOR i=0L,Nfiles-1 DO BEGIN
    st=mrdfits(files[i],1,h)
    FOR j=0L,Nages-1 DO BEGIN
    	all_ssps[i,j]=ptr_new(st.(j+1))
    ENDFOR
ENDFOR

;interpolate on requested age an metalicities
Ntemplates=NagesWanted*NZWanted
templates=ptrarr(Ntemplates)
FOR i=0L,Ntemplates-1 DO BEGIN
	ij=index2ij([i],[NagesWanted,NZWanted])
	iage=ij[0,0]
	imet=ij[0,1]
	age=age_values_wanted[iage]
	met=metalicity_values_wanted[imet]
	;==== CAUTION, this is nearest neighbour interpolation
	dist_age=abs(age-ages)
	dist_Z=abs(met-Z)
	ind_age=where(dist_age EQ min(dist_age))
	ind_Z=where(dist_Z EQ min(dist_Z))
	templates[i]=all_ssps[ind_Z[0],ind_age[0]]
	;print,ind_age[0],ind_Z[0]
	;==== This is nearest neighbour on age and interpolation on Z
	ssps_used=all_ssps[*,ind_age[0]]
	xx=indgen(Nz)

  toto=interpol(xx,Z,met)
  u1=abs(toto-round(toto))
  u2=abs(toto-fix(toto))
  ;print,round(toto),toto,fix(toto)
  ;print,1.-u1,1.-u2

  vec=(*ssps_used[round(toto)])*(1.-u1) + (*ssps_used[fix(toto)])*(1.-u2)
  ind=where(vec LT 0,count)
  IF count NE 0 THEN BEGIN
  	stop
  ENDIF
  templates[i]=ptr_new(vec)

	;stop
ENDFOR

age_values=age_values_wanted
metalicity_values=metalicity_values_wanted

;print,age_values
;print,metalicity_values

;stop

angstroem2mic=1.e6/1.e10
template_wav=wavs*angstroem2mic   ;mic

;do_interpol=1L

;IF do_interpol THEN BEGIN
	;stop
	lambir=dustem_get_wavelengths()

	ind_extrapol=where(lambir GT max(template_wav) OR lambir LT min(template_wav),count_extrapol)

	FOR i=0L,Ntemplates-1 DO BEGIN
		;ind=where(lambir LE max(template_wav) AND lambir GE min(template_wav),count)

		vec=interpol(*templates[i],template_wav,lambir)
		vec[ind_extrapol]=0.   ;set extrapolated values to 0 (otherwise would be slightly negative)
		templates[i]=ptr_new(vec)

	ENDFOR
	template_wav=lambir
;ENDIF

the_end:

RETURN,templates

END