OV-GSO-DC / dustem-wrapper_idl

Blame view

src/idl/dustem_read_cb19_stellar_templates.pro 4.66 KB
Edit Raw Normal View History
a22638d5   Jean-Philippe Bernard   first commit 
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
 
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=!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
    templates[i]=ptr_new((*ssps_used[round(toto)])*(1.-u1) + (*ssps_used[fix(toto)])*(1.-u2))
	;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

the_end:

RETURN,templates

END