OV-GSO-DC / dustem-wrapper_idl

Blame view

src/idl/dustem_brute_force_fit.pro 2.83 KB
Edit Raw Normal View History
6ff7d98f   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
 
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=1e21 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