cleaning/improving/filling gaps

Ilyes Choubani
1 parent 8bab7d4b
Showing 9 changed files with 39 additions and 67 deletions Show diff stats
src/idl/dustem_activate_plugins.pro
src/idl/dustem_check_data.pro
src/idl/dustem_fit_sed_readme.pro
src/idl/dustem_init.pro
src/idl/dustem_mpfit_data.pro
src/idl/dustem_mpfit_run.pro
src/idl/dustem_restore_system_variables.pro
src/idl/dustem_save_system_variables.pro
src/idl/dustem_set_data.pro
@@ -69,7 +69,6 @@ FOR i=0L,n_elements(param_descs)-1 DO BEGIN
  
    IF parameter_type EQ 'PLUGIN' THEN BEGIN   ; Selecting the plugins
  
-        stop
         ftn = strmid(param_descs(i),0) ; String containing the name of the plugin and the keyword used (ie: dustem_create_continuum_2)    
         ii = strsplit(ftn,'_',count=countx) & ii = ii(countx-1)-1 ; Locating the last underscore to automate the extraction of the plugin's keyword
         ftn = strmid(ftn,0,ii) ; String containing the name of the plugin without the associated keyword
@@ -331,7 +331,7 @@ IF KEYWORD_SET(data_ext) THEN begin
     values = 0.5*atan(uext.values(uexpts),qext.values(qexpts))/degtorad
     polar_variance_iqu2ippsi,ext.values(sedpts),qexsed.values(qexpts),uext.values(uexpts),(data_ext(testx).sigextII)(extpts),(data_ext(testqex).sigextQQ)(qexpts),(data_ext(testuex).sigextUU)(uexpts),(data_ext(testx).sigextIQ)(extpts),(data_ext(testx).sigextIU)(extpts),(data_ext(testx).sigextQU)(extpts),variance_smallp,variance_psi,variance_largep
     sigma = sqrt(variance_psi)
-    psi_em={instru_names:instru_names,filt_names:filt_names,wav:wav,values:values,sigma:sigma}
+    psi_ext={instru_names:instru_names,filt_names:filt_names,wav:wav,values:values,sigma:sigma}
  
     goto, pexzone  ; the variances for polextPSI_EXT were already computed. 
  
@@ -334,11 +334,8 @@ dustem_init_parinfo,pd,iv,up_limited=ulimed,lo_limited=llimed,up_limits=ulims,lo
 ;== INITIALIZE ANY PLUGINS
 dustem_init_plugins, pd,fpd 
  
-stop
 dustem_init_fixed_params,fpd,fiv
  
-
-stop
 ;== RUN THE FIT
 tol=1.e-16
 use_Nitermax=5 ;maximum number of iterations. 
@@ -379,7 +376,7 @@ if keyword_set(restore) then begin
    file=restore
    dustem_restore_system_variables,file
 end
-   
+ 
 ;== PLOT THE FIT RESULTS RESTORED FROM .SAV FILE
 ; tit='Spectral Energy Distribution (Saved)'
 ; ytit=textoidl('I_\nu (MJy/sr) for N_H=10^{20} H/cm^2')
@@ -122,7 +122,7 @@ defsysv,&#39;!indef&#39;, exists=exists_indef
 IF NOT exists_indef THEN defsysv,'!indef',-32768.
  
 ;define !dustem_fit
