diff --git a/src/idl/dustem_fit.pro b/src/idl/dustem_fit.pro index 1c6cc7a..a43b4a5 100755 --- a/src/idl/dustem_fit.pro +++ b/src/idl/dustem_fit.pro @@ -71,17 +71,19 @@ defsysv, '!dustem_data', { $ ;Data to fit ; Define global !fit_rchi2_weight with the same structure and field order like !dustem_data tagnames = tag_names(!dustem_data) -instr="defsysv, '!fit_rchi2_weight', {" +;instr="defsysv, '!fit_rchi2_weight', {" +instr="fit_rchi2_weight_string={" for i = 0, n_elements(tagnames)-1 do instr+=tagnames(i)+': rchi2_weight.'+tagnames(i)+', ' -instr+='empty: 0. } +instr+='empty: 0. }' toto=execute(instr) +defsysv, '!fit_rchi2_weight',ptr_new(fit_rchi2_weight_string) ;=== Set weight for reduced chi2 ponderation between emission, extinction, polarized extinction and emission === ; Calculate total according to fileds in !dustem_data structure total = 0 -for i = 0, n_tags(!dustem_data)-1 do total = total + !fit_rchi2_weight.(i) +for i = 0, n_tags(!dustem_data)-1 do total = total + (*!fit_rchi2_weight).(i) ; values are relative to total -for i = 0, n_tags(!dustem_data)-1 do !fit_rchi2_weight.(i) /= total +for i = 0, n_tags(!dustem_data)-1 do (*!fit_rchi2_weight).(i) /= total dustem_init,mode=use_mode diff --git a/src/idl/dustem_init.pro b/src/idl/dustem_init.pro index da6b0f1..84a3de7 100644 --- a/src/idl/dustem_init.pro +++ b/src/idl/dustem_init.pro @@ -211,13 +211,15 @@ if !run_pol then begin endif else tagnames=tag_names(*!dustem_data) -instr="defsysv, '!fit_rchi2_weight', {" +;instr="defsysv, '!fit_rchi2_weight', {" +instr="fit_rchi2_weight_str={" FOR i = 0, n_elements(tagnames)-1 DO BEGIN instr+=tagnames(i)+': rchi2_weight.'+tagnames(i) IF i NE n_elements(tagnames)-1 THEN instr+=',' ENDFOR instr+='}' toto=execute(instr) +defsysv, '!fit_rchi2_weight',ptr_new(fit_rchi2_weight_str) ;#########system variables necessary for cgwindow plotting########## defsysv, '!dustemcgwin_id', { $ ;IDs of windows to plot diff --git a/src/idl/dustem_mpfit_data.pro b/src/idl/dustem_mpfit_data.pro index efd2cc4..2fedefd 100755 --- a/src/idl/dustem_mpfit_data.pro +++ b/src/idl/dustem_mpfit_data.pro @@ -103,8 +103,8 @@ FOR j=0,n_elements(fields)-1 DO BEGIN ; 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) ) + if (*!fit_rchi2_weight).(i) gt 0 then begin + weight = sqrt(1./(*!fit_rchi2_weight).(i) ) endif else begin weight = 1.d99 endelse diff --git a/src/idl/dustem_mpfit_run.pro b/src/idl/dustem_mpfit_run.pro index 94db7e8..6d37535 100755 --- a/src/idl/dustem_mpfit_run.pro +++ b/src/idl/dustem_mpfit_run.pro @@ -102,7 +102,7 @@ FOR i_tag=0,count_data_tag-1 DO BEGIN 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 + 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 @@ -129,7 +129,7 @@ FOR i_tag=0,count_data_tag-1 DO BEGIN 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 + 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 @@ -169,7 +169,7 @@ FOR i_tag=0,count_data_tag-1 DO BEGIN 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 + 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 @@ -181,7 +181,7 @@ FOR i_tag=0,count_data_tag-1 DO BEGIN 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 + 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 @@ -206,7 +206,7 @@ FOR i_tag=0,count_data_tag-1 DO BEGIN 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 + 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 @@ -218,7 +218,7 @@ FOR i_tag=0,count_data_tag-1 DO 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 + 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 @@ -237,34 +237,34 @@ DOF = n_elements(p_min) total_chi2 = 0 total_rchi2 = 0 total_n = 0 -if !fit_rchi2_weight.ext ne 0 then begin +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 +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.qsed ne 0 then begin + 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 + if (*!fit_rchi2_weight).used ne 0 then begin total_n += n_used total_chi2 += chi2_used total_rchi2 += rchi2_used endif - if !fit_rchi2_weight.qext ne 0 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 + if (*!fit_rchi2_weight).uext ne 0 then begin total_n += n_uext total_chi2 += chi2_uext total_rchi2 += rchi2_uext @@ -273,16 +273,16 @@ if !run_pol then begin endif -total_wchi2 = chi2_ext * !fit_rchi2_weight.ext $ - + chi2_sed * !fit_rchi2_weight.sed +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 $ ; + 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 + + chi2_qsed * (*!fit_rchi2_weight).qsed $ + + chi2_used * (*!fit_rchi2_weight).used $ + + chi2_qext * (*!fit_rchi2_weight).qext $ + + chi2_uext * (*!fit_rchi2_weight).uext endif ;Using this expression for now for the weighting. @@ -296,8 +296,8 @@ 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 +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) 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. diff --git a/src/idl/dustem_plot_fit_extsed.pro b/src/idl/dustem_plot_fit_extsed.pro index 8621f68..c023ed0 100755 --- a/src/idl/dustem_plot_fit_extsed.pro +++ b/src/idl/dustem_plot_fit_extsed.pro @@ -128,7 +128,7 @@ dustem_ext = interpol(st.ext.ext_tot,st.ext.wav,(*!dustem_data.ext).wav) 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 +wrchi2_ext = rchi2_ext * (*!fit_rchi2_weight).ext print,"chi2 EXT = ",chi2_ext," rchi2 EXT = ",rchi2_ext;," weighted rchi2 EXT=",wrchi2_ext ;xr=[0.3,40] ;yr=[1.e-26,1.e-21] diff --git a/src/idl/dustem_set_data.pro b/src/idl/dustem_set_data.pro index fc891cb..bbd9858 100755 --- a/src/idl/dustem_set_data.pro +++ b/src/idl/dustem_set_data.pro @@ -313,30 +313,30 @@ if !run_pol then begin endif IF keyword_set(rchi2_weight) THEN BEGIN - !fit_rchi2_weight.sed = rchi2_weight.sed - !fit_rchi2_weight.ext = rchi2_weight.ext + (*!fit_rchi2_weight).sed = rchi2_weight.sed + (*!fit_rchi2_weight).ext = rchi2_weight.ext IF !run_pol THEN BEGIN - IF isa(qsed) THEN !fit_rchi2_weight.qsed = rchi2_weight.qsed - IF isa(used) THEN !fit_rchi2_weight.used = rchi2_weight.used + IF isa(qsed) THEN (*!fit_rchi2_weight).qsed = rchi2_weight.qsed + IF isa(used) THEN (*!fit_rchi2_weight).used = rchi2_weight.used - IF isa(qext) THEN !fit_rchi2_weight.qext = rchi2_weight.qext - IF isa(uext) THEN !fit_rchi2_weight.uext = rchi2_weight.uext + IF isa(qext) THEN (*!fit_rchi2_weight).qext = rchi2_weight.qext + IF isa(uext) THEN (*!fit_rchi2_weight).uext = rchi2_weight.uext ENDIF ENDIF ELSE BEGIN - !fit_rchi2_weight.sed=1. - !fit_rchi2_weight.ext=1. + (*!fit_rchi2_weight).sed=1. + (*!fit_rchi2_weight).ext=1. If !run_pol THEN BEGIN - IF isa(qsed) THEN !fit_rchi2_weight.qsed=1. - IF isa(used) THEN !fit_rchi2_weight.used=1. + IF isa(qsed) THEN (*!fit_rchi2_weight).qsed=1. + IF isa(used) THEN (*!fit_rchi2_weight).used=1. - IF isa(qext) THEN !fit_rchi2_weight.qext=1. - IF isa(uext) THEN !fit_rchi2_weight.uext=1. + IF isa(qext) THEN (*!fit_rchi2_weight).qext=1. + IF isa(uext) THEN (*!fit_rchi2_weight).uext=1. ENDIF diff --git a/src/idl/dustem_show_file_dependencies.pro b/src/idl/dustem_show_file_dependencies.pro index ac19f90..7b19eaf 100644 --- a/src/idl/dustem_show_file_dependencies.pro +++ b/src/idl/dustem_show_file_dependencies.pro @@ -52,14 +52,11 @@ dustem_run_example,'MC10' dustem_make_polarization_sed_example,model='G17_MODELA' dustem_fit_intensity_example,Nitermax=1,fits_save_and_restore='/tmp/mysavefile1.fits' dustem_fit_polarization_example,Nitermax=1,fits_save_and_restore='/tmp/mysavefile2.fits' - -dustem_fit_sed_ext_pol_example,Nitermax=1,fits_save_and_restore='/tmp/mysavefile2.fits' - -dustem_fit_ext_pol_example,Nitermax=1,fits_save_and_restore='/tmp/mysavefile2.fits' -dustem_fit_intensity_mbb_example,Nitermax=1 -dustem_fit_modified_bb_example,Nitermax=1 - -;dustem_fit_intensity_example +dustem_fit_ext_pol_example,Nitermax=1,fits_save_and_restore='/tmp/mysavefile3.fits' +dustem_fit_sed_ext_pol_example,Nitermax=1,fits_save_and_restore='/tmp/mysavefile4.fits' +dustem_fit_intensity_mbb_example,Nitermax=1,postscript='/tmp/dustemwrap_postcript_example.ps' +;The following is obsolete. No need to be there. +;dustem_fit_modified_bb_example,Nitermax=1 ; END LIST OF TOP-LEVEL ROUTINES help,/function,output=ff -- libgit2 0.21.2