dustem_brute_force_fit.pro
5.59 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
FUNCTION dustem_brute_force_fit,sed $
,table_name $
,filters $
,normalize=normalize $
,params_hit=params_hit $
,params_uncertainties=params_uncertainties $
,params_min=params_min $
,params_max=params_max $
,fact=fact $
,chi2=chi2 $
,rchi2=rchi2 $
,best_sed=best_sed $
,best_grid_sed=best_grid_sed $
,show_sed=show_sed $
,nostop=nostop $
,reset=reset
;+
; NAME:
; dustem_brute_force_fit
; PURPOSE:
; does brut force fit from a pre-computed model grid
; CATEGORY:
; DustEM
; CALLING SEQUENCE:
; res=dustem_brute_force_fit(sed,table_name,filters[,/normalize][,fact=][,chi2=][,rchi2=][rchi2][best_sed=][,/show_sed][,/nostop][,/reset])
; INPUTS:
; sed = sed to be fit. Should contain all filters initialy provided to the routine
; table_name = name of the fits file to be used as a grid (see dustem_make_sed_table.pro for creation)
; filters = filters to be used in the grid table (only used if !dustem_grid does not exist or if /reset)
; OPTIONAL INPUT PARAMETERS:
; None
; OUTPUTS:
; res = best values of the free parameters. Not that linear parameters to the sed should be multiplied by fact.
; OPTIONAL OUTPUT PARAMETERS:
; fact = best scaling factor. Such that best_sed*fact is for NH=1.e20 H/cm2 (as in the grid). IF sed is already normalized to NH=1e20 H/cm2, no need to set /normalize
; chi2 = chi2 of the best fit
; rchi2 = reduced chi2 of the best fit
; best_grid_sed = best SED found in the grid
; best_sed = best fit SED matching input SED (such that best_sed=best_grid_sed*fact)
; params_hit= 1 if parameter hits the top, -1 if hits the bottom, 0 otherwise
; params_uncertainties = 90% confidence level uncertainty on parameters
; params_min = 90% lower limit of parameters
; params_max = 90% upper limit of parameters
; ACCEPTED KEY-WORDS:
; help = if set, print this help
; normalize = if set, the best match is done based on the shape of the SEDs only, instead of their absolute values. fact will be computed.
; show_sed = if set, plot the seds
; nostop = if set, do not stop
; COMMON BLOCKS:
; None
; SIDE EFFECTS:
; None
; RESTRICTIONS:
; The DustEMWrap IDL code must be installed
; PROCEDURE:
; Scans the grid SEDs to find the best fit (lowest chi2)
; EXAMPLES
;
; MODIFICATION HISTORY:
; Written by JPB Jan 2024
; 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_brute_force_fit'
params=0.
goto,the_end
ENDIF
;stop
;sed must be normalized to NH=1e20 H/cm2
;==== test for the existence of a loaded grid and initialize if needed.
defsysv,'!dustem_grid',0,exist=exist
IF exist EQ 0 or keyword_set(reset)THEN BEGIN
dustem_define_grid,table_name,filters
ENDIF
seds=!dustem_grid.seds
Nseds=!dustem_grid.Nsed
Nfilters=!dustem_grid.Nfilters
Nparams=!dustem_grid.Nparams
chi2s=dblarr(Nseds)
grid_seds=dblarr(Nseds,Nfilters)
total_seds=dblarr(Nseds)
grid_param_values=dblarr(Nseds,Nparams)
;=== compute sed grid Array
FOR i=0L,Nseds-1 DO BEGIN
FOR j=0L,Nfilters-1 DO BEGIN
grid_seds[i,j]=seds[i].(j+Nparams)
ENDFOR
total_seds[i]=seds[i].total
ENDFOR
;=== compute the parameter values Array
FOR i=0L,Nseds-1 DO BEGIN
FOR j=0L,Nparams-1 DO BEGIN
grid_param_values[i,j]=seds[i].(j)
ENDFOR
ENDFOR
;stop
;=== loop over grid lines to compute chi2s
;we should actually store the total of grid sed for each grid sed into the grid fits file (!)
;we could also store the total of observed seds in the observed sed save fits file (!)
facts=dblarr(Nseds)
facts[*]=1.d0
IF keyword_set(normalize) THEN BEGIN
FOR i=0L,Nseds-1 DO BEGIN
;facts[i]=total(sed.stokesI)/total(grid_seds[i,*])
facts[i]=total(sed.stokesI)/total_seds[i]
this_sed=sed.stokesI/facts[i]
this_var=sed.sigmaII/facts[i]^2
;chi2s[i]=total((sed.stokesI-grid_seds[i,*])^2/sed.sigmaII)
chi2s[i]=total((this_sed-grid_seds[i,*])^2/this_var)
ENDFOR
ENDIF ELSE BEGIN
FOR i=0L,Nseds-1 DO BEGIN
chi2s[i]=total((sed.stokesI-grid_seds[i,*])^2/sed.sigmaII)
ENDFOR
ENDELSE
;stop
chi2=min(chi2s)
ind=where(chi2s EQ chi2,count)
;==== should check bumping to table edges there
params_hit=intarr(Nparams)
FOR i=0L,Nparams-1 DO BEGIN
IF params[i] EQ minvalues[i] THEN params_hit[i]=-1
IF params[i] EQ maxvalues[i] THEN params_hit[i]=1
ENDFOR
params=grid_param_values[ind[0],*] ;these are the best fit parameters
fact=facts[ind[0]] ; This is the corresponding scaling factor
Nfreedom=Nseds-Nparams-1 ;This is the degree of freedom
rchi2=chi2/Nfreedom ;This is the reduced chi2
best_grid_sed=reform(grid_seds[ind[0],*]) ;This is the best SED fround in the grid.
best_sed=best_grid_sed*fact ;This is the best scaled grid SED matching the input SED.
;==== should compute parameter uncertainties there
conf_level=90./100.
dchi2=delta_chi2(Nfreedom,conf_level)
ind=where(chi2s LE chi2+dchi2,count)
pmins=fltarr(Nparams)
pmaxs=fltarr(Nparams)
FOR i=0L,Nparams-1 DO BEGIN
params_min[i]=min(params[ind,i])
params_max[i]=max(params[ind,i])
ENDFOR
params_uncertainties=params_max-params_min
IF keyword_set(show_sed) THEN BEGIN
xtit='Wavelength [mic]'
ytit='SED (model=red,data=blue)'
;should also plot uncertainties
minv=la_min([sed.stokesI,best_sed])
maxv=la_max([sed.stokesI,best_sed])
cgplot,sed.wave,sed.stokesI,psym='Filled Circle',color='blue',/ylog,/ysty,xtit=xtit,ytit=ytit,/xlog,yrange=[minv,maxv]
cgoplot,sed.wave,best_sed,psym='Filled Square',color='red'
if not keyword_set(nostop) THEN stop
ENDIF
the_end:
RETURN,params
END