Blame view

src/idl/dustem_mpfit_run.pro 16.7 KB
40597848   Ilyes Choubani   Plotting of numbe...
1
FUNCTION dustem_mpfit_run,x,p_min,niter=niter_completed,_EXTRA=_extra
427f1205   Jean-Michel Glorian   version 4.2 merged
2
3
4

;+
; NAME:
9ded3be9   Ilyes Choubani   Display update. F...
5
;    dustem_mpfit_run 
a2745344   Annie Hughes   updated inline do...
6
;
427f1205   Jean-Michel Glorian   version 4.2 merged
7
8
; PURPOSE:
;    function used to fit dustem model with dustem_mpfit_sed
a2745344   Annie Hughes   updated inline do...
9
;
427f1205   Jean-Michel Glorian   version 4.2 merged
10
11
; CATEGORY:
;    Dustem
a2745344   Annie Hughes   updated inline do...
12
;
427f1205   Jean-Michel Glorian   version 4.2 merged
13
; CALLING SEQUENCE:
36e8b879   Jean-Michel Glorian   update from deborah
14
;    res=dustem_mpfit_run_vg(x,params[,_EXTRA=extra][,cc_sed=cc_sed][,/help])
a2745344   Annie Hughes   updated inline do...
15
;
427f1205   Jean-Michel Glorian   version 4.2 merged
16
17
18
; INPUTS:
;    x         = necessary for interface but unused
;    params    = Dustem Model parameters normalized to initial values
a2745344   Annie Hughes   updated inline do...
19
;
427f1205   Jean-Michel Glorian   version 4.2 merged
20
21
22
; OPTIONAL INPUT PARAMETERS:
;    cc_sed    = SED color correction factors
;    _extra    = extra parameters for the plot routine
a2745344   Annie Hughes   updated inline do...
23
;
427f1205   Jean-Michel Glorian   version 4.2 merged
24
25
; OUTPUTS:
;    None
a2745344   Annie Hughes   updated inline do...
26
;
427f1205   Jean-Michel Glorian   version 4.2 merged
27
28
; OPTIONAL OUTPUT PARAMETERS:
;    None
a2745344   Annie Hughes   updated inline do...
29
;
427f1205   Jean-Michel Glorian   version 4.2 merged
30
31
; ACCEPTED KEY-WORDS:
;    help      = If set, print this help
a2745344   Annie Hughes   updated inline do...
32
;
427f1205   Jean-Michel Glorian   version 4.2 merged
33
34
; COMMON BLOCKS:
;    None
a2745344   Annie Hughes   updated inline do...
35
;
427f1205   Jean-Michel Glorian   version 4.2 merged
36
37
; SIDE EFFECTS:
;    SED and model are plotted if !dustem_show_plot
a2745344   Annie Hughes   updated inline do...
38
;
427f1205   Jean-Michel Glorian   version 4.2 merged
39
40
; RESTRICTIONS:
;    The dustem idl wrapper must be installed
a2745344   Annie Hughes   updated inline do...
41
;
427f1205   Jean-Michel Glorian   version 4.2 merged
42
43
; PROCEDURE:
;    None
a2745344   Annie Hughes   updated inline do...
44
;
427f1205   Jean-Michel Glorian   version 4.2 merged
45
; EXAMPLES
36e8b879   Jean-Michel Glorian   update from deborah
46
;    res=dustem_mpfit_run_vg(x,params[,_EXTRA=extra][,cc_sed=cc_sed][,/help])
a2745344   Annie Hughes   updated inline do...
47
48
;
; MODIFICATION HISTORY:
427f1205   Jean-Michel Glorian   version 4.2 merged
49
50
;    Adapted to extinction and polarization by Vincent Guillet (IAS)
;    Written by J.-Ph. Bernard
a2745344   Annie Hughes   updated inline do...
51
52
;    Evolution details on the DustEMWrap gitlab.
;    See http://dustemwrap.irap.omp.eu/ for FAQ and help.  
427f1205   Jean-Michel Glorian   version 4.2 merged
53
54
55
56
;-

;note that x is not used.

427f1205   Jean-Michel Glorian   version 4.2 merged
57
IF keyword_set(help) THEN BEGIN
9be94157   Jean-Philippe Bernard   improved
58
  doc_library,'dustem_mpfit_run'
427f1205   Jean-Michel Glorian   version 4.2 merged
59
60
61
  sed=0.
  goto,the_end
