replaced system variable fit_rchi2_weight by a pointer

Jean-Philippe Bernard
1 parent 7fb01a01
Showing 7 changed files with 49 additions and 48 deletions Show diff stats
src/idl/dustem_fit.pro
src/idl/dustem_init.pro
src/idl/dustem_mpfit_data.pro
src/idl/dustem_mpfit_run.pro
src/idl/dustem_plot_fit_extsed.pro
src/idl/dustem_set_data.pro
src/idl/dustem_show_file_dependencies.pro
@@ -71,17 +71,19 @@ defsysv, &#39;!dustem_data&#39;, {    $	;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
  
@@ -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                  
@@ -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
@@ -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,&quot;Nb of data points = &quot;, 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.
@@ -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]
@@ -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
  
@@ -52,14 +52,11 @@ dustem_run_example,&#39;MC10&#39;
 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