dustem_brute_force_fit.pro
2.83 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
FUNCTION dustem_brute_force_fit,sed,table_name,filters,normalize=normalize,fact=fact,chi2=chi2,rchi2=rchi2,show_sed=show_sed
;stop
;sed must be normalized to NH=1e20 H/cm2
defsysv,'!dustem_grid',0,exist=exist
IF exist EQ 0 THEN BEGIN
dustem_define_grid,table_name,filters
ENDIF
;help,!dustem_grid,/str
;** Structure <33fa4168>, 6 tags, length=544, data length=540, refs=2:
; NSED LONG 8
; NFILTERS LONG 8
; NPARAMS LONG 3
; FILTERS STRING Array[8]
; PARAMS STRING Array[3]
; SEDS STRUCT -> <Anonymous> Array[8]
;help,!dustem_grid.seds,/str
;IDL> help,!dustem_grid.seds,/str
;** Structure <33f4de68>, 11 tags, length=44, data length=44, refs=3:
; P1 FLOAT 0.500000
; P2 FLOAT 0.00100000
; P3 FLOAT 0.00100000
; INIRCAM11 FLOAT 1.90915e-05
; INIRCAM16 FLOAT 2.40584e-05
; INIRCAM19 FLOAT 0.00566705
; INIRCAM21 FLOAT 0.00130508
; IMIRI2 FLOAT 0.0194735
; IMIRI3 FLOAT 0.0153818
; IMIRI6 FLOAT 0.0844443
; IMIRI11 FLOAT 0.0264565
seds=!dustem_grid.seds
Nseds=!dustem_grid.Nsed
Nfilters=!dustem_grid.Nfilters
Nparams=!dustem_grid.Nparams
chi2s=dblarr(Nseds)
grid_seds=fltarr(Nseds,Nfilters)
grid_param_values=fltarr(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
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 to compute chi2
;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,*])
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)
IF keyword_set(show_sed) THEN BEGIN
xtit='Wavelength [mic]'
ytit='SED (model=red,data=blue)'
cgplot,sed.wave,sed.stokesI,psym='Filled Circle',color='blue',/ylog,yrange=[0.01,10],/ysty,xtit=xtit,ytit=ytit,/xlog
cgoplot,sed.wave,grid_seds[ind[0],*],psym='Filled Square',color='red'
ENDIF
params=grid_param_values[ind[0],*]
fact=facts[ind[0]]
Nfreedom=Nseds-Nparams-1
rchi2=chi2/Nfreedom
RETURN,params
END