OV-GSO-DC / dustem-wrapper_idl

Commit 45da29c6ff2c1aade6ae88296092c9a212cb6e41

Authored by Ilyes Choubani
1 parent c8399ce3
Exists in master

updating mpfit_run

Showing 1 changed file with 62 additions and 49 deletions   Show diff stats
Inline Side-by-side
src/idl/dustem_mpfit_run.pro
  Diff comments   View file @45da29c
... ... @@ -51,10 +51,19 @@ p_dim=p_min*(*(*!dustem_fit).param_init_values)
51 51 IF !run_ionfrac gt 0. THEN toto = dustem_create_ionfrac()
52 52  
53 53 ;RUN THE MODEL AND GET THE SPECTRUM
54   -dustem_activate_plugins,p_dim/(*(*!dustem_fit).param_init_values)
55   -st = dustem_run(p_dim)
  54 +
  55 +stop
  56 +dustem_activate_plugins,p_dim/(*(*!dustem_fit).param_init_values),st ;st is an output...
  57 +
  58 +stop
  59 +
  60 +st1 = dustem_run(p_dim)
  61 +
  62 +stop
  63 +
56 64 st_model=dustem_read_all(!dustem_soft_dir)
57 65  
  66 +
58 67 chi2_sed = 0
59 68 chi2_ext = 0
60 69 chi2_qsed =0
... ... @@ -89,11 +98,13 @@ tagnames=tag_names(!dustem_data)
89 98 dustem_out = [0]
90 99 ind_data_tag=dustem_data_tag(count=count_data_tag)
91 100  
  101 +;stop
92 102 FOR i_tag=0,count_data_tag-1 DO BEGIN
93 103  
94 104 CASE tagnames(ind_data_tag(i_tag)) OF
95 105  
96 106 'SED' : BEGIN
  107 +
97 108 dustem_sed = dustem_compute_sed(p_dim,st)
98 109 chi2_sed = total(((*!dustem_data.sed).values-dustem_sed)^2/(*!dustem_data.sed).sigma^2)
99 110 n_sed = n_elements((*!dustem_data.sed).values)
... ... @@ -252,12 +263,56 @@ FOR i_tag=0,count_data_tag-1 DO BEGIN
252 263  
253 264 END
254 265  
255   -;!!!!!!!!!!!!!YOU NEED TO RUN 'COMPUTE_SED' AGAIN because you have no idea if the sed block will be run.
  266 + 'QSED': BEGIN
  267 +
  268 + toto = dustem_compute_stokes(p_dim,st,dustem_qsed,dustem_used,Q_spec,U_spec)
  269 +
  270 + ;taking care of the color correction for spass ; SHOULD NOT LEAVE THIS ON
  271 + ;dustem_qsed(-1)=dustem_qsed(-2)
  272 +; testspass1 = ((*!dustem_data.qsed).filt_names)(-1) eq 'SPASS1'
  273 +; if testspass1 then dustem_qsed(-1)=dustem_qsed(-2)
  274 +;
  275 + chi2_qsed =total(((*!dustem_data.qsed).values-dustem_qsed)^2/(*!dustem_data.qsed).sigma^2)
  276 + n_qsed = n_elements((*!dustem_data.qsed).values)
  277 + rchi2_qsed = chi2_qsed / n_qsed
  278 + wrchi2_qsed = rchi2_qsed * !fit_rchi2_weight.qsed
  279 + if isa(!dustem_data.qsed) then dustem_out = [ dustem_out, dustem_qsed ] ; classic command
  280 + message,"chi2 QSED = "+strtrim(chi2_qsed,2)+" rchi2 QSED = "+strtrim(rchi2_qsed,2),/continue;," weighted rchi2 POLSED=",wrchi2_polsed
  281 + print,"chi2 per wl: ",((*!dustem_data.qsed).values-dustem_qsed)^2/(*!dustem_data.qsed).sigma^2
  282 +
  283 + END
256 284  
257   - 'POLFRAC': BEGIN
  285 + 'USED': BEGIN
  286 + ;THIS IMPLIES THAT QSED HAS BEEN RAN BEFORE
  287 + ;toto = dustem_compute_stokes(p_dim,st,dustem_qsed,dustem_used,Q_spec,U_spec)
  288 + chi2_used =total(((*!dustem_data.used).values-dustem_used)^2/(*!dustem_data.used).sigma^2)
  289 + n_used = n_elements((*!dustem_data.used).values)
  290 + rchi2_used = chi2_used / n_used
  291 + wrchi2_used = rchi2_used * !fit_rchi2_weight.used
  292 + if isa(!dustem_data.used) then dustem_out = [ dustem_out, dustem_used] ; classic command
  293 + message,"chi2 USED = "+strtrim(chi2_used,2)+" rchi2 USED = "+strtrim(rchi2_used,2),/continue;," weighted rchi2 POLSED=",wrchi2_polsed
  294 + print,"chi2 per wl: ",((*!dustem_data.used).values-dustem_used)^2/(*!dustem_data.used).sigma^2
  295 +
  296 + END
  297 +
  298 + 'POLSED': BEGIN
  299 +
  300 + ;PLUGIN IS ADDED INSIDE DUSTEM_COMPUTE_SED
  301 +
  302 + IF !run_pol and !run_lin THEN BEGIN
  303 +
  304 + dustem_polsed = dustem_compute_polsed(p_dim,st);dustem_compute_polsed(p_dim,st,cont=cont,freefree=freefree,synchrotron=synchrotron)
  305 +
  306 + ENDIF
  307 + END
  308 +
  309 +
  310 + 'POLFRAC': BEGIN
  311 +
258 312 if !run_lin then begin
  313 +
259 314 if ~isa(dustem_sed) then dustem_sed = dustem_compute_sed(p_dim,st)
260   - dustem_polsed = dustem_compute_polsed(p_dim,st) ; this quantity hasn't been computed yet that is why this line is present here.
  315 + if ~isa(dustem_polsed) then dustem_polsed = dustem_compute_polsed(p_dim,st) ; this quantity hasn't been computed yet that is why this line is present here.
261 316 ind_filt=where((*!dustem_data.polfrac).filt_names NE 'SPECTRUM',count_filt) ;locating the filters of polfrac=smallp in SED xcat
262 317 ind_fspec=where((*!dustem_data.polfrac).filt_names EQ 'SPECTRUM',count_fspec) ;locating the spectrum points in polfrac data. Later testing will have to be on the wavelengths
263 318 ;
... ... @@ -309,53 +364,10 @@ FOR i_tag=0,count_data_tag-1 DO BEGIN
309 364 dustem_polfrac = dustem_polsed_x/dustem_sed_x
310 365  
311 366 ENDIF
  367 +
312 368 ;stop
313 369 END
314 370  
315   -
316   -;-----------------------------------------------------------------
317   - 'POLSED': BEGIN
318   -
319   - ;PLUGIN IS ADDED INSIDE DUSTEM_COMPUTE_SED
320   -
321   - IF !run_lin THEN BEGIN
322   -
323   - dustem_polsed = dustem_compute_polsed(p_dim,st);dustem_compute_polsed(p_dim,st,cont=cont,freefree=freefree,synchrotron=synchrotron)
324   -
325   - ENDIF
326   - END
327   -
328   - 'QSED': BEGIN
329   -
330   - toto = dustem_compute_stokes(p_dim,st,dustem_qsed,dustem_used,Q_spec,U_spec)
331   -
332   - ;taking care of the color correction for spass ; SHOULD NOT LEAVE THIS ON
333   - ;dustem_qsed(-1)=dustem_qsed(-2)
334   -; testspass1 = ((*!dustem_data.qsed).filt_names)(-1) eq 'SPASS1'
335   -; if testspass1 then dustem_qsed(-1)=dustem_qsed(-2)
336   -;
337   - chi2_qsed =total(((*!dustem_data.qsed).values-dustem_qsed)^2/(*!dustem_data.qsed).sigma^2)
338   - n_qsed = n_elements((*!dustem_data.qsed).values)
339   - rchi2_qsed = chi2_qsed / n_qsed
340   - wrchi2_qsed = rchi2_qsed * !fit_rchi2_weight.qsed
341   - if isa(!dustem_data.qsed) then dustem_out = [ dustem_out, dustem_qsed ] ; classic command
342   - message,"chi2 QSED = "+strtrim(chi2_qsed,2)+" rchi2 QSED = "+strtrim(rchi2_qsed,2),/continue;," weighted rchi2 POLSED=",wrchi2_polsed
343   - print,"chi2 per wl: ",((*!dustem_data.qsed).values-dustem_qsed)^2/(*!dustem_data.qsed).sigma^2
344   -
345   - END
346   -
347   - 'USED': BEGIN
348   - ;THIS IMPLIES THAT QSED HAS BEEN RAN BEFORE
349   - chi2_used =total(((*!dustem_data.used).values-dustem_used)^2/(*!dustem_data.used).sigma^2)
350   - n_used = n_elements((*!dustem_data.used).values)
351   - rchi2_used = chi2_used / n_used
352   - wrchi2_used = rchi2_used * !fit_rchi2_weight.used
353   - if isa(!dustem_data.used) then dustem_out = [ dustem_out, dustem_used] ; classic command
354   - message,"chi2 USED = "+strtrim(chi2_used,2)+" rchi2 USED = "+strtrim(rchi2_used,2),/continue;," weighted rchi2 POLSED=",wrchi2_polsed
355   - print,"chi2 per wl: ",((*!dustem_data.used).values-dustem_used)^2/(*!dustem_data.used).sigma^2
356   - dustemwrap_plot,p_dim,st,dustem_sed,dustem_qsed,dustem_used,dustem_polsed,dustem_polfrac,_extra=_extra
357   -
358   - END
359 371  
360 372 'QEXT': BEGIN ; Q - EXTINCTION CROSS SECTION
361 373  
... ... @@ -453,6 +465,7 @@ FOR i_tag=0,count_data_tag-1 DO BEGIN
453 465 ENDCASE
454 466  
455 467 ENDFOR
  468 +dustemwrap_plot,p_dim,st,dustem_sed,dustem_qsed,dustem_used,dustem_polsed,dustem_polfrac,_extra=_extra
456 469  
457 470 DOF = n_elements(p_min)
458 471 total_chi2 = 0
... ...