-dustem_fit_st={data:ptr_new(),wavelength:ptr_new(), $
+dustem_fit_st={data:ptr_new(), $ ;,wavelength:ptr_new(), ;because wavelength arrays in the different datasets (data) can be different
                param_descs:ptr_new(),param_init_values:ptr_new(),param_func:ptr_new(), $
                fixed_param_descs:ptr_new(),fixed_param_init_values:ptr_new(), $
                chi2:0.D0,rchi2:0.D0, $
@@ -63,7 +63,7 @@ if !run_pol then begin
  
     ind_data_tag=dustem_data_tag(count=count_data_tag)
     tagnames=tag_names(!dustem_data)
-    testpol = tagnames EQ 'POLSED' or tagnames EQ 'POLEXT' or tagnames EQ 'POLFRAC' or tagnames EQ 'POLSED' or tagnames EQ 'PSI_EM' or tagnames EQ 'PSI_EXT'
+    testpol = tagnames EQ 'POLSED' or tagnames EQ 'POLEXT' or tagnames EQ 'POLFRAC' or tagnames EQ 'PSI_EM' or tagnames EQ 'PSI_EXT'
     data_ind=where(~testpol,ctestpol)
     if ctestpol ne 0 then begin
         match,ind_data_tag,data_ind,ind1,ind2 ;not using match2 because when matching indices, if the first is 0 it doesn't count which is wrong.
@@ -54,12 +54,6 @@ IF !run_ionfrac gt 0. THEN toto = dustem_create_ionfrac()
  
 dustem_activate_plugins,p_dim/(*(*!dustem_fit).param_init_values),st ;st is an output...
  
-;stop
-
-;st1 = dustem_run(p_dim) 
-    
-;stop
-
 st_model=dustem_read_all(!dustem_soft_dir)
  
  
@@ -125,7 +119,7 @@ FOR i_tag=0,count_data_tag-1 DO BEGIN
             ;I think this test isn't the correct one because when we want to show the SED data and not fit it we have an issue
             ;with the current hiding/showing formalism!!!!  
  
-	        if isa(!dustem_data.sed) then dustem_out = [ dustem_out, dustem_sed ]
+	        dustem_out = [ dustem_out, dustem_sed ]
  
 		END
  
@@ -149,23 +143,23 @@ FOR i_tag=0,count_data_tag-1 DO BEGIN
 	        	;print,"chi2 per wl: ",((*!dustem_data.ext).values-dustem_ext)^2/(*!dustem_data.ext).sigma^2
  
 			; Plot if needed
-			IF !dustem_show_plot NE 0 THEN BEGIN
-
-				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
-        		
-				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
+; 			IF !dustem_show_plot NE 0 THEN BEGIN
+; 
+; 				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
+;         		
+; 				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
  
 			;dustem_plot_ext,st,/print_Rv ;commented this recently
  
@@ -277,7 +271,7 @@ FOR i_tag=0,count_data_tag-1 DO BEGIN
 			    n_qsed = n_elements((*!dustem_data.qsed).values)
 			    rchi2_qsed  = chi2_qsed / n_qsed
 	        	wrchi2_qsed  = rchi2_qsed * !fit_rchi2_weight.qsed 
-	            if isa(!dustem_data.qsed) then dustem_out = [ dustem_out, dustem_qsed ] ; classic command
+	            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
  
@@ -290,7 +284,7 @@ FOR i_tag=0,count_data_tag-1 DO BEGIN
  			    n_used = n_elements((*!dustem_data.used).values)
  			    rchi2_used  = chi2_used / n_used
   	        	wrchi2_used  = rchi2_used * !fit_rchi2_weight.used
- 	            if isa(!dustem_data.used) then dustem_out = [ dustem_out, dustem_used] ; classic command
+ 	            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
  
@@ -312,9 +306,6 @@ FOR i_tag=0,count_data_tag-1 DO BEGIN
                    IF !run_pol and !run_lin THEN BEGIN
                        degtorad = !pi/180
                        if ~isa(dustem_used) or ~isa(dustem_qsed) then toto = dustem_compute_stokes(p_dim,st,dustem_qsed,dustem_used,Q_spec,U_spec,PSI_spec,dustem_psi_em)
-                       ;dustem_psi_em = 0.5*atan(dustem_used,dustem_qsed)/degtorad 
-                       
-                      
                    ENDIF		
     	END
  
@@ -50,16 +50,12 @@ restore,file,verbose=verbose
  
 defsysv,'!dustem_dat',saved_dustem_dat
 defsysv,'!dustem_data',saved_dustem_data
+defsysv,'!dustem_show',saved_dustem_show
 defsysv,'!DUSTEM_DO_CC',saved_DUSTEM_DO_CC
 defsysv,'!DUSTEM_F90_EXEC',saved_DUSTEM_F90_EXEC
 defsysv,'!dustem_filters',saved_dustem_filters
 defsysv,'!dustem_fit',saved_dustem_fit
 defsysv,'!DUSTEM_F_HI',saved_DUSTEM_F_HI
-;defsysv,'!DUSTEM_IDL_CONTINUUM',saved_DUSTEM_IDL_CONTINUUM
-;defsysv,'!DUSTEM_IDL_FREEFREE',saved_DUSTEM_IDL_FREEFREE
-;defsysv,'!DUSTEM_IDL_SYNCHROTRON',saved_DUSTEM_IDL_SYNCHROTRON
-
-;defsysv, '!DUSTEM_COMPOSITE_STELLAR_POPULATION',saved_DUSTEM_COMPOSITE_STELLAR_POPULATION
 defsysv,'!DUSTEM_INPUTS',saved_DUSTEM_INPUTS
 defsysv,'!DUSTEM_INSTRUMENT_DESCRIPTION',saved_DUSTEM_INSTRUMENT_DESCRIPTION
 defsysv,'!DUSTEM_NEVER_DO_CC',saved_DUSTEM_NEVER_DO_CC
@@ -82,8 +78,6 @@ defsysv,&#39;!RUN_RRF&#39;,saved_RUN_RRF
 defsysv,'!RUN_TLS',saved_RUN_TLS
 defsysv,'!RUN_UNIV',saved_RUN_UNIV
 defsysv,'!DUSTEM_REDSHIFT',saved_DUSTEM_REDSHIFT
-;defsysv,'!DUSTEM_ISRF_STAR_ADD',saved_DUSTEM_ISRF_STAR_ADD
-defsysv, '!DUSTEM_SCOPE',saved_DUSTEM_SCOPE
 defsysv, '!DUSTEM_PLUGIN',saved_DUSTEM_PLUGIN
  
 the_end:
@@ -80,14 +80,12 @@ ENDIF
  
 saved_dustem_dat=!dustem_dat
 saved_dustem_data=!dustem_data
+saved_dustem_show=!dustem_show
 saved_dustem_do_cc=!DUSTEM_DO_CC
 saved_dustem_f90_exec=!DUSTEM_F90_EXEC
 saved_dustem_filters=!dustem_filters
 saved_dustem_fit=!dustem_fit
 saved_dustem_f_hi=!DUSTEM_F_HI
-;saved_dustem_IDL_CONTINUUM=!DUSTEM_IDL_CONTINUUM ;Probably not needed - further testing required
-;saved_DUSTEM_IDL_FREEFREE=!DUSTEM_IDL_FREEFREE ;Probably not needed - further testing required
-;saved_DUSTEM_IDL_SYNCHROTRON=!DUSTEM_IDL_SYNCHROTRON ;Probably not needed - further testing required
 saved_DUSTEM_INPUTS=!DUSTEM_INPUTS
 saved_DUSTEM_INSTRUMENT_DESCRIPTION=!DUSTEM_INSTRUMENT_DESCRIPTION
 saved_DUSTEM_NEVER_DO_CC=!DUSTEM_NEVER_DO_CC
@@ -110,21 +108,17 @@ saved_RUN_RRF=!RUN_RRF
 saved_RUN_TLS=!RUN_TLS
 saved_RUN_UNIV=!RUN_UNIV
 saved_DUSTEM_REDSHIFT=!DUSTEM_REDSHIFT
-;saved_DUSTEM_COMPOSITE_STELLAR_POPULATION=!DUSTEM_COMPOSITE_STELLAR_POPULATION
-saved_DUSTEM_SCOPE=!DUSTEM_SCOPE
 saved_DUSTEM_PLUGIN=!DUSTEM_PLUGIN
  
 save,saved_dustem_dat,$
 	saved_dustem_data, $
+	saved_dustem_show, $
 	saved_dustem_fit, $
 	saved_dustem_do_cc, $
 	saved_dustem_f90_exec, $
     saved_dustem_filters, $
     saved_dustem_fit, $
     saved_dustem_f_hi, $
-    ;saved_dustem_IDL_CONTINUUM, $ ;Probably not needed - further testing required
-	;saved_DUSTEM_IDL_FREEFREE, $ ;Probably not needed - further testing required
-	;saved_DUSTEM_IDL_SYNCHROTRON, $ ;Probably not needed - further testing required
 	saved_DUSTEM_INPUTS, $
 	saved_DUSTEM_INSTRUMENT_DESCRIPTION, $
 	saved_DUSTEM_NEVER_DO_CC, $
@@ -147,9 +141,8 @@ save,saved_dustem_dat,$
 	saved_RUN_TLS, $
 	saved_RUN_UNIV, $
 	saved_DUSTEM_REDSHIFT, $
-        saved_DUSTEM_SCOPE, $
-        saved_DUSTEM_PLUGIN, $
-	;saved_DUSTEM_COMPOSITE_STELLAR_POPULATION, $
+    saved_DUSTEM_PLUGIN, $
+	
 	file=file
  
 the_end:
-FUNCTION dustem_set_data,m_fit=m_fit,m_show=m_show,x_fit=x_fit,x_show=x_show,sed=sed,ext=ext,polsed=polsed,polext=polext,used=used,qsed=qsed,polfrac=polfrac,psi_em=psi_em,psi_ext=psi_ext,Qext=Qext,Uext=Uext,rchi2_weight=rchi2_weight,f_HI=f_HI,help=help
+FUNCTION dustem_set_data,m_fit=m_fit,m_show=m_show,x_fit=x_fit,x_show=x_show,rchi2_weight=rchi2_weight,f_HI=f_HI,help=help
  
 ;+
 ; NAME:
@@ -53,25 +53,21 @@ IF keyword_set(help) THEN BEGIN
   goto,the_end
 ENDIF
  
-;Cleaninig data/showing structure 
+;Cleaninig fitting/showing structure(s) 
 for i=0L,n_tags(!dustem_data)-1 do begin
     !dustem_data.(i) = ptr_new()
     !dustem_show.(i) = ptr_new()
 endfor
  
-
-;I have decided to do the generation of data in check_data. 
-
+;#FITTING STRUCTURE
  
 xx=dustem_check_data(data_sed=m_fit,data_ext=x_fit,sed=sed,ext=ext,polext=polext,polsed=polsed,polfrac=polfrac,psi_em=psi_em,qsed=qsed,used=used,qext=qext,uext=uext,psi_ext=psi_ext)    
 ;outputs whatever the input structure has, into organized seperate structures formatted for the fitting process.
  
 ;FITTING STRUCTURE INITIALIZATION
-
 IF isa(sed) THEN !dustem_data.sed = ptr_new(sed)
 IF isa(ext) THEN !dustem_data.ext = ptr_new(ext)
  
- 
 ;####################
 ;OLD VARIABLE. KEPT IT.
 ; If f_HI is specified, multiply the data by f_HI
@@ -84,8 +80,6 @@ endif else f_HI = 1.
 defsysv,'!dustem_f_HI',f_HI
 ;####################
  
-;NOTA BENE: The user is bound to fit Q and U together.
-
 if !run_pol then begin
  
     if keyword_set(m_fit) then begin
@@ -139,8 +133,8 @@ if !run_pol then begin
     endif
  
     if keyword_set(x_fit) then begin
+    
         ;Two tests that determine what the user will eventually fit
-        
         testexstks = isa(qext) and isa(uext)
  
         IF ~testexstks then message, 'Extinction stokes parameters/Polarization data is not set correctly. Please check the input structure. Aborting...'
@@ -168,15 +162,15 @@ if !run_pol then begin
     endif
 endif
  
+;#SHOWING STRUCTURE
  
 if not keyword_set(m_show) then m_show=m_fit
  
 ;Setting of the structure that will be displayed
 yy=dustem_check_data(data_sed=m_show,data_ext=x_show,sed=sed,ext=ext,polext=polext,polsed=polsed,polfrac=polfrac,qsed=qsed,used=used,qext=qext,uext=uext,psi_em=psi_em,psi_ext=psi_ext)
-;stop
+
 IF isa(sed) THEN !dustem_show.sed = ptr_new(sed)
 IF isa(ext) THEN !dustem_show.ext = ptr_new(ext)
-
  
 ;####################
 ;OLD VARIABLE. KEPT IT.
@@ -303,6 +297,10 @@ ENDIF ELSE BEGIN
  
 ENDELSE
  
+;Setting the data and wavelength tags in '!dustem_fit'
+((*!dustem_fit).data)=ptr_new(!dustem_data)
+
+
 the_end: