dustem_brute_force_fit.pro
4.51 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
FUNCTION dustem_brute_force_fit,sed,table_name,filters,normalize=normalize,fact=fact,chi2=chi2,rchi2=rchi2,show_sed=show_sed
;+
; 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=normalize,fact=fact,chi2=chi2,rchi2=rchi2,show_sed=show_sed)
; INPUTS:
; None
; OPTIONAL INPUT PARAMETERS:
; key = input parameter number
; First parameter: E(B-V) for extinction applied to the template, as used in Phangs
; Following 13*6 parameters are weights by which to multiply the MILES templates [Msun/pc^2]
; val = corresponding input parameter value
; OUTPUTS:
; stellar_brightness = stellar spectrum spectrum resampled on dustem wavelengths [MJy/sr]
; OPTIONAL OUTPUT PARAMETERS:
; scope = scope of the plugin
; paramdefault = default values of parameters
; paramtag = plugin parameter names as strings
; ACCEPTED KEY-WORDS:
; help = if set, print this help
; method = SSP used. can be EMILES or CB19. default='EMILES'
; COMMON BLOCKS:
; None
; SIDE EFFECTS:
; None
; RESTRICTIONS:
; The DustEMWrap IDL code must be installed
; PROCEDURE:
; This is a DustEMWrap plugin for phangs
; It differs from dustem_plugin_emiles_stellar_continuum.pro only by the way the extinction is applied to the output spectrum
; EXAMPLES
; dustem_init
; vec=dustem_plugin_phangs_stellar_continuum(scope=scope)
; MODIFICATION HISTORY:
; Written by JPB June 2023
; 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_plugin_phangs_stellar_continuum'
output=0.
goto,the_end
ENDIF
;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)
params=grid_param_values[ind[0],*]
fact=facts[ind[0]]
Nfreedom=Nseds-Nparams-1
rchi2=chi2/Nfreedom
IF keyword_set(show_sed) THEN BEGIN
xtit='Wavelength [mic]'
ytit='SED (model=red,data=blue)'
;should also plot uncertainties
cgplot,sed.wave,sed.stokesI,psym='Filled Circle',color='blue',/ylog,/ysty,xtit=xtit,ytit=ytit,/xlog
cgoplot,sed.wave,grid_seds[ind[0],*]*fact,psym='Filled Square',color='red'
ENDIF
RETURN,params
END