Blame view

src/idl/dustem_mpfit_run.pro 28.9 KB
427f1205   Jean-Michel Glorian   version 4.2 merged
1
2
3
4
FUNCTION dustem_mpfit_run,x,p_min,_EXTRA=extra

;+
; NAME:
36e8b879   Jean-Michel Glorian   update from deborah
5
;    dustem_mpfit_run_vg
427f1205   Jean-Michel Glorian   version 4.2 merged
6
7
8
9
10
; PURPOSE:
;    function used to fit dustem model with dustem_mpfit_sed
; CATEGORY:
;    Dustem
; CALLING SEQUENCE:
36e8b879   Jean-Michel Glorian   update from deborah
11
;    res=dustem_mpfit_run_vg(x,params[,_EXTRA=extra][,cc_sed=cc_sed][,/help])
427f1205   Jean-Michel Glorian   version 4.2 merged
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
; INPUTS:
;    x         = necessary for interface but unused
;    params    = Dustem Model parameters normalized to initial values
; OPTIONAL INPUT PARAMETERS:
;    cc_sed    = SED color correction factors
;    _extra    = extra parameters for the plot routine
; OUTPUTS:
;    None
; OPTIONAL OUTPUT PARAMETERS:
;    None
; ACCEPTED KEY-WORDS:
;    help      = If set, print this help
; COMMON BLOCKS:
;    None
; SIDE EFFECTS:
;    SED and model are plotted if !dustem_show_plot
; RESTRICTIONS:
;    The dustem idl wrapper must be installed
; PROCEDURE:
;    None
; EXAMPLES
36e8b879   Jean-Michel Glorian   update from deborah
33
;    res=dustem_mpfit_run_vg(x,params[,_EXTRA=extra][,cc_sed=cc_sed][,/help])
4750086c   Ilyes Choubani   nouvelle philosph...
34
;  ODIFICATION HISTORY:
427f1205   Jean-Michel Glorian   version 4.2 merged
35
36
37
38
39
40
41
42
;    Adapted to extinction and polarization by Vincent Guillet (IAS)
;    Written by J.-Ph. Bernard
;    see evolution details on the dustem cvs maintained at CESR
;    Contact J.-Ph. Bernard (Jean-Philippe.Bernard@cesr.fr) in case of problems.
;-

;note that x is not used.

427f1205   Jean-Michel Glorian   version 4.2 merged
43
IF keyword_set(help) THEN BEGIN
9be94157   Jean-Philippe Bernard   improved
44
  doc_library,'dustem_mpfit_run'
427f1205   Jean-Michel Glorian   version 4.2 merged
45
46
47
  sed=0.
  goto,the_end
ENDIF
427f1205   Jean-Michel Glorian   version 4.2 merged
48
49

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

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

427f1205   Jean-Michel Glorian   version 4.2 merged
53
;RUN THE MODEL AND GET THE SPECTRUM
759a527d   Ilyes Choubani   general update
54
dustem_activate_plugins,p_dim/(*(*!dustem_fit).param_init_values)
427f1205   Jean-Michel Glorian   version 4.2 merged
55
st = dustem_run(p_dim)
36e8b879   Jean-Michel Glorian   update from deborah
56
st_model=dustem_read_all(!dustem_soft_dir)
427f1205   Jean-Michel Glorian   version 4.2 merged
57

759a527d   Ilyes Choubani   general update
58

4750086c   Ilyes Choubani   nouvelle philosph...
59
60
61
 
;problem that needs fixing - case of the composite ISRF
; st=st if you want to make the isrf thingy work. This is because dustem_activate_plugins requires the st structure but in the case of the composite ISRF the problem is that the wrapper's first read isrf is a MATHIS one so the first abundances the wrapper gets from dustem are the product of dust heated by the MATHIS ISRF
452c334e   Ilyes Choubani   Implementation Of...
62
63


427f1205   Jean-Michel Glorian   version 4.2 merged
64
65
66
chi2_sed      = 0  
chi2_ext      = 0  
chi2_polext   = 0 
36e8b879   Jean-Michel Glorian   update from deborah
67
chi2_polsed   = 0 
452c334e   Ilyes Choubani   Implementation Of...
68
69
70
71
72
chi2_polfrac   = 0
chi2_qsed      =0
chi2_used      =0 
chi2_qext      =0
chi2_uext      =0 
427f1205   Jean-Michel Glorian   version 4.2 merged
73
74
75
rchi2_sed     = 0
rchi2_ext     = 0
rchi2_polext  = 0
36e8b879   Jean-Michel Glorian   update from deborah
76
77
rchi2_polsed  = 0
rchi2_polfrac  = 0
452c334e   Ilyes Choubani   Implementation Of...
78
79
80
81
rchi2_qsed     =0
rchi2_used     =0
rchi2_qext     =0
rchi2_uext     =0
427f1205   Jean-Michel Glorian   version 4.2 merged
82
83
84
wrchi2_sed    = 0
wrchi2_ext    = 0
wrchi2_polext = 0
36e8b879   Jean-Michel Glorian   update from deborah
85
86
wrchi2_polsed = 0
wrchi2_polfrac = 0
452c334e   Ilyes Choubani   Implementation Of...
87
88
89
90
wrchi2_qsed    =0
wrich2_used    =0
wrchi2_qext    =0
wrich2_uext    =0
427f1205   Jean-Michel Glorian   version 4.2 merged
91
92
93
n_sed         = 0
n_ext         = 0
n_polext      = 0
36e8b879   Jean-Michel Glorian   update from deborah
94
95
n_polsed      = 0
n_polfrac      = 0
452c334e   Ilyes Choubani   Implementation Of...
96
97
98
99
n_qsed         =0
n_used         =0
n_qext         =0
n_uext         =0
36e8b879   Jean-Michel Glorian   update from deborah
100
101

win = 0
427f1205   Jean-Michel Glorian   version 4.2 merged
102
103
104

