dustem_mpfit_data.pro
5.07 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
FUNCTION dustem_mpfit_data,tol=tol,Nitermax=Nitermax,gtol=gtol,function_name=function_name,xtol=xtol,_extra=_extra
;+
; NAME:
; dustem_mpfit_data
;
; PURPOSE:
; call MPFIT to fit observational datas using ponderated reduced chi2
; between extinction, emission, and polarisation as defined in !fit_rchi2_weight
;
; CATEGORY:
; Dustem
;
; CALLING SEQUENCE:
; res=dustem_mpfit_data(tol=tol,Nitermax=Nitermax,gtol=gtol)
;
; INPUTS:
; None
;
; OPTIONAL INPUT PARAMETERS:
;
; OUTPUTS:
; None
;
; OPTIONAL OUTPUT PARAMETERS:
;
; ACCEPTED KEY-WORDS:
; help = If set, print this help
;
; COMMON BLOCKS:
; None
;
; SIDE EFFECTS:
; None
;
; RESTRICTIONS:
; The DustEM fortran code must be installed
; The DustEMWrap IDL code must be installed
;
; PROCEDURE:
; None
;
; EXAMPLES
;
; MODIFICATION HISTORY:
; Written by V. Guillet (2012)
; 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_mpfit_data'
goto,the_end
ENDIF
disp_str='====================================================='
IF not keyword_set(function_name) THEN BEGIN
func_name='dustem_mpfit_run'
ENDIF ELSE BEGIN
func_name=function_name
ENDELSE
;functargs={plot_it:1}
use_tol=1.e-14 ;Default tolerance value
use_Nitermax=3 ;Default maximum number of iterations
IF keyword_set(tol) THEN use_tol=tol
IF keyword_set(Nitermax) THEN use_Nitermax=Nitermax
; Build array concatenating observational data
if !run_pol then begin
ind_data_tag=dustem_data_tag(count=count_data_tag)
tagnames=tag_names(*!dustem_data)
testpol = tagnames EQ 'POLSED' or tagnames EQ 'POLEXT' or tagnames EQ 'POLFRAC' or tagnames EQ 'PSI_EM' or tagnames EQ 'PSI_EXT' or tagnames EQ 'FPOLEXT'
data_ind=where(~testpol,ctestpol)
if ctestpol ne 0 then begin
match,ind_data_tag,data_ind,ind1,ind2 ;not using match2 because when matching indices, if the first is 0 it doesn't count which is wrong.
;we also only need the points that match.
ind_data_tag = ind_data_tag[ind1]
count_data_tag = n_elements(ind_data_tag)
endif
endif else ind_data_tag=dustem_data_tag(count=count_data_tag)
fields=['wav','values','sigma']
FOR j=0,n_elements(fields)-1 DO BEGIN
instruc = fields(j)+' = ([0' ; we add 0 first to allow for ',' after it
for i=0, count_data_tag-1 do begin $ ; Add each data set
; Adapt sigma to the reduced chi2 ponderation defined in !fit_rchi2_weight between sed, ext, polext and polsed
; Total_reduced_chi2_driving_mpfit = sum_i ( chi2_i * (weight_i/n_i) )
; n : the total number of observational variables in the data
; n_i : the number of obs. variable in data i
; w_i : the weight for data i from !fit_rchi2_weight for the ponderation of reduced chi2
; ex : w_sed = 0.8, w_ext = 0.2 , i.e. rchi2_sed counts 4 times as much as rchi2_ext
; sed : 314 pts, ext : 52 pts, i.e. 6 times much more constraints in sed than in ext
; sigma_sed *= sqrt(314./366/0.8) = 1.03
; sigma_ext *= sqrt(52./366./0.2) = 0.84
; By increasing sigma by a factor w, we diminish the weight of the corresponding point by a factor w^2
if count_data_tag gt 1 and j eq 2 then begin ; sigma ponderation = sqrt(ni/wi)
if (*!fit_rchi2_weight).(i) gt 0 then begin
weight = sqrt(1./(*!fit_rchi2_weight).(i) )
endif else begin
weight = 1.d99
endelse
endif else begin
weight = 1.
endelse
;print,"WEIGHT",i,j,weight,n_elements((*!dustem_data.(i)).sigma),"/",n_elements(wav)
instruc += ', ((*(*!dustem_data).(' + strtrim(ind_data_tag(i),2) + ')).'+fields(j)+')*'+strtrim(weight,2)
;print,"INSTRUC",instruc
endfor
instruc += '])(1:*)' ; from 1 to exclude first 0
;print,'-> INSTRUC: ',instruc
toto=execute(instruc) ; Variables wav, values & sigma are created
ENDFOR
;stop
IF tag_exist(_extra,'show_plot') THEN !dustem_show_plot = _extra.show_plot
nocatch=!dustem_nocatch
p_min = dustem_mpfitfun(func_name,wav,values,sigma, $
parinfo=(*!dustem_parinfo), perror=perror, yfit=yfit, $
bestnorm=bestnorm, dof=dof,quiet=quiet, status=status, xtol=xtol, niter=niter_completed, ftol=use_tol, $
maxiter=use_Nitermax,gtol=gtol,errmsg=errmsg,functargs=_extra,nocatch=nocatch)
;stop
;!dustem_iter.act=niter_completed
message,errmsg,/info
message,'mpfitfun Status='+strtrim(status,2),/info
p_dim=p_min*(*(*!dustem_fit).param_init_values)
!dustem_end = 1
;PRINT RESULTS
print,''
print,''
print,'<o)))>< <o)))>< <o)))>< <o)))>< <o)))>< <o)))>< <o)))>< <o)))><'
print,'<o)))>< <o)))>< <o)))>< <o)))>< <o)))>< <o)))>< <o)))>< <o)))><'
print,''
print,'FINISHED DUSTEMWRAP FITTING...'
print,''
print,'The final values of the free parameters (pd) are : '
print, p_dim
print,''
print,'<o)))>< <o)))>< <o)))>< <o)))>< <o)))>< <o)))>< <o)))>< <o)))><'
print,'<o)))>< <o)))>< <o)))>< <o)))>< <o)))>< <o)))>< <o)))>< <o)))><'
print,''
print,''
IF !dustem_noobj THEN BEGIN
dustemwrap_plot_noobj,p_dim,stp,_extra=_extra
ENDIF ELSE BEGIN
dustemwrap_plot,p_dim,stp,_extra=_extra
ENDELSE
the_end:
RETURN, p_dim
END