ENDIF
427f1205   Jean-Michel Glorian   version 4.2 merged
62
63

p_dim=p_min*(*(*!dustem_fit).param_init_values)
427f1205   Jean-Michel Glorian   version 4.2 merged
64

427f1205   Jean-Michel Glorian   version 4.2 merged
65
IF !run_ionfrac gt 0. THEN toto = dustem_create_ionfrac()
452c334e   Ilyes Choubani   Implementation Of...
66

427f1205   Jean-Michel Glorian   version 4.2 merged
67
;RUN THE MODEL AND GET THE SPECTRUM
45da29c6   Ilyes Choubani   updating mpfit_run
68

2a8e29d5   Ilyes Choubani   small corrections...
69
dustem_activate_plugins,p_dim/(*(*!dustem_fit).param_init_values),st=st ;st is an output...
36e8b879   Jean-Michel Glorian   update from deborah
70
st_model=dustem_read_all(!dustem_soft_dir)
427f1205   Jean-Michel Glorian   version 4.2 merged
71

427f1205   Jean-Michel Glorian   version 4.2 merged
72
73
chi2_sed      = 0  
chi2_ext      = 0  
452c334e   Ilyes Choubani   Implementation Of...
74
75
76
77
chi2_qsed      =0
chi2_used      =0 
chi2_qext      =0
chi2_uext      =0 
18e4331f   Ilyes Choubani   general update (f...
78

427f1205   Jean-Michel Glorian   version 4.2 merged
79
80
rchi2_sed     = 0
rchi2_ext     = 0
452c334e   Ilyes Choubani   Implementation Of...
81
82
83
84
rchi2_qsed     =0
rchi2_used     =0
rchi2_qext     =0
rchi2_uext     =0
18e4331f   Ilyes Choubani   general update (f...
85

427f1205   Jean-Michel Glorian   version 4.2 merged
86
87
wrchi2_sed    = 0
wrchi2_ext    = 0
452c334e   Ilyes Choubani   Implementation Of...
88
89
90
91
wrchi2_qsed    =0
wrich2_used    =0
wrchi2_qext    =0
wrich2_uext    =0
18e4331f   Ilyes Choubani   general update (f...
92

427f1205   Jean-Michel Glorian   version 4.2 merged
93
94
n_sed         = 0
n_ext         = 0
452c334e   Ilyes Choubani   Implementation Of...
95
96
97
98
n_qsed         =0
n_used         =0
n_qext         =0
n_uext         =0
36e8b879   Jean-Michel Glorian   update from deborah
99
100

win = 0
427f1205   Jean-Michel Glorian   version 4.2 merged
101

0068116a   Ilyes Choubani   General update + ...
102
tagnames=tag_names((*!dustem_data))
427f1205   Jean-Michel Glorian   version 4.2 merged
103
dustem_out = [0]
36e8b879   Jean-Michel Glorian   update from deborah
104
105
ind_data_tag=dustem_data_tag(count=count_data_tag)

b5ccb706   Jean-Philippe Bernard   improved to fit p...
106
FOR i_tag=0,count_data_tag-1 DO BEGIN
427f1205   Jean-Michel Glorian   version 4.2 merged
107

36e8b879   Jean-Michel Glorian   update from deborah
108
    CASE tagnames(ind_data_tag(i_tag)) OF
266ae799   Ilyes Choubani   General update
109
        
b74c4452   Ilyes Choubani   Treating extincti...
110
		'SED'	: BEGIN 
283845da   Ilyes Choubani   Small corrections...
111
	        dustem_sed = dustem_compute_sed(p_dim,st=st,SED_spec=SED_spec)
2c22f184   Jean-Philippe Bernard   fixed bug when so...
112
113
          ;help,*(*!dustem_fit).CURRENT_PARAM_VALUES
          ;stop
0068116a   Ilyes Choubani   General update + ...
114
	        chi2_sed   = total(((*(*!dustem_data).sed).values-dustem_sed)^2/(*(*!dustem_data).sed).sigma^2)
2c22f184   Jean-Philippe Bernard   fixed bug when so...
115
			    n_sed = n_elements((*(*!dustem_data).sed).values)
b5ccb706   Jean-Philippe Bernard   improved to fit p...
116
	       	rchi2_sed  = chi2_sed / n_sed
d0f68bb6   Jean-Philippe Bernard   replaced system v...
117
	       	wrchi2_sed  = rchi2_sed * (*!fit_rchi2_weight).sed
b5ccb706   Jean-Philippe Bernard   improved to fit p...
118
	       	message,"chi2 SED = "+strtrim(chi2_sed,2)+"  rchi2 SED = "+strtrim(rchi2_sed,2),/continue;," weighted rchi2 SED=",wrchi2_sed
0068116a   Ilyes Choubani   General update + ...
119
120
121
122
123
124
125
126
			;print,'wl ',(*(*!dustem_data).sed).wav
			;print,"Delta SED per wl: ",(*(*!dustem_data).sed).values-dustem_sed
			;print,"error per wl: ",(*(*!dustem_data).sed).sigma
			;print,"chi2 per wl: ",((*(*!dustem_data).sed).values-dustem_sed)^2/(*(*!dustem_data).sed).sigma^2
	        alpha = total(((*(*!dustem_data).sed).values*dustem_sed)/(*(*!dustem_data).sed).sigma^2) $
	              		/ total(((*(*!dustem_data).sed).values)^2/(*(*!dustem_data).sed).sigma^2)
	        beta = total(((*(*!dustem_data).sed).values*dustem_sed)/(*(*!dustem_data).sed).sigma^2) $
	             		/ total(dustem_sed^2/(*(*!dustem_data).sed).sigma^2)
427f1205   Jean-Michel Glorian   version 4.2 merged
127
		
b5ccb706   Jean-Philippe Bernard   improved to fit p...
128
	        print,"Best SED chi2 obtained if f_HI = ",alpha, " or if we multiply model by ",beta
b74c4452   Ilyes Choubani   Treating extincti...
129
130
131
132
133
             
            ;These two lines were put here so that the chi2 is updated everytime we plot
            ;Moving it to the bottom so that all chi2 are taken into account  
;  	        (*!dustem_fit).chi2=chi2_sed
;  	        (*!dustem_fit).rchi2=rchi2_sed    
fcb6eade   Ilyes Choubani   general update - ...
134
	        
5a2643a4   Ilyes Choubani   cleaning/improvin...
135
	        dustem_out = [ dustem_out, dustem_sed ]
266ae799   Ilyes Choubani   General update
136
            
36e8b879   Jean-Michel Glorian   update from deborah
137
138
		END

b74c4452   Ilyes Choubani   Treating extincti...
139
		'EXT' 	: BEGIN     		
283845da   Ilyes Choubani   Small corrections...
140
            dustem_ext = dustem_compute_ext(p_dim,st=st,EXT_spec=EXT_spec)   
0068116a   Ilyes Choubani   General update + ...
141
142
	        chi2_ext   = total(((*(*!dustem_data).ext).values-dustem_ext)^2/(*(*!dustem_data).ext).sigma^2)
			n_ext = n_elements((*(*!dustem_data).ext).values)
b5ccb706   Jean-Philippe Bernard   improved to fit p...
143
			rchi2_ext  = chi2_ext / n_ext
d0f68bb6   Jean-Philippe Bernard   replaced system v...
144
	        wrchi2_ext  = rchi2_ext * (*!fit_rchi2_weight).ext
b74c4452   Ilyes Choubani   Treating extincti...
145
	        print,"chi2 EXT = ",chi2_ext,"  rchi2 EXT = ",rchi2_ext;," weighted rchi2 EXT=",wrchi2_ext
0068116a   Ilyes Choubani   General update + ...
146
        	;print,"chi2 per wl: ",((*(*!dustem_data).ext).values-dustem_ext)^2/(*(*!dustem_data).ext).sigma^2
b5ccb706   Jean-Philippe Bernard   improved to fit p...
147
148
149
150
151

		 	dustem_out = [ dustem_out, dustem_ext ]

		 	END

b74c4452   Ilyes Choubani   Treating extincti...
152
		'POLEXT': BEGIN 
b5ccb706   Jean-Philippe Bernard   improved to fit p...
153

759a527d   Ilyes Choubani   general update
154
			
b74c4452   Ilyes Choubani   Treating extincti...
155
			;Haven't investigated !run_circ=1 yet... 
759a527d   Ilyes Choubani   general update
156
			;---------------the procedure within this block has not been updated for a vert long time
b5ccb706   Jean-Philippe Bernard   improved to fit p...
157
158
159
160
161
			IF !dustem_show_plot NE 0 and !run_circ then begin
	                   win+=1
	                 
				dustem_plot_circ,st,win=win,xr=[0.1,5],aligned=st_model.align.grains.aligned,yr=[-0.03,0.03]
			endif
759a527d   Ilyes Choubani   general update
162
163
164
			;--------------------------------------------------------------------------------------
			
			
b5ccb706   Jean-Philippe Bernard   improved to fit p...
165
			IF !run_lin then begin
3f93c873   Ilyes Choubani   correcting "compu...
166
                     
283845da   Ilyes Choubani   Small corrections...
167
                    dustem_polext = dustem_compute_polext(p_dim,st=st,POLEXT_spec=POLEXT_spec,SPEXT_spec=SPEXT_spec,dustem_fpolext=dustem_fpolext,dustem_ext=dustem_ext)
36e8b879   Jean-Michel Glorian   update from deborah
168

b5ccb706   Jean-Philippe Bernard   improved to fit p...
169
			ENDIF
b5ccb706   Jean-Philippe Bernard   improved to fit p...
170
			END
427f1205   Jean-Michel Glorian   version 4.2 merged
171

b74c4452   Ilyes Choubani   Treating extincti...
172
        'QSED': BEGIN 
45da29c6   Ilyes Choubani   updating mpfit_run
173
        		    
283845da   Ilyes Choubani   Small corrections...
174
		        stokes = dustem_compute_stokes(p_dim,st=st,Q_spec=Q_spec,U_spec=U_spec,PSI_spec=PSI_spec,dustem_psi_em=dustem_psi_em) 
3f93c873   Ilyes Choubani   correcting "compu...
175
    		    dustem_qsed = stokes[0]  
45da29c6   Ilyes Choubani   updating mpfit_run
176
177
		        ;taking care of the color correction for spass ; SHOULD NOT LEAVE THIS ON
		        ;dustem_qsed(-1)=dustem_qsed(-2)
0068116a   Ilyes Choubani   General update + ...
178
; 		        testspass1 = ((*(*!dustem_data).qsed).filt_names)(-1) eq 'SPASS1' 
45da29c6   Ilyes Choubani   updating mpfit_run
179
180
;  			    if testspass1 then dustem_qsed(-1)=dustem_qsed(-2)
;  			    
0068116a   Ilyes Choubani   General update + ...
181
182
			    chi2_qsed = total(((*(*!dustem_data).qsed).values-dustem_qsed)^2/(*(*!dustem_data).qsed).sigma^2) 
			    n_qsed = n_elements((*(*!dustem_data).qsed).values)
45da29c6   Ilyes Choubani   updating mpfit_run
183
			    rchi2_qsed  = chi2_qsed / n_qsed
d0f68bb6   Jean-Philippe Bernard   replaced system v...
184
	        	wrchi2_qsed  = rchi2_qsed * (*!fit_rchi2_weight).qsed 
5a2643a4   Ilyes Choubani   cleaning/improvin...
185
	            dustem_out = [ dustem_out, dustem_qsed ] ; classic command
45da29c6   Ilyes Choubani   updating mpfit_run
186
	            message,"chi2 QSED = "+strtrim(chi2_qsed,2)+"  rchi2 QSED = "+strtrim(rchi2_qsed,2),/continue;," weighted rchi2 POLSED=",wrchi2_polsed
0068116a   Ilyes Choubani   General update + ...
187
	        	print,"chi2 per wl: ",((*(*!dustem_data).qsed).values-dustem_qsed)^2/(*(*!dustem_data).qsed).sigma^2
45da29c6   Ilyes Choubani   updating mpfit_run
188
189
		        
		    END   
1355825c   Ilyes Choubani   General update
190

b74c4452   Ilyes Choubani   Treating extincti...
191
        'USED': BEGIN 
3f93c873   Ilyes Choubani   correcting "compu...
192
    		    dustem_used = stokes[1] 
0068116a   Ilyes Choubani   General update + ...
193
194
 			    chi2_used  = total(((*(*!dustem_data).used).values-dustem_used)^2/(*(*!dustem_data).used).sigma^2) 
 			    n_used = n_elements((*(*!dustem_data).used).values)
45da29c6   Ilyes Choubani   updating mpfit_run
195
 			    rchi2_used  = chi2_used / n_used
d0f68bb6   Jean-Philippe Bernard   replaced system v...
196
  	        	wrchi2_used  = rchi2_used * (*!fit_rchi2_weight).used
5a2643a4   Ilyes Choubani   cleaning/improvin...
197
 	            dustem_out = [ dustem_out, dustem_used] ; classic command
45da29c6   Ilyes Choubani   updating mpfit_run
198
 	            message,"chi2 USED = "+strtrim(chi2_used,2)+"  rchi2 USED = "+strtrim(rchi2_used,2),/continue;," weighted rchi2 POLSED=",wrchi2_polsed
0068116a   Ilyes Choubani   General update + ...
199
  	        	print,"chi2 per wl: ",((*(*!dustem_data).used).values-dustem_used)^2/(*(*!dustem_data).used).sigma^2
45da29c6   Ilyes Choubani   updating mpfit_run
200
201
202
                
  		END
  		
b74c4452   Ilyes Choubani   Treating extincti...
203
        'POLSED': BEGIN
45da29c6   Ilyes Choubani   updating mpfit_run
204
205
 			
             	  IF !run_pol and !run_lin THEN BEGIN
3f93c873   Ilyes Choubani   correcting "compu...
206
207
                        ;the keyword for dustem_sed here is used so that the function uses it to compute the polarization fraction instead of
                        ;running dustem_compute_sed again in the procedure. Which it used to do.  
283845da   Ilyes Choubani   Small corrections...
208
            	       	dustem_polsed = dustem_compute_polsed(p_dim,st=st,P_spec=P_spec,SP_spec=SP_spec,dustem_polfrac=dustem_polfrac,dustem_sed=dustem_sed)
bffa7018   Ilyes Choubani   update
209
            	      	
45da29c6   Ilyes Choubani   updating mpfit_run
210
211
         		  ENDIF
		END
67bd858a   Ilyes Choubani   Replication of th...
212
        	
452c334e   Ilyes Choubani   Implementation Of...
213
 		
b74c4452   Ilyes Choubani   Treating extincti...
214
  		'QEXT': BEGIN ;UPDATED  
52a3c8dc   Jean-Philippe Bernard   mise a jour
215
                ;stop
283845da   Ilyes Choubani   Small corrections...
216
                Stokext = dustem_compute_stokext(p_dim,st=st,QEXT_spec=QEXT_spec,UEXT_spec=UEXT_spec,PSIEXT_spec=PSIEXT_spec,dustem_psi_ext=dustem_psi_ext) 
3f93c873   Ilyes Choubani   correcting "compu...
217
                dustem_qext = stokext[0]
0068116a   Ilyes Choubani   General update + ...
218
219
                chi2_qext   =total(((*(*!dustem_data).qext).values-dustem_qext)^2/(*(*!dustem_data).qext).sigma^2) 
   			    n_qext = n_elements((*(*!dustem_data).qext).values)
452c334e   Ilyes Choubani   Implementation Of...
220
   			    rchi2_qext  = chi2_qext / n_qext
d0f68bb6   Jean-Philippe Bernard   replaced system v...
221
  	        	wrchi2_qext  = rchi2_qext * (*!fit_rchi2_weight).qext
452c334e   Ilyes Choubani   Implementation Of...
222
223
   	            dustem_out = [ dustem_out, dustem_qext] ; classic command
   	            message,"chi2 QEXT = "+strtrim(chi2_qext,2)+"  rchi2 QEXT = "+strtrim(rchi2_qext,2),/continue;," weighted rchi2 POLSED=",wrchi2_polsed
0068116a   Ilyes Choubani   General update + ...
224
  	        	print,"chi2 per wl: ",((*(*!dustem_data).qext).values-dustem_qext)^2/(*(*!dustem_data).qext).sigma^2 
452c334e   Ilyes Choubani   Implementation Of...
225
226
227
     	           
  		END
;  		
b74c4452   Ilyes Choubani   Treating extincti...
228
  		'UEXT': BEGIN ;UPDATED
3f93c873   Ilyes Choubani   correcting "compu...
229
          		dustem_uext = stokext[1]
0068116a   Ilyes Choubani   General update + ...
230
231
 			    chi2_uext   =total(((*(*!dustem_data).uext).values-dustem_uext)^2/(*(*!dustem_data).uext).sigma^2) 
   			    n_uext = n_elements((*(*!dustem_data).uext).values)
452c334e   Ilyes Choubani   Implementation Of...
232
   			    rchi2_uext  = chi2_uext / n_uext
d0f68bb6   Jean-Philippe Bernard   replaced system v...
233
  	        	wrchi2_uext  = rchi2_uext * (*!fit_rchi2_weight).uext
452c334e   Ilyes Choubani   Implementation Of...
234
235
   	            dustem_out = [ dustem_out, dustem_uext] ; classic command
   	            message,"chi2 UEXT = "+strtrim(chi2_uext,2)+"  rchi2 UEXT = "+strtrim(rchi2_uext,2),/continue;," weighted rchi2 POLSED=",wrchi2_polsed
0068116a   Ilyes Choubani   General update + ...
236
  	        	print,"chi2 per wl: ",((*(*!dustem_data).uext).values-dustem_uext)^2/(*(*!dustem_data).uext).sigma^2
b74c4452   Ilyes Choubani   Treating extincti...
237
            
452c334e   Ilyes Choubani   Implementation Of...
238
239
 	           
  		END
67bd858a   Ilyes Choubani   Replication of th...
240
      	
b74c4452   Ilyes Choubani   Treating extincti...
241
242
	    ELSE : ;Do nothing
	    
427f1205   Jean-Michel Glorian   version 4.2 merged
243
244
    ENDCASE

b5ccb706   Jean-Philippe Bernard   improved to fit p...
245
ENDFOR
b74c4452   Ilyes Choubani   Treating extincti...
246

9eb53891   Jean-Philippe Bernard   modified call to ...
247
;=== computing chi^2
36e8b879   Jean-Michel Glorian   update from deborah
248
249
250
251
DOF = n_elements(p_min)
total_chi2 = 0
total_rchi2 = 0
total_n = 0
d0f68bb6   Jean-Philippe Bernard   replaced system v...
252
if (*!fit_rchi2_weight).ext ne 0 then begin
36e8b879   Jean-Michel Glorian   update from deborah
253
254
255
256
	total_n += n_ext
	total_chi2 +=  chi2_ext
	total_rchi2 += rchi2_ext
endif
d0f68bb6   Jean-Philippe Bernard   replaced system v...
257
if (*!fit_rchi2_weight).sed ne 0 then begin
36e8b879   Jean-Michel Glorian   update from deborah
258
259
260
261
262
	total_n += n_sed
	total_chi2 += chi2_sed
	total_rchi2 += rchi2_sed
endif
if !run_pol then begin
9ded3be9   Ilyes Choubani   Display update. F...
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
    
    if tag_exist(*!fit_rchi2_weight,'qsed') then begin ;takes into account the activation$
    ;of the !run_pol keyword without the user knowing it. 
    ;ie: the use of a model with polarization and without polarization data or the user wanting to fit polarization.
        
        if (*!fit_rchi2_weight).qsed ne 0 then begin
        	total_n += n_qsed
        	total_chi2 += chi2_qsed
        	total_rchi2 += rchi2_qsed
        endif
        if (*!fit_rchi2_weight).used ne 0 then begin
        	total_n += n_used
        	total_chi2 += chi2_used
        	total_rchi2 += rchi2_used
        endif
        
1355825c   Ilyes Choubani   General update
279
    endif
9ded3be9   Ilyes Choubani   Display update. F...
280
281
282
283
284
285
286
287
288
289
290
291
292
293
    
    if tag_exist(*!fit_rchi2_weight,'qext') then begin
        
        if (*!fit_rchi2_weight).qext ne 0 then begin
        	total_n += n_qext
        	total_chi2 += chi2_qext
        	total_rchi2 += rchi2_qext
        endif
        if (*!fit_rchi2_weight).uext ne 0 then begin
        	total_n += n_uext
        	total_chi2 += chi2_uext
        	total_rchi2 += rchi2_uext
        endif
        
1355825c   Ilyes Choubani   General update
294
    endif
452c334e   Ilyes Choubani   Implementation Of...
295
296
endif

d0f68bb6   Jean-Philippe Bernard   replaced system v...
297
298
total_wchi2 = chi2_ext * (*!fit_rchi2_weight).ext $
            + chi2_sed * (*!fit_rchi2_weight).sed 
36e8b879   Jean-Michel Glorian   update from deborah
299
if !run_pol then begin
9ded3be9   Ilyes Choubani   Display update. F...
300
301
302
303
304
305
306
307
308
309
310
311
312
313

    if tag_exist(*!fit_rchi2_weight,'qsed') THEN total_wchi2+=chi2_qsed * (*!fit_rchi2_weight).qsed $
    +chi2_used * (*!fit_rchi2_weight).used
    
    if tag_exist(*!fit_rchi2_weight,'qext') THEN total_wchi2+=chi2_qext * (*!fit_rchi2_weight).qext $
    +chi2_used * (*!fit_rchi2_weight).used

; 	total_wchi2 	+=  $;chi2_polext  * !fit_rchi2_weight.polext $
; ;              		 + chi2_polsed  * !fit_rchi2_weight.polsed $
; ;              		 + chi2_polfrac * !fit_rchi2_weight.polfrac $
;              		 + chi2_qsed * (*!fit_rchi2_weight).qsed $
;              		 + chi2_used * (*!fit_rchi2_weight).used $
;              		 + chi2_qext * (*!fit_rchi2_weight).qext $
;              		 + chi2_uext * (*!fit_rchi2_weight).uext  
36e8b879   Jean-Michel Glorian   update from deborah
314
315
endif

b74c4452   Ilyes Choubani   Treating extincti...
316
317
318
;Using this expression for now for the weighting. 
(*!dustem_fit).chi2 = total_wchi2
(*!dustem_fit).rchi2 = total_wchi2/(total_n-DOF)
9ded3be9   Ilyes Choubani   Display update. F...
319

b74c4452   Ilyes Choubani   Treating extincti...
320

9ccf7615   Jean-Philippe Bernard   modified to fit u...
321
message,'======================================================================',/continue
36e8b879   Jean-Michel Glorian   update from deborah
322
323
324
print,"Total weighted chi2 driving mpfitfun = ", total_wchi2
print,"Nb of data points = ", total_n
print,"DOF = ", DOF
18e4331f   Ilyes Choubani   general update (f...
325
326
print,"Reduced chi2 EXT/SED/POLEXT/POLSED/POLFRAC",rchi2_ext,rchi2_sed;,rchi2_polext;,rchi2_polsed,rchi2_polfrac
print,"Weighted reduced chi2 EXT/SED/POLEXT/POLSED/POLFRAC",wrchi2_ext,wrchi2_sed;,wrchi2_polext;,wrchi2_polsed,wrchi2_polfrac
d0f68bb6   Jean-Philippe Bernard   replaced system v...
327
328
if !run_pol then 	print,"Weights",(*!fit_rchi2_weight).ext,(*!fit_rchi2_weight).sed $;,!fit_rchi2_weight.polext $;,!fit_rchi2_weight.polsed,!fit_rchi2_weight.polfrac $ 
else 			print,"Weights",(*!fit_rchi2_weight).ext,(*!fit_rchi2_weight).sed
36e8b879   Jean-Michel Glorian   update from deborah
329
330
print,"Mean weighted reduced chi2 = ", total_wchi2/(total_n-DOF)
print,"Unweighted reduced Total chi2 = ", total_chi2/(total_n-DOF)
b74c4452   Ilyes Choubani   Treating extincti...
331
print,"Unweighted mean reduced chi2 = ", total_rchi2/4*total_n/(total_n-DOF);total_rchi2/4*total_n/(total_n-DOF) original ; I do not understand this line. Maybe redo some computation.
9ccf7615   Jean-Philippe Bernard   modified to fit u...
332
message,'======================================================================',/continue
386ad2de   Jean-Philippe Bernard   cant remeber what...
333
334
;the following is needed under fawlty (not IDL) which message,/reset does not work like in IDL
!Error_State.msg=''
9eb53891   Jean-Philippe Bernard   modified call to ...
335

dc84fd94   Ilyes Choubani   Small corrections...
336
337
338
339
If !dustem_show_plot THEN BEGIN
    ;=== plotting the fit
    if !dustem_iter.act NE !dustem_iter.prv or !dustem_iter.prv EQ 1 then begin
        IF !dustem_noobj EQ 0 THEN BEGIN
478a4d49   Ilyes Choubani   Correcting some p...
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
          dustemwrap_plot,  p_dim, $
                            st = st, $
                            dustem_sed = dustem_sed, $
                            SED_spec = SED_spec, $
                            dustem_qsed = dustem_qsed, $
                            Q_spec = Q_spec, $
                            dustem_used = dustem_used, $
                            U_spec = U_spec, $
                            dustem_polsed = dustem_polsed,   $
                            P_spec = P_spec, $
                            dustem_polfrac = dustem_polfrac, $
                            SP_spec = SP_spec, $
                            dustem_psi_em = dustem_psi_em, $
                            PSI_spec = PSI_spec, $
                            dustem_ext = dustem_ext, $
                            EXT_spec = EXT_spec, $
                            dustem_qext = dustem_qext,       $
                            QEXT_spec = QEXT_spec, $
                            dustem_uext = dustem_uext, $
                            UEXT_spec = UEXT_spec, $
                            dustem_polext = dustem_polext, $
                            POLEXT_spec = POLEXT_spec, $
                            dustem_fpolext = dustem_fpolext, $
                            SPEXT_spec = SPEXT_spec,         $
                            dustem_psi_ext = dustem_psi_ext, $
                            PSIEXT_spec = PSIEXT_spec,       $
                             _extra=_extra
dc84fd94   Ilyes Choubani   Small corrections...
367
        ENDIF ELSE BEGIN
478a4d49   Ilyes Choubani   Correcting some p...
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
          dustemwrap_plot_noobj, p_dim,$
                                 st = st,$
                                 dustem_sed = dustem_sed, $
                                 SED_spec = SED_spec, $
                                 dustem_qsed = dustem_qsed, $
                                 Q_spec = Q_spec, $
                                 dustem_used = dustem_used, $
                                 U_spec = U_spec, $
                                 dustem_polsed = dustem_polsed, $ 
                                 P_spec = P_spec, $
                                 dustem_polfrac = dustem_polfrac, $
                                 SP_spec = SP_spec, $
                                 dustem_psi_em = dustem_psi_em, $
                                 PSI_spec = PSI_spec, $
                                 dustem_ext = dustem_ext, $
                                 EXT_spec = EXT_spec, $
                                 dustem_qext = dustem_qext, $
                                 QEXT_spec = QEXT_spec, $
                                 dustem_uext = dustem_uext, $
                                 UEXT_spec = UEXT_spec, $
                                 dustem_polext = dustem_polext, $
                                 POLEXT_spec = POLEXT_spec, $
                                 dustem_fpolext = dustem_fpolext, $
                                 SPEXT_spec = SPEXT_spec, $
                                 dustem_psi_ext = dustem_psi_ext, $
                                 PSIEXT_spec = PSIEXT_spec, $
                                 _extra = _extra
dc84fd94   Ilyes Choubani   Small corrections...
395
396
397
398
399
400
401
        ENDELSE
        if !dustem_iter.act EQ 1 then begin
            !dustem_iter.prv=0
            !dustem_iter.act=0 
        endif else !dustem_iter.prv=!dustem_iter.act
    endif
ENDIF
36e8b879   Jean-Michel Glorian   update from deborah
402
403
404
405
406
407
408
409
410

;print,"Total reduced chi2 = ", (wrchi2_sed+wrchi2_ext+wrchi2_polext+wrchi2_polsed+wrchi2_polfrac)/(n_sed+n_ext+n_polext+n_polsed+n_polfrac-DOF)
; Reduced chi2 is like if there were the same nb of data points for each curve
;print,"Mean reduced chi2 = ", (wrchi2_sed+wrchi2_ext+wrchi2_polext+wrchi2_polsed+wrchi2_polfrac) * total_n/(total_n-DOF)

; Plot size distribution
win+=1
xr=[0.3,5e3]
yr=[0.1,100]
9ccf7615   Jean-Philippe Bernard   modified to fit u...
411
;IF !dustem_show_plot NE 0 THEN dustem_plot_sdist,xr=xr,yr=yr,ps=filename,win=win;,/donotclose
427f1205   Jean-Michel Glorian   version 4.2 merged
412
413
414
415

; Suppress the artificial [0] at the beginning 
dustem_out = dustem_out[1:*]

9ccf7615   Jean-Philippe Bernard   modified to fit u...
416
417
418
419
420
421
422
;IF defined(extra) THEN BEGIN
;  IF keyword_set(extra.plot_it) THEN BEGIN
;;    message,'params='+strtrim(p_min,2),/info,/continue
;;    stop
;    print,'dustem_mpfit_run_vg: Current parameters values:',p_dim
;  ENDIF
;ENDIF
427f1205   Jean-Michel Glorian   version 4.2 merged
423
424
425
426
427
428

the_end:

RETURN,dustem_out

END