tagnames=tag_names(!dustem_data)
dustem_out = [0]
36e8b879   Jean-Michel Glorian   update from deborah
105
106
107
ind_data_tag=dustem_data_tag(count=count_data_tag)

;if exist('./noplot') then !dustem_show_plot = 0 else if !dustem_show_plot eq 0 then !dustem_show_plot = 1
427f1205   Jean-Michel Glorian   version 4.2 merged
108

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

6730c3f8   Jean-Philippe Bernard   modified to be co...
111
112
;stop

36e8b879   Jean-Michel Glorian   update from deborah
113
114
;print,i_tag,ind_data_tag(i_tag)
;print,tagnames(ind_data_tag(i_tag))
4fd64cbb   Ilyes Choubani   dustem_fit_sed_po...
115
116


36e8b879   Jean-Michel Glorian   update from deborah
117
    CASE tagnames(ind_data_tag(i_tag)) OF
427f1205   Jean-Michel Glorian   version 4.2 merged
118

b5ccb706   Jean-Philippe Bernard   improved to fit p...
119
		'SED'	: BEGIN
452c334e   Ilyes Choubani   Implementation Of...
120
            
452c334e   Ilyes Choubani   Implementation Of...
121
	        dustem_sed = dustem_compute_sed(p_dim,st);dustem_compute_sed(p_dim,st,cont=cont,freefree=freefree,synchrotron=synchrotron)
b5ccb706   Jean-Philippe Bernard   improved to fit p...
122
123
124
125
126
127
128
129
130
131
132
133
134
	        chi2_sed   = total(((*!dustem_data.sed).values-dustem_sed)^2/(*!dustem_data.sed).sigma^2)
			n_sed = n_elements((*!dustem_data.sed).values)
	       	rchi2_sed  = chi2_sed / n_sed
	       	wrchi2_sed  = rchi2_sed * !fit_rchi2_weight.sed
	       	message,"chi2 SED = "+strtrim(chi2_sed,2)+"  rchi2 SED = "+strtrim(rchi2_sed,2),/continue;," weighted rchi2 SED=",wrchi2_sed
			;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
135
		
b5ccb706   Jean-Philippe Bernard   improved to fit p...
136
	        print,"Best SED chi2 obtained if f_HI = ",alpha, " or if we multiply model by ",beta
759a527d   Ilyes Choubani   general update
137
138
	        (*!dustem_fit).chi2=chi2_sed
	        (*!dustem_fit).rchi2=rchi2_sed
36e8b879   Jean-Michel Glorian   update from deborah
139

b5ccb706   Jean-Philippe Bernard   improved to fit p...
140
141
			; Plot if needed
	        IF !dustem_show_plot NE 0 THEN BEGIN
452c334e   Ilyes Choubani   Implementation Of...
142
143
144
	            ;title='DUSTEM WRAP('+strtrim(string(win),2)+')'
	            window,win,title='DUSTEM WRAP (SED)' 
	            
452c334e   Ilyes Choubani   Implementation Of...
145
146
	            ;cgwindow,WMULTI=[2,2,1],/CURRENT
                
452c334e   Ilyes Choubani   Implementation Of...
147
148
149
	            dustem_plot_fit_sed,st,dustem_sed,_extra=extra,res=p_min*(*(*!dustem_fit).param_init_values),chi2=(*!dustem_fit).chi2,rchi2=(*!dustem_fit).rchi2
	           
	            
9be94157   Jean-Philippe Bernard   improved
150
	            ;stop
b5ccb706   Jean-Philippe Bernard   improved to fit p...
151
152
153
		 	ENDIF

	        dustem_out = [ dustem_out, dustem_sed ]
36e8b879   Jean-Michel Glorian   update from deborah
154
155
156

		END

b5ccb706   Jean-Philippe Bernard   improved to fit p...
157
		'EXT' 	: BEGIN 
452c334e   Ilyes Choubani   Implementation Of...
158
159
    		
    		;ADDING PLUGIN TO SPECTRUM----------------
759a527d   Ilyes Choubani   general update
160
    		IF ptr_valid(*!dustem_plugin) THEN BEGIN
452c334e   Ilyes Choubani   Implementation Of...
161
    		for i=0L,n_tags(*!dustem_scope)-1 do begin
4fd64cbb   Ilyes Choubani   dustem_fit_sed_po...
162
            if total(strsplit(*(*!dustem_scope).(i),'+',/extract) eq 'ADD_EXT') then st.ext.ext_tot+=(*(*!dustem_plugin).(i))
452c334e   Ilyes Choubani   Implementation Of...
163
            endfor
759a527d   Ilyes Choubani   general update
164
            endif
452c334e   Ilyes Choubani   Implementation Of...
165
166
        	;------------------------------------------
        	
b5ccb706   Jean-Philippe Bernard   improved to fit p...
167
			dustem_ext = interpol(st.ext.ext_tot,st.ext.wav,(*!dustem_data.ext).wav) 
452c334e   Ilyes Choubani   Implementation Of...
168
            
b5ccb706   Jean-Philippe Bernard   improved to fit p...
169
170
171
172
173
174
	        chi2_ext   = total(((*!dustem_data.ext).values-dustem_ext)^2/(*!dustem_data.ext).sigma^2)
			n_ext = n_elements((*!dustem_data.ext).values)
			rchi2_ext  = chi2_ext / n_ext
	        	wrchi2_ext  = rchi2_ext * !fit_rchi2_weight.ext
	        	print,"chi2 EXT = ",chi2_ext,"  rchi2 EXT = ",rchi2_ext;," weighted rchi2 EXT=",wrchi2_ext
	        	;print,"chi2 per wl: ",((*!dustem_data.ext).values-dustem_ext)^2/(*!dustem_data.ext).sigma^2
36e8b879   Jean-Michel Glorian   update from deborah
175

