dustem_read_grid_table2arrays.pro
3.48 KB
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
PRO dustem_read_grid_table2arrays,input_table_name, $
Ngrains,model,Ncomb,NFparams,INCSP,use_polar,use_double, $
parameter_description,parameter_pmin,parameter_pmax,parameter_pnval,parameter_plog,fparameter_description,fparameter_value, $
param_values,Nparams,seds,filters,Nfilters,wavs,Nwavs, $
I_tot_vec,ext_tot_vec,I_grains_vec,ext_grains_vec,grain_keywords=grain_keywords
fst=file_info(input_table_name)
IF fst.exists NE 1 THEN BEGIN
message,'input table '+input_table_name+' not found',/continue
stop
ENDIF
st0=mrdfits(input_table_name,0,h0) ;This gets the main header
Ngrains=sxpar(h0,'NGRAIN')
grain_keywords=strarr(Ngrains)
FOR i=0L,Ngrains-1 DO BEGIN
grain_keywords[i]=sxpar(h0,'KEYWG'+strtrim(i+1,2),count=count)
IF count EQ 0 THEN grain_keywords[i]='?'
ENDFOR
model=strtrim(sxpar(h0,'MODEL'),2)
Ncomb=sxpar(h0,'NCOMB')
Nparams=sxpar(h0,'NPARAM')
NFparams=sxpar(h0,'NFPARAM')
INCSP=sxpar(h0,'INCSP') ;weither grid includes emission and extinction spectra
use_polar=sxpar(h0,'POLAR')
use_double=sxpar(h0,'DOUBLE')
parameter_description=strarr(Nparams)
parameter_pmin=fltarr(Nparams)
parameter_pmax=fltarr(Nparams)
parameter_pnval=lonarr(Nparams)
parameter_plog=intarr(Nparams)
fparameter_description=strarr(NFparams)
fparameter_value=fltarr(NFparams)
FOR i=0L,Nparams-1 DO BEGIN
parameter_description[i]=sxpar(h0,'PNAME'+strtrim(i+1,2))
parameter_pmin[i]=sxpar(h0,'PMIN'+strtrim(i+1,2))
parameter_pmax[i]=sxpar(h0,'PMAX'+strtrim(i+1,2))
parameter_pnval[i]=sxpar(h0,'PNVAL'+strtrim(i+1,2))
parameter_plog[i]=sxpar(h0,'PLOG'+strtrim(i+1,2))
ENDFOR
FOR i=0L,NFparams-1 DO BEGIN
fparameter_description[i]=sxpar(h0,'FPNAME'+strtrim(i+1,2))
fparameter_value[i]=sxpar(h0,'FPVALUE'+strtrim(i+1,2))
ENDFOR
;=== read parameter values
param_values=mrdfits(input_table_name,1,h1)
Nparams=N_tags(param_values)
;=== read seds
seds=mrdfits(input_table_name,2,h2)
filters=tag_names(seds[0])
Nfilters=n_elements(filters)-1 ;removing the TOTAL column
filters=strmid(filters[0:Nfilters-1],1,10000)
wavs=mrdfits(input_table_name,3,h3)
I_tot=mrdfits(input_table_name,4,h4) ;in MJy/sr
ext_tot=mrdfits(input_table_name,5,h5) ;in MJy/sr
Nwavs=n_elements(wavs.wav)
I_tot_vec=fltarr(Nwavs,Ncomb)
ext_tot_vec=fltarr(Nwavs,Ncomb)
;Here we assume igrain=0, only to create the emission and extinction arrays
I_grains=replicate(I_tot[0],[Ngrains,Ncomb]) ;in MJy/sr for NH=1e20 H/cm2
ext_grains=replicate(ext_tot[0],[Ngrains,Ncomb])
I_grains_vec=fltarr(Nwavs,Ngrains,Ncomb) ;in MJy/sr for NH=1e20 H/cm2
ext_grains_vec=fltarr(Nwavs,Ngrains,Ncomb)
; search for iu, the first fits extension of the grid with per-grain emission spectrum
found_ext=0
kkk=1L
WHILE found_ext EQ 0 DO BEGIN
str=sxpar(h0,'UNIT'+strtrim(kkk,2))
IF str EQ 'Grain1 Emission Spectrum [MJy/sr for NH=1e20 H/cm2]' THEN BEGIN
found_ext=1
iu=kkk
ENDIF ELSE BEGIN
kkk=kkk+1
ENDELSE
ENDWHILE
;stop
FOR i=0L,Ngrains-1 DO BEGIN
aa=mrdfits(input_table_name,iu,hg) ;in MJy/sr for NH=1e20 H/cm2
I_grains[i,*]=aa
iu=iu+1
ENDFOR
;=== This actually assumes that extinction spectra will always be just after emission spectra
FOR i=0L,Ngrains-1 DO BEGIN
aa=mrdfits(input_table_name,iu,hg)
ext_grains[i,*]=aa
iu=iu+1
ENDFOR
;fill in old wavelengths vectors from structures in fits file
FOR k=0L,Nwavs-1 DO BEGIN
I_tot_vec[k,*]=I_tot.(k)
ext_tot_vec[k,*]=ext_tot.(k)
I_grains_vec[k,*,*]=I_grains.(k)
ext_grains_vec[k,*,*]=ext_grains.(k)
ENDFOR
;removing un-needed old variables
I_grains=0
ext_grains=0
ext_tot=0
I_tot=0
;seds=0
END