dustem_str_tools.pro
6.58 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
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
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
FUNCTION ADD_MOD, sf, tag, nsz, unit=unit
; adds a model tag to a existing SMDAT structure
; SF (I): input SMDAT structure
; TAG (I): tag of model to be added
; NSZ (I): int array(2) of sizes [nr of data pts, nr of grain types]
tag = STRLOWCASE(STRTRIM(tag,2))
if n_elements(nsz) GT 1 then begin
n1 = nsz(0) & n2 = nsz(1)
if n_elements(nsz) EQ 3 then n3 = nsz(2) else n3 = 0
endif else begin
print,'(F) ADD_MOD: array of sizes must > 1D'
endelse
if TAG EQ 'emis' then begin
if n_elements(unit) EQ 0 then unit='x(microns) SED(erg/s/cm2/sr)'
s1 = { UNIT : unit, $
X : dblarr(n1), $ ; wave
Y : dblarr(n1,n2), $ ; SED(wave, grain type) (index ntype is total)
YP : dblarr(n1,n2) } ; polarized SED(wave, grain type) (index ntype is total)
ENDIF else if TAG EQ 'ext' then begin
if n_elements(unit) EQ 0 then unit='x(microns) sigma(cm2/H)'
s1 = {UNIT : unit, $
X : dblarr(n1), $ ; wave
Y : dblarr(n1,n2), $ ; sigma_ext(wave, grain type) (index ntype is total)
ABS : dblarr(n1,n2), $ ; sigma_abs(wave, grain type)
SCA : dblarr(n1,n2), $ ; sigma_sca(wave, grain type)
ABS_P : dblarr(n1,n2), $ ; absorption sigma_pol(wave, grain type)
SCA_P : dblarr(n1,n2), $ ; scattering sigma_pol(wave, grain type)
ALB : dblarr(n1,n2), $ ; alb(wave, grain type)
XR : 0d, $ ; ref wave
YR_ABS : dblarr(n2), $ ; tau_abs/NH @ XR
YR_SCA : dblarr(n2), $ ; tau_sca/NH @ XR
RV : 0d }
; ENDIF else if TAG EQ 'pol' then begin
; if n_elements(unit) EQ 0 then unit='x(microns) SED(erg/s/cm2/sr)'
; s1 = {UNIT : unit, $
; X: dblarr(n1), $ ; wave
; Y: dblarr(n1,n2), $ ; polarized SED(wave, grain type) (index ntype is total)
; POL: dblarr(n1,n2) } ; sigma_pol(wave, grain type) (ABS + SCAT)
ENDIF else if tag EQ 'sdist' then begin
if n3 EQ 0 then begin
print,'(F) ADD_MOD: 3d dimension missing '
return,0
endif
if n_elements(unit) EQ 0 then unit='x(cm) a^4*dn/da(cm3/H)'
s1 = {UNIT : unit, $
XTOT : dblarr(n1), $ ; grain size
YTOT : dblarr(n1,n2), $ ; size distribution
XI : dblarr(n3,n2), $ ; grain size per type
YI : dblarr(n3,n2) } ; grain size dist per type
ENDIF
sm = CREATE_STRUCT( CREATE_STRUCT('M_'+tag,s1), sf )
return, sm
END
FUNCTION STR_INST, n1, n2=n2, n3=n3
; returns structure containing flux and CC in band
; N1 (I): number of data points
; N2 (I): nr of grain types or different models
; N3 (I): nr of points in transmission (band flux data)
if n_elements(n2) EQ 0 then n2 = 1
if n_elements(n3) NE 0 then begin
strct = { NAME : strarr(n1), $
X : dblarr(n1), $
YD : dblarr(n1), $
ERR : dblarr(n1), $
YM : dblarr(n1,n2), $
FLX : dblarr(n1,n2), $
CC : dblarr(n1,n2), $
RR : dblarr(n1), $
ISEL : intarr(n1)+1, $
UNIT : '', $
NPAR : 0, $
CHI2 : 0d }
s1 = CREATE_STRUCT('TRANS', {X : dblarr(n1,n3), Y : dblarr(n1,n3) })
strct = CREATE_STRUCT( strct, s1 )
endif else if n_elements(BAND) EQ 0 then begin
strct = { X : dblarr(n1), $
YD : dblarr(n1), $
ERR : dblarr(n1) , $
YM : dblarr(n1,n2), $
ISEL : intarr(n1)+1, $
UNIT : '', $
NPAR : 0, $
CHI2 : 0d }
endif
RETURN, strct
END
FUNCTION CHI2, y, model, npar, err=err
; returns the chi-square value of fit MODEL to data Y
;
; NPAR (I): nr of parameters in the fit
; ERR (I): error of each data point. Default is ERR=0.1*Y
ny = n_elements(y)
if n_elements(err) EQ 0 OR TOTAL(err) EQ 0 then begin
err = 0.1*y
print,'(W) CHI2: error missing, set to 10 %'
endif
ndof = ny - npar ; nr of degrees of freedom
chi = TOTAL( ((y-model)/err)^2 )
if ndof GT 0 then chi = chi / ndof
return, chi
END
FUNCTION FIL_CHI2, sf, ntype=ntype
; fills in the Chi-square fields of an input SMDAT structure
; (see format in GET_BAND_FLUX)
;
; SF (I): input SMDAT structure
; NTYPE (I): index of SED to be used
if ntype EQ 0 then begin
tt = SIZE(sf.sed.y)
ntype = tt(2)
endif
ntag = N_TAGS(sf)
stag = TAG_NAMES(sf)
itg = WHERE( STRPOS(stag,'I_') GE 0, ctg )
if ctg EQ 0 then begin
print,'(F) FIL_CHI2: no instrument to fill in'
return,0
endif
yy=0d & mm=0d & ee=0d & ift=0
for k = 0, ctg-1 do begin
i1 = WHERE( sf.(itg(k)).isel EQ 1 AND sf.(itg(k)).err GT 0, c1)
if c1 GT 0 then begin
sf.(itg(k)).chi2 = $
CHI2(sf.(itg(k)).yd(i1), sf.(itg(k)).ym(i1,ntype-1), sf.(itg(k)).npar, err=sf.(itg(k)).err(i1))
endif
yy = [ yy, sf.(itg(k)).yd ]
ee = [ ee, sf.(itg(k)).err]
mm = [ mm, sf.(itg(k)).ym(*,ntype-1) ]
ift = [ ift, sf.(itg(k)).isel ]
endfor
yy=yy(1:*) & ee=ee(1:*) & mm=mm(1:*) & ift=ift(1:*)
it = WHERE( ift EQ 1 AND ee GT 0, ct)
if ct GT 0 then sf.chi2 = CHI2( yy(it), mm(it), sf.npar, err=ee(it))
return, sf
END
FUNCTION ADD_INST, sf, tag, nsz, ntrans=ntrans
; fills in data for existing instruments in SMDAT structure
; or add a new instrument to SMDAT structure
;
; SF (I): SMDAT input structure
; TAG (I): name of instrument to add in
; NSZ (I): array(2) of sizes [nr of data pts, nr of grain types]
; NTRANS (I): nr of points for transmission (bad flux data)
stag = TAG_NAMES(sf)
ntag = N_TAGS(sf)
tag = STRUPCASE('I_'+STRTRIM(tag,2))
itg = WHERE( stag EQ tag, ctg )
if n_elements(nsz) EQ 2 then begin
nx = nsz(0) & ntype = nsz(1)
endif else begin
print,'(F) ADD_MOD: array of sizes must be 2D'
endelse
if ctg EQ 0 then begin
if n_elements(ntrans) EQ 0 then st = CREATE_STRUCT( tag, STR_INST(nx,n2=ntype) ) $
else st = CREATE_STRUCT( tag, STR_INST(nx,n2=ntype,n3=ntrans) )
sf = CREATE_STRUCT( sf, st )
stag = TAG_NAMES(sf)
itg = WHERE( stag EQ tag )
sf.(itg).isel = intarr(nx) + 1
endif else begin
print,'(F) ADD_INST: INST already exists in structure'
endelse
return,sf
END