b5ccb706   Jean-Philippe Bernard   improved to fit p...
176
177
			; Plot if needed
			IF !dustem_show_plot NE 0 THEN BEGIN
36e8b879   Jean-Michel Glorian   update from deborah
178

b5ccb706   Jean-Philippe Bernard   improved to fit p...
179
180
181
182
183
184
				win+=1
	                	;dustem_plot_ext,st,yr=[1.e-6,1.e0],/ysty,xr=[0.5,40],/xsty,win=win
				xr=[0.3,40]
				yr=[1.e-26,1.e-21]
				dustem_plot_ext,st,yr=yr,/ysty,xr=xr,/xsty,ps=filename,win=win,mergePAH=mergePAH,/donotclose,multi=1
				dustem_plot_ext,st,yr=yr,/ysty,xr=xr,/xsty,ps=filename,win=win,mergePAH=mergePAH,multi=2
452c334e   Ilyes Choubani   Implementation Of...
185
        		
b5ccb706   Jean-Philippe Bernard   improved to fit p...
186
187
188
189
190
191
192
193
194
				win+=1
	                	;dustem_plot_ext,st,/UV,yr=[0,2.5],/ysty,xr=[0,10],/xsty,win=win
				xr=[0.2,10]
				yr=[0,3e-21]
				dustem_plot_ext,st,/UV,yr=yr,/ysty,xr=xr,/xsty,ps=filename,win=win,mergePAH=mergePAH,/donotclose,multi=1
				dustem_plot_ext,st,/UV,yr=yr,/ysty,xr=xr,/xsty,ps=filename,win=win,mergePAH=mergePAH,multi=2

			ENDIF
			
759a527d   Ilyes Choubani   general update
195
			;dustem_plot_ext,st,/print_Rv ;commented this recently
b5ccb706   Jean-Philippe Bernard   improved to fit p...
196
197
198
199
200
201
202

		 	dustem_out = [ dustem_out, dustem_ext ]

		 	END

		'POLEXT': BEGIN

759a527d   Ilyes Choubani   general update
203
204
			
			;---------------the procedure within this block has not been updated for a vert long time
b5ccb706   Jean-Philippe Bernard   improved to fit p...
205
206
207
208
209
			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
210
211
212
			;--------------------------------------------------------------------------------------
			
			
b5ccb706   Jean-Philippe Bernard   improved to fit p...
213
214

			IF !run_lin then begin
759a527d   Ilyes Choubani   general update
215
                    ;compute_polext needs updating. It's been months. 
452c334e   Ilyes Choubani   Implementation Of...
216
                    dustem_polext = dustem_compute_polext(p_dim,st,dustem_qext,dustem_uext,Q_ext,U_ext)
b5ccb706   Jean-Philippe Bernard   improved to fit p...
217
	        		chi2_polext   = total(((*!dustem_data.polext).values-dustem_polext)^2/(*!dustem_data.polext).sigma^2)
452c334e   Ilyes Choubani   Implementation Of...
218
    				n_polext = n_elements((*!dustem_data.polext).values)
b5ccb706   Jean-Philippe Bernard   improved to fit p...
219
220
221
222
223
224
225
226
227
	        		rchi2_polext  = chi2_polext / n_polext
	        		wrchi2_polext  = rchi2_polext * !fit_rchi2_weight.polext
	        		print,"chi2 POLEXT = ",chi2_polext,"  rchi2 POLEXT = ",rchi2_polext;," weighted rchi2 POLEXT=",wrchi2_polext
	        		;print,"chi2 per wl: ",((*!dustem_data.polext).values-dustem_polext)^2/(*!dustem_data.polext).sigma^2

	        		; Fit in (pmax,lmax,K)
	        		x=st.polext.wav
	        		logx=alog(x)
	        		logy=alog(st.polext.ext_tot)
4750086c   Ilyes Choubani   nouvelle philosph...
228
	        		
b5ccb706   Jean-Philippe Bernard   improved to fit p...
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
	        		; Whittet+1992 (Table 1)
	        		Uband = 0.36
	        		Bband = 0.43
	        		Vband = 0.55
	        		Rband = 0.63
	        		Iband = 0.78
	        		Jband = 1.21
	        		Hband = 1.64
	        		Kband = 2.04
	        		bands = [Uband,Bband,Vband,Rband,Iband,Jband,Hband]
	        		nbands = n_elements(bands)
	        		jband = indgen(nbands)*0
	        		for k=0,nbands-1 do begin
	                		xx=min(abs(st.polext.wav-bands(k)),j)
	                		jband(k) = j
	        		endfor
	        		;j=where(st.polext.ext_tot gt 0 and x gt 0.2 and x lt 1)
	        		res = poly_fit(logx(jband),logy(jband),2)
			
	        		print
	        		print,"Serkowski Fit with free K over ",string(st.polext(jband).wav,format='(100F5.2)')
	        		print,"-------------------------"
	        		K    = -res(2)
	        		lmax = exp(res(1)/(2*K))
	        		pmax = exp(res(0)+K*alog(lmax)^2)
	        		print,"* lmax = ",lmax," micron"
	        		print,"* K    = ",K
	        		print,"* pmax = ",pmax," for NH=1e21"
	        		print
			
				; Plot if needed
				IF !dustem_show_plot NE 0 THEN BEGIN

452c334e   Ilyes Choubani   Implementation Of...
262
263
264
265
266
;  	       			; Polarization fraction in the UV-IR
;  	                win+=1
; 	                           
;  	       			dustem_plot_polar,st,p_dim,aligned=st_model.align.grains.aligned,yr=[0,0.15],/ysty,xr=[0.1,50],/xsty,win=win;multi=0,ps=filename
; ; 					
b5ccb706   Jean-Philippe Bernard   improved to fit p...
267
	   				; Polarisation curve (Serkowski)
452c334e   Ilyes Choubani   Implementation Of...
268
269
270
271
 					win+=1
 					;xr=[0.3,1e4]
 					xr=[1e-4,3.5]
 					yr=[1e-30,5e-24]
b7f936fc   Ilyes Choubani   to JP: initial fi...
272
273
274
275
276
 					
 					dustem_plot_polext,st,p_dim,dustem_polext,aligned=st_model.align.grains.aligned,win=win
 					
;  					dustem_plot_polar,st,p_dim,aligned=st_model.align.grains.aligned,/UV,yr=yr,/ysty,xr=xr,/xsty,ps=filename,win=win,/donotclose,multi=1
;  					dustem_plot_polar,st,p_dim,aligned=st_model.align.grains.aligned,/UV,yr=yr,/ysty,xr=xr,/xsty,ps=filename,win=win,multi=2
b5ccb706   Jean-Philippe Bernard   improved to fit p...
277
		
452c334e   Ilyes Choubani   Implementation Of...
278
279
280
;  					win+=1
;  					dustem_plot_fpol,win=win,xr=[10,1000],yr=[0,1.1]
;  					
b5ccb706   Jean-Philippe Bernard   improved to fit p...
281
282
	                                
	            ENDIF
36e8b879   Jean-Michel Glorian   update from deborah
283

b5ccb706   Jean-Philippe Bernard   improved to fit p...
284
285
286
				;dustem_plot_polar,st,/print_ratio,cont=cont
	                     
				dustem_out = [ dustem_out, dustem_polext ]
36e8b879   Jean-Michel Glorian   update from deborah
287

b5ccb706   Jean-Philippe Bernard   improved to fit p...
288
			ENDIF
36e8b879   Jean-Michel Glorian   update from deborah
289

b5ccb706   Jean-Philippe Bernard   improved to fit p...
290
			END
427f1205   Jean-Michel Glorian   version 4.2 merged
291

452c334e   Ilyes Choubani   Implementation Of...
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
;USE OF POLFRAC DEPRECATED---------------------------------------

; 		'POLFRAC': BEGIN 
;             
;             ;PLUGIN IS ADDED INSIDE DUSTEM_COMPUTE_POLSED AND DUSTEM_COMPUTE_SED
;             
; 			IF !run_lin then begin
; 
; 	            dustem_polsed = dustem_compute_polsed(p_dim,st);dustem_compute_polsed(p_dim,st.polsed,cont=cont,freefree=freefree,synchrotron=synchrotron)
; 				dustem_sed_tmp = dustem_compute_sed(p_dim,st);dustem_compute_sed(p_dim,st,cont=cont,freefree=freefree,synchrotron=synchrotron)
; 				dustem_sed = dustem_polsed
; 				ind_filt=where((*!dustem_data.polsed).filt_names NE 'SPECTRUM',count_filt)
; 				for i = 0, count_filt-1 do begin
; 					j=where((*!dustem_data.polsed).filt_names(ind_filt(i)) eq (*!dustem_data.sed).filt_names,count)
; 					if count eq 1 then dustem_sed(i) = dustem_sed_tmp(j(0))
; 				endfor 
; 				dustem_polfrac = dustem_polsed/dustem_sed
; 
; 			; We normalize at 353GHz
; 			;x=min(((*!dustem_data.polfrac).wav-850)^2,j353)
; 			;dustem_polfrac = dustem_polfrac/dustem_polfrac(j353)*(*!dustem_data.polfrac).values(j353)
; 
; 	        	chi2_polfrac   = total(((*!dustem_data.polfrac).values-dustem_polfrac)^2/(*!dustem_data.polfrac).sigma^2)
; 				n_polfrac = n_elements((*!dustem_data.polfrac).values)
; 	        	rchi2_polfrac  = chi2_polfrac / n_polfrac
; 	        	wrchi2_polfrac  = rchi2_polfrac * !fit_rchi2_weight.polfrac
; 	        	print,"chi2 POLFRAC = ",chi2_polfrac,"  rchi2 POLFRAC = ",rchi2_polfrac;," weighted rchi2 POLFRAC=",wrchi2_polfrac
; 	        	;print,"chi2 per wl: ",((*!dustem_data.polfrac).values-dustem_polfrac)^2/(*!dustem_data.polfrac).sigma^2
; 
; 	            dustem_out = [ dustem_out, dustem_polfrac ]
; 
; 			; Plot if needed
; 				IF !dustem_show_plot NE 0 THEN BEGIN
; 
; 				; Also show prediction for polarized P/I
; 					win+=1
; 	       			dustem_plot_polar,st,aligned=st_model.align.grains.aligned,/Pfrac,win=win,yr=[0,0.25],/ysty,xr=[10,5d3];dustem_plot_polar,st,aligned=st_model.align.grains.aligned,/Pfrac,win=win,cont=cont,freefree=freefree,synchrotron=synchrotron,yr=[0,0.25],/ysty,xr=[10,5d3]
; 
; 		 		ENDIF
; 
; 		 	ENDIF
; 
; 			END


;-----------------------------------------------------------------
b5ccb706   Jean-Philippe Bernard   improved to fit p...
338
		'POLSED': BEGIN
452c334e   Ilyes Choubani   Implementation Of...
339
340
341
		    
		    ;PLUGIN IS ADDED INSIDE DUSTEM_COMPUTE_SED
			
b5ccb706   Jean-Philippe Bernard   improved to fit p...
342
343
			IF !run_lin THEN BEGIN

4fd64cbb   Ilyes Choubani   dustem_fit_sed_po...
344
	        	dustem_polsed = dustem_compute_polsed(p_dim,st);dustem_compute_polsed(p_dim,st,cont=cont,freefree=freefree,synchrotron=synchrotron)
452c334e   Ilyes Choubani   Implementation Of...
345
	        
b5ccb706   Jean-Philippe Bernard   improved to fit p...
346
347
348
349
350
351
352
	        	chi2_polsed   = total(((*!dustem_data.polsed).values-dustem_polsed)^2/(*!dustem_data.polsed).sigma^2)
				n_polsed = n_elements((*!dustem_data.polsed).values)
	        	rchi2_polsed  = chi2_polsed / n_polsed
	        	wrchi2_polsed  = rchi2_polsed * !fit_rchi2_weight.polsed
	        	message,"chi2 POLSED = "+strtrim(chi2_polsed,2)+"  rchi2 POLSED = "+strtrim(rchi2_polsed,2),/continue;," weighted rchi2 POLSED=",wrchi2_polsed
	        	print,"chi2 per wl: ",((*!dustem_data.polsed).values-dustem_polsed)^2/(*!dustem_data.polsed).sigma^2
	        	;print,"SNR per wl: ",(*!dustem_data.polsed).values/(*!dustem_data.polsed).sigma
452c334e   Ilyes Choubani   Implementation Of...
353
354
355
356
              
	            dustem_out = [ dustem_out, dustem_polsed ] ; classic command
		        ;-----------------------------------------------
	         
b5ccb706   Jean-Philippe Bernard   improved to fit p...
357
358
				; Plot if needed
				IF !dustem_show_plot NE 0 THEN BEGIN
b5ccb706   Jean-Philippe Bernard   improved to fit p...
359
					win+=1
452c334e   Ilyes Choubani   Implementation Of...
360
361
; 					
					xr=[10,8e3]
b5ccb706   Jean-Philippe Bernard   improved to fit p...
362
					yr=[1e-34,1e-28]
452c334e   Ilyes Choubani   Implementation Of...
363
					
759a527d   Ilyes Choubani   general update
364
365
366
					dustem_plot_polsed,st,p_dim,dustem_polsed,aligned=st_model.align.grains.aligned,win=win
; 					dustem_plot_polar,st,p_dim,aligned=st_model.align.grains.aligned,/SED,ps=filename,win=win,cont=cont,freefree=freefree,synchrotron=synchrotron,xr=xr,yr=yr,/xsty,/donotclose,multi=1
; 					dustem_plot_polar,st,p_dim,aligned=st_model.align.grains.aligned,/SED,ps=filename,win=win,cont=cont,freefree=freefree,synchrotron=synchrotron,xr=xr,yr=yr,/xsty,multi=2
452c334e   Ilyes Choubani   Implementation Of...
367
368
369

                    					
					
b5ccb706   Jean-Philippe Bernard   improved to fit p...
370
				ENDIF
b5ccb706   Jean-Philippe Bernard   improved to fit p...
371
			ENDIF
427f1205   Jean-Michel Glorian   version 4.2 merged
372
		END
452c334e   Ilyes Choubani   Implementation Of...
373
374
375
        
        'QSED': BEGIN
		    
607060e5   Ilyes Choubani   test version
376
		        toto = dustem_compute_stokes(p_dim,st,dustem_qsed,dustem_used,Q_spec,U_spec) 
452c334e   Ilyes Choubani   Implementation Of...
377
378
379
380
381
382
383
384
385
386
387
			    chi2_qsed   =total(((*!dustem_data.qsed).values-dustem_qsed)^2/(*!dustem_data.qsed).sigma^2) 
			    n_qsed = n_elements((*!dustem_data.qsed).values)
			    rchi2_qsed  = chi2_qsed / n_qsed
	        	wrchi2_qsed  = rchi2_qsed * !fit_rchi2_weight.qsed 
	            dustem_out = [ dustem_out, dustem_qsed ] ; classic command
	            message,"chi2 QSED = "+strtrim(chi2_qsed,2)+"  rchi2 QSED = "+strtrim(rchi2_qsed,2),/continue;," weighted rchi2 POLSED=",wrchi2_polsed
	        	print,"chi2 per wl: ",((*!dustem_data.qsed).values-dustem_qsed)^2/(*!dustem_data.qsed).sigma^2
		        
		      
		        
		        if !dustem_show_plot NE 0 then begin
d012e324   Ilyes Choubani   FIX to the use of...
388
389
		        fact = 1.00E-20;(3.e10/(1.e-4*(st.polsed.wav)))/1.d40
    	        fact1 = 1.00E-20;(3.e10/(1.e-4*((*!dustem_data.qsed).wav)))/1.d40 ; for data plotting   
452c334e   Ilyes Choubani   Implementation Of...
390
	            
d012e324   Ilyes Choubani   FIX to the use of...
391
392
393
; 	            fact=1/(4.*!pi)/(3.e8/1.e-6/(st.polsed).wav)*1.00e17
;                 fact1=1.E-20
; 	            
452c334e   Ilyes Choubani   Implementation Of...
394
395
396
		        win+=1
                window  ,win ,title='DUSTEM WRAP (QSED)'
		        titstq='Stokes Q Polarization Intensity'
d012e324   Ilyes Choubani   FIX to the use of...
397
                ytitstq=textoidl('Q_\nu (W/m^2/HZ/sr) for N_H=10^{20} H/cm^2')
452c334e   Ilyes Choubani   Implementation Of...
398
                xtit=textoidl('\lambda (\mum)')
607060e5   Ilyes Choubani   test version
399
                xr=[10,2e4]
4fd64cbb   Ilyes Choubani   dustem_fit_sed_po...
400
                normpol=Q_spec*fact
452c334e   Ilyes Choubani   Implementation Of...
401
402
                normpol1=dustem_qsed*fact1
                ;normpol1=(*!dustem_data.qsed).values*fact1
d012e324   Ilyes Choubani   FIX to the use of...
403
404
                yr=[1e-28,1e-17]
                ;yr=[1e-34,1e-23]
4fd64cbb   Ilyes Choubani   dustem_fit_sed_po...
405
                indq=where(Q_spec lt 0, countq)
759a527d   Ilyes Choubani   general update
406
407
                ;indq=where((*!dustem_data.qsed).values lt 0, countq)
                
452c334e   Ilyes Choubani   Implementation Of...
408
409
410
                ylog=1
                if countq ne 0 then begin
                ylog=0
759a527d   Ilyes Choubani   general update
411
                yr=[-2e-31,5e-31] 
452c334e   Ilyes Choubani   Implementation Of...
412
413
414
                endif
                
                
d012e324   Ilyes Choubani   FIX to the use of...
415
                
4fd64cbb   Ilyes Choubani   dustem_fit_sed_po...
416
417
                cgplot,st.polsed.wav,Q_spec*fact,xtit='',ytit=ytitstq,tit=titstq,ylog=ylog,/xlog,xr=xr,yr=yr,/ys,/xs,position=[0.17,0.35,0.93,0.93],xtickformat='(A1)'
                cgoplot,st.polsed.wav,Q_spec*fact,color='black'
452c334e   Ilyes Choubani   Implementation Of...
418
419
420
421
                cgoplot,(*!dustem_data.qsed).wav,(*!dustem_data.qsed).values*fact1,color='Tomato',psym=16,symsize=1,thick=2
                err_bar,(*!dustem_data.qsed).wav,(*!dustem_data.qsed).values*fact1,yrms=3.*((*!dustem_data.qsed).sigma)/2.*fact1,color=cgColor('Tomato')
                cgoplot,(*!dustem_data.qsed).wav,dustem_qsed*fact1,psym=6,color='red',symsize=2
                
4fd64cbb   Ilyes Choubani   dustem_fit_sed_po...
422
423
                cgplot,st.polsed.wav,Q_spec*fact/normpol,xtit=xtit,ytit='Normalized',tit='',/xlog,xr=xr,/ys,/xs,yr=[0,2],ylog=0,position=[0.17,0.14,0.93,0.35],/noerase,yticks=2,ymino=5,xticklen=0.1    
                cgoplot,st.polsed.wav,Q_spec*fact/normpol,color='black'
452c334e   Ilyes Choubani   Implementation Of...
424
425
426
427
428
429
430
431
432
                cgoplot,(*!dustem_data.qsed).wav,(*!dustem_data.qsed).values*fact1/normpol1,color='Tomato',psym=16,symsize=1,thick=2
                err_bar,(*!dustem_data.qsed).wav,(*!dustem_data.qsed).values*fact1/normpol1,yrms=3.*((*!dustem_data.qsed).sigma)/2.*fact1/normpol1,color=cgColor('Tomato')
                        
                endif
		        
 		END

        'USED': BEGIN
		    
607060e5   Ilyes Choubani   test version
433
                toto = dustem_compute_stokes(p_dim,st,dustem_qsed,dustem_used,Q_spec,U_spec)
452c334e   Ilyes Choubani   Implementation Of...
434
435
436
437
438
439
440
441
442
443
			    chi2_used   =total(((*!dustem_data.used).values-dustem_used)^2/(*!dustem_data.used).sigma^2) 
 			    n_used = n_elements((*!dustem_data.used).values)
 			    rchi2_used  = chi2_used / n_used
 	        	wrchi2_used  = rchi2_used * !fit_rchi2_weight.used
 	            dustem_out = [ dustem_out, dustem_used] ; classic command
 	            message,"chi2 USED = "+strtrim(chi2_used,2)+"  rchi2 USED = "+strtrim(rchi2_used,2),/continue;," weighted rchi2 POLSED=",wrchi2_polsed
 	        	print,"chi2 per wl: ",((*!dustem_data.used).values-dustem_used)^2/(*!dustem_data.used).sigma^2
 
	            
	            if !dustem_show_plot ne 0 then begin
d012e324   Ilyes Choubani   FIX to the use of...
444
445
	            fact = 1.00E-20;(3.e10/(1.e-4*(st.polsed.wav)))/1.d40
    	        fact1 =1.00E-20;(3.e10/(1.e-4*((*!dustem_data.used).wav)))/1.d40 ; for data plotting   
452c334e   Ilyes Choubani   Implementation Of...
446
447
448
449
	     
	            win+=1
                window  ,win ,title='DUSTEM WRAP (USED)'                 
                titstu='Stokes U Polarization Intensity'
d012e324   Ilyes Choubani   FIX to the use of...
450
                ytitstu=textoidl('U_\nu (W/m^2/Hz/sr) for N_H=10^{20} H/cm^2')
452c334e   Ilyes Choubani   Implementation Of...
451
                xtit=textoidl('\lambda (\mum)')
ba82c5cf   Ilyes Choubani   improved with mod...
452
    	        xr=[10,6e4]
4fd64cbb   Ilyes Choubani   dustem_fit_sed_po...
453
    	        normpol2=U_spec*fact
452c334e   Ilyes Choubani   Implementation Of...
454
455
                normpol3=dustem_used*fact1
                ;yr=[min(normpol2)-min(normpol2)/10,max(normpol2)+max(normpol2)*10]	
759a527d   Ilyes Choubani   general update
456
457
458
                 
                 ;conditions here are on computed spectra and they should be on the data itself.
                 ;This even impairs the visualization of the data. Correct this!!!!!
452c334e   Ilyes Choubani   Implementation Of...
459
                  
d012e324   Ilyes Choubani   FIX to the use of...
460
                 yr=[1e-28,1e-17];[1e-34,1e-23]
4fd64cbb   Ilyes Choubani   dustem_fit_sed_po...
461
                 indu=where(U_spec lt 0, countu)
759a527d   Ilyes Choubani   general update
462
                 ;indu=where((*!dustem_data.used).values lt 0, countu)
452c334e   Ilyes Choubani   Implementation Of...
463
464
465
                 ylog=1
                 if countu ne 0 then begin
                 ylog=0
759a527d   Ilyes Choubani   general update
466
                 yr=[-2e-31,5e-31]  
452c334e   Ilyes Choubani   Implementation Of...
467
468
                 endif 
                  
4fd64cbb   Ilyes Choubani   dustem_fit_sed_po...
469
470
                cgplot,st.polsed.wav,U_spec*fact,xtit='',ytit=ytitstu,tit=titstu,ylog=ylog,/xlog,xr=xr,yr=yr,/ys,/xs,position=[0.17,0.35,0.93,0.93],xtickformat='(A1)',color='black'
                cgoplot,st.polsed.wav,U_spec*fact,color='black'
452c334e   Ilyes Choubani   Implementation Of...
471
472
473
474
475
476
                cgoplot,(*!dustem_data.used).wav,(*!dustem_data.used).values*fact1,color='Teal',psym=16,symsize=1,thick=2
                err_bar,(*!dustem_data.used).wav,(*!dustem_data.used).values*fact1,yrms=3.*((*!dustem_data.used).sigma)/2.*fact1,color=cgColor('Teal')
                cgoplot,(*!dustem_data.used).wav,dustem_used*fact1,psym=6,color='red',symsize=2
                
                
                
4fd64cbb   Ilyes Choubani   dustem_fit_sed_po...
477
478
                cgplot,st.polsed.wav,U_spec*fact/normpol2,xtit=xtit,ytit='Normalized',tit='',/xlog,xr=xr,/ys,/xs,yr=[0,2],ylog=0,position=[0.17,0.14,0.93,0.35],/noerase,yticks=2,ymino=5,xticklen=0.1    
                cgoplot,st.polsed.wav,U_spec*fact/normpol2,color='black'
452c334e   Ilyes Choubani   Implementation Of...
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
                cgoplot,(*!dustem_data.used).wav,(*!dustem_data.used).values*fact1/normpol3,color='Teal',psym=16,symsize=1,thick=2
                err_bar,(*!dustem_data.used).wav,(*!dustem_data.used).values*fact1/normpol3,yrms=3.*((*!dustem_data.used).sigma)/2.*fact1/normpol3,color=cgColor('Teal')
                                
                ;cgwindow, 'plot', 
                	 
	            endif
 		END
 		
 		
 		
 		
  		'QEXT': BEGIN ; Q - EXTINCTION CROSS SECTION
                   
 			    chi2_qext   =total(((*!dustem_data.qext).values-dustem_qext)^2/(*!dustem_data.qext).sigma^2) 
   			    n_qext = n_elements((*!dustem_data.qext).values)
   			    rchi2_qext  = chi2_qext / n_qext
  	        	wrchi2_qext  = rchi2_qext * !fit_rchi2_weight.qext
   	            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
  	        	print,"chi2 per wl: ",((*!dustem_data.qext).values-dustem_qext)^2/(*!dustem_data.qext).sigma^2
                if !dustem_show_plot ne 0 then begin
	            
	            win+=1
                window  ,win ,title='DUSTEM WRAP (QEXT)'	                 
                titqex='Q-polarized Cross-section'
                ytitqex=textoidl('\sigma_{Q}/N_H (cm^2/H)')
                xtit=textoidl('\lambda^{-1} (\mum^{-1})')
    	        ;xr=[0.3,1e4]
    	        xr=[1e-4,3.5]
    	          
                yr=[1e-34,5e-24]
                indqex=where(Q_ext lt 0, countqex)
                ylog=1
                if countqex ne 0 then begin
                ylog=0
                yr=[-9e-27,5e-27]  
                endif 
                normpol = st.polext.ext_tot * 0. + 1.d21  
                cgplot,1/st.polext.wav,Q_ext,xtit='',ytit=ytitqex,tit=titqex,ylog=ylog,/xlog,xr=xr,yr=yr,/ys,/xs,position=[0.17,0.35,0.93,0.93],xtickformat='(A1)',color='black'
                cgoplot,1/st.polext.wav,Q_ext/normpol,color='black'
                cgoplot,1/(*!dustem_data.qext).wav,(*!dustem_data.qext).values/normpol,color='Indian Red',psym=16,symsize=1,thick=2
                err_bar,1/(*!dustem_data.qext).wav,(*!dustem_data.qext).values/normpol,yrms=3.*((*!dustem_data.qext).sigma)/2/normpol,color=cgColor('Indian Red')
                cgoplot,1/(*!dustem_data.qext).wav,dustem_qext/normpol,psym=6,color='red',symsize=2
                
                
                
                cgplot,1/st.polsed.wav,Q_ext/Q_ext,xtit=xtit,ytit='Normalized',tit='',/xlog,xr=xr,/ys,/xs,yr=[0,2],ylog=0,position=[0.17,0.14,0.93,0.35],/noerase,yticks=2,ymino=5,xticklen=0.1    
                cgoplot,1/st.polsed.wav,Q_ext/Q_ext,color='black'
                cgoplot,1/(*!dustem_data.qext).wav,(*!dustem_data.qext).values/dustem_qext,color='Indian Red',psym=16,symsize=1,thick=2
                err_bar,1/(*!dustem_data.qext).wav,(*!dustem_data.qext).values/dustem_qext,yrms=3.*((*!dustem_data.qext).sigma)/2.,color=cgColor('Indian Red')
                    
            	            endif                  
     	           
  		END
;  		
  		'UEXT': BEGIN
  
 			    
 			    chi2_uext   =total(((*!dustem_data.uext).values-dustem_uext)^2/(*!dustem_data.uext).sigma^2) 
   			    n_uext = n_elements((*!dustem_data.uext).values)
   			    rchi2_uext  = chi2_uext / n_uext
  	        	wrchi2_uext  = rchi2_uext * !fit_rchi2_weight.uext
   	            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
  	        	print,"chi2 per wl: ",((*!dustem_data.uext).values-dustem_uext)^2/(*!dustem_data.uext).sigma^2
                if !dustem_show_plot ne 0 then begin
	            
	            win+=1
                window  ,win ,title='DUSTEM WRAP (UEXT)'	                 
                tituex='U-polarized Cross-section'
                ytituex=textoidl('\sigma_{U}/N_H (cm^2/H)')
                xtit=textoidl('\lambda^{-1} (\mum^{-1})')
    	        ;xr=[0.3,1e4]
    	        xr=[1e-4,3.5]
    	          
                yr=[1e-34,5e-24]
                induex=where(U_ext lt 0, countuex)
                ylog=1
                if countuex ne 0 then begin
                ylog=0
                yr=[-9e-27,5e-27]  
                endif 
                normpol = st.polext.ext_tot * 0. + 1.d21  
                cgplot,1/st.polext.wav,U_ext,xtit='',ytit=ytituex,tit=tituex,ylog=ylog,/xlog,xr=xr,yr=yr,/ys,/xs,position=[0.17,0.35,0.93,0.93],xtickformat='(A1)',color='black'
                cgoplot,1/st.polext.wav,U_ext/normpol,color='black'
                cgoplot,1/(*!dustem_data.uext).wav,(*!dustem_data.uext).values/normpol,color='Steel Blue',psym=16,symsize=1,thick=2
                err_bar,1/(*!dustem_data.uext).wav,(*!dustem_data.uext).values/normpol,yrms=3.*((*!dustem_data.uext).sigma)/2/normpol,color=cgColor('Steel Blue')
                cgoplot,1/(*!dustem_data.uext).wav,dustem_uext/normpol,psym=6,color='red',symsize=2
                
                
                
                cgplot,1/st.polsed.wav,U_ext/U_ext,xtit=xtit,ytit='Normalized',tit='',/xlog,xr=xr,/ys,/xs,yr=[0,2],ylog=0,position=[0.17,0.14,0.93,0.35],/noerase,yticks=2,ymino=5,xticklen=0.1    
                cgoplot,1/st.polsed.wav,U_ext/U_ext,color='black'
                cgoplot,1/(*!dustem_data.uext).wav,(*!dustem_data.uext).values/dustem_uext,color='Steel Blue',psym=16,symsize=1,thick=2
                err_bar,1/(*!dustem_data.uext).wav,(*!dustem_data.uext).values/dustem_uext,yrms=3.*((*!dustem_data.uext).sigma)/2/normpol,color=cgColor('Steel Blue')
                    
            	            endif                  
     	        
                
                  
 	           
  		END
 		
 		
	    ELSE : stop
427f1205   Jean-Michel Glorian   version 4.2 merged
584
585
586

    ENDCASE

b5ccb706   Jean-Philippe Bernard   improved to fit p...
587
ENDFOR
427f1205   Jean-Michel Glorian   version 4.2 merged
588

36e8b879   Jean-Michel Glorian   update from deborah
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
DOF = n_elements(p_min)
total_chi2 = 0
total_rchi2 = 0
total_n = 0
if !fit_rchi2_weight.ext ne 0 then begin
	total_n += n_ext
	total_chi2 +=  chi2_ext
	total_rchi2 += rchi2_ext
endif
if !fit_rchi2_weight.sed ne 0 then begin
	total_n += n_sed
	total_chi2 += chi2_sed
	total_rchi2 += rchi2_sed
endif
if !run_pol then begin
if !fit_rchi2_weight.polext ne 0 then begin
	total_n += n_polext
	total_chi2 += chi2_polext
	total_rchi2 += rchi2_polext
endif
if !fit_rchi2_weight.polsed ne 0 then begin
452c334e   Ilyes Choubani   Implementation Of...
610
611
612
	total_n += n_polsed 
	total_chi2 += chi2_polsed 
	total_rchi2 += rchi2_polsed 
36e8b879   Jean-Michel Glorian   update from deborah
613
614
615
616
617
618
endif
if !fit_rchi2_weight.polfrac ne 0 then begin
	total_n += n_polfrac
	total_chi2 += chi2_polfrac
	total_rchi2 += rchi2_polfrac
endif
452c334e   Ilyes Choubani   Implementation Of...
619
620
621
622
623
624
625
626
627
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
36e8b879   Jean-Michel Glorian   update from deborah
628
endif
452c334e   Ilyes Choubani   Implementation Of...
629
630
631
632
633
634
635
636
637
638
639
640
641
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

endif

36e8b879   Jean-Michel Glorian   update from deborah
642
643
644
645
646
647

total_wchi2 = chi2_ext * !fit_rchi2_weight.ext $
            + chi2_sed * !fit_rchi2_weight.sed 
if !run_pol then begin
	total_wchi2 	+= chi2_polext  * !fit_rchi2_weight.polext $
             		 + chi2_polsed  * !fit_rchi2_weight.polsed $
452c334e   Ilyes Choubani   Implementation Of...
648
649
650
651
652
             		 + 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
653
654
endif

9ccf7615   Jean-Philippe Bernard   modified to fit u...
655
message,'======================================================================',/continue
36e8b879   Jean-Michel Glorian   update from deborah
656
657
658
659
660
661
662
663
664
print,"Total weighted chi2 driving mpfitfun = ", total_wchi2
print,"Nb of data points = ", total_n
print,"DOF = ", DOF
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
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
print,"Mean weighted reduced chi2 = ", total_wchi2/(total_n-DOF)
print,"Unweighted reduced Total chi2 = ", total_chi2/(total_n-DOF)
452c334e   Ilyes Choubani   Implementation Of...
665
print,"Unweighted mean reduced chi2 = ", total_rchi2/4*total_n/(total_n-DOF);total_rchi2/4*total_n/(total_n-DOF) original
9ccf7615   Jean-Philippe Bernard   modified to fit u...
666
message,'======================================================================',/continue
36e8b879   Jean-Michel Glorian   update from deborah
667
668
669
670
671
672
673
674
675

;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...
676
;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
677
678
679
680

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

9ccf7615   Jean-Philippe Bernard   modified to fit u...
681
682
683
684
685
686
687
;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
688
689
690
691
692
693

the_end:

RETURN,dustem_out

END