improved

Ilyes Choubani
1 parent cba7ef49
Showing 8 changed files with 43 additions and 48 deletions Show diff stats
src/idl/dustem_fit_sed_extsed_readme.pro
src/idl/dustem_fit_sed_polsed_readme.pro
src/idl/dustem_fit_sed_readme.pro
src/idl/dustem_init_plugins.pro
src/idl/dustem_plot_ext.pro
src/idl/dustem_set_params.pro
src/idl/dustem_write_all.pro
src/idl/dustem_write_all_web3p8.pro
@@ -93,22 +93,22 @@ st=dustem_set_data(ext=extsed,sed=spec,rchi2_weight=rchi2_weight,f_HI=f_HI)
 ;=== Set which parameters you want to fit
 pd = [ $
              'dustem_create_continuum_2',$
-             'dustem_create_powerlaw_1',$
+             ;'dustem_create_powerlaw_1',$
              '(*!dustem_params).gas.G0', $      ;G0
              '(*!dustem_params).grains(0).mdust_o_mh',$  ;PAH0 mass fraction
              '(*!dustem_params).grains(1).mdust_o_mh',$     ;PAH1 mass fraction
              '(*!dustem_params).grains(2).mdust_o_mh' $    ;amCBEx
               ]                   
  
-p_truth=[0.002,0.5,1.,7.8000E-04,7.8000E-04,1.6500E-04]    
+p_truth=[0.002,1.,7.8000E-04,7.8000E-04,1.6500E-04]    
  
 ;=== set initial value of parameters you want to fit
 ;iv=p_truth   ;exact solution
-iv=p_truth+[0.01,0.07,0.,3.e-5,1.e-5,-2.e-5]    ;shifted from solution
+iv=p_truth+[0.01,0.,3.e-5,1.e-5,-2.e-5]    ;shifted from solution
  
-ulimed=[0, 0 ,0   , 0   ,0   ,0 ]
-llimed=[1 ,1 ,1   , 1   ,1   ,1 ]
-llims =[0. ,0., 0.  ,0.   ,0.  ,0.]
+ulimed=[0, 0 ,0   , 0   ,0   ]
+llimed=[1 ,1 ,1   , 1   ,1   ]
+llims =[0. ,0., 0.  ,0. ,0.  ]
  
  
 ; 
@@ -124,9 +124,9 @@ dustem_init_parinfo,pd,iv,up_limited=ulimed,lo_limited=llimed,up_limits=ulims,lo
  
  
 ;=== RUN fit
-tol=1.e-5
-xtol=1.e-5
-Nitermax=6    ;maximum number of iteration. This is the criterium which will stop the fit procedure
+tol=1.e-12
+xtol=1.e-12
+Nitermax=3   ;maximum number of iteration. This is the criterium which will stop the fit procedure
 IF keyword_set(itermax) THEN Nitermax=itermax
  
 xrange=[1.,2.e4]
@@ -67,14 +67,14 @@ sed=dustem_initialize_sed(Nfilt)
 sed.filter=filters
 sed.wave=dustem_filter2wav(filters)
 sed.instru=dustem_filter2instru(filters)
-polsed=sed
-;polext=sed 
+polsed=sed ; comment if polext is to be used
+;polext=sed ;comment if polsed is to be used
  
-;stokes fake parameters
+;stokes fake parameters - comment if stokes extinction parameters are to be used
 qsed=sed
 used=sed
  
-; ;stokes extinction fake parameters
+; ;stokes extinction fake parameters - comment if stokes parameters are to be used
 ; qext=sed
 ; uext=sed
  
@@ -136,18 +136,19 @@ help,!run_pol
 sed.spec=dustem_compute_sed(p_truth,sst,_extra=extra)   ;sst is not set on purpose, for dustem_compute_sed to compute the sed using fortran code
 sed.error=(sed.spec*0.1);>0.1
  
+
 polsed.spec=dustem_compute_polsed(p_truth,sstp,_extra=extra,dustem_qsed,dustem_used)   ;sst is not set on purpose, for dustem_compute_sed to compute the sed using fortran code
 polsed.error=(polsed.spec*0.1)
  
  
-;creating fake polext
+; ;creating fake polext
+;   
+; polext.spec=dustem_compute_polext(p_truth,ssttp,dustem_qext,dustem_uext)
+; polext.error=(polext.spec*0.1)
  
-;polext.spec=dustem_compute_polext(p_truth,ssttp,dustem_qext,dustem_uext)
-;polext.error=(polext.spec*0.1)
-;  
  
  
-;SET qsed and used data
+;SET qsed and used data - comment if extinction polarization data is to be used
 qsed.spec=dustem_qsed
 qsed.error=qsed.spec*0.1
  
@@ -156,20 +157,20 @@ used.spec=dustem_used
 used.error=used.spec*0.1
  
  
-; ;SET qext and uext data
+; ;SET qext and uext data - comment if polarizarion data is to be used
 ; qext.spec=dustem_qext
 ; qext.error=qext.spec*0.1
 ; 
 ; 
 ; uext.spec=dustem_uext
 ; uext.error=uext.spec*0.1
-; 
+; ; 
  
 ;=== SET THE OBSERVATION STRUCTURE
  
 ind=where(sed.wave GE 0.1)
  
-st=dustem_set_data(polsed=polsed[ind],sed=sed[ind],qsed=qsed[ind],used=used[ind])
+st=dustem_set_data(polsed=polsed[ind],sed=sed[ind],qsed=qsed[ind],used=used[ind]) ; comment if second line is to be used
 ;st=dustem_set_data(sed=sed[ind],polext=polext[ind],qext=qext[ind],uext=uext[ind])
  
 ;=== RUN fit
@@ -202,7 +203,7 @@ message,&#39;Saved &#39;+file_out,/continue
 stop
  
 ;======================================
-;====You can exit IDL here and re-enter
+;====You can exit IDL here and re-enter   -   THIS PART APPEARS TO BE FAULTY. NEEDS INVESTIGATION 
 ;======================================
 file='/tmp/DUSTEM_polsed_fit_example.sav'
 dustem_restore_system_variables,file
@@ -228,7 +229,7 @@ xtit=textoidl(&#39;\lambda (\mum)&#39;)
 loadct,13
 IF keyword_set(postcript) THEN BEGIN
   set_plot,'PS'
-  ps_file=!dustem_wrap_soft_dir+'/Docs/Figures/'+'Last_dustem_spinning_fit.ps'
+  ps_file=!dustem_wrap_soft_dir+'/Docs/'+'Last_dustem_spinning_fit.ps';'/Docs/Figures/'+'Last_dustem_spinning_fit.ps'
   device,filename=ps_file,/color
 ENDIF
  
@@ -248,13 +248,13 @@ st=dustem_set_data(sed=spec)
 ;== SET THE FITTED PARAMETERS
 dustem_init_parinfo,pd,iv,up_limited=ulimed,lo_limited=llimed,up_limits=ulims,lo_limits=llims
  
+
 dustem_init_fixed_params,fpd,fiv
  
 dustem_init_plugins, pd ;initializing the scopes of the plugins to be read later
  
  
  
-
 ;=== RUN fit
 tol=1.e-14
 use_Nitermax=2               ;maximum number of iteration. This is the criterium which will stop the fit procedure
@@ -17,7 +17,7 @@ if fi eq &#39;dustem&#39; then begin
 fii=strtrim(strmid(pd(i),14,4),2)
  
 j=0L
-if i eq 0 then j+=1
+if i eq 0 then j+=2
  
 trl=fii eq strtrim(strmid(pd(i+j-1),14,4),2)
  
@@ -28,42 +28,38 @@ case fii of ;TESTING ON THE 4 FIRST LETTERS OF THE PLUGIN&#39;S NAME
 'cont': begin
  
 plugin_name='CONTINUUM' ; THIS STRING SHOULD MATCH THE NAME OF THE PLUGIN IN THE SCOPE LINE OF THE PLUGIN FUNCTION
-if test && (~trl) then varvar=create_struct(plugin_name,'ADD_SED') else begin
-if (~trl) then varvar=create_struct(varvar,plugin_name,'ADD_SED')
-endelse 
+if test && (~trl) then varvar=create_struct(plugin_name,'ADD_SED')
+if (~test) && (~trl) then varvar=create_struct(varvar,plugin_name,'ADD_SED')
+print, varvar
 end
  
 'sync': begin
  
 plugin_name='SYNCHROTRON'
-if test && (~trl) then varvar=create_struct(plugin_name,'ADD_SED+ADD_POLSED') else begin
-if (~trl) then varvar=create_struct(varvar,plugin_name,'ADD_SED+ADD_POLSED')
-endelse
+if test && (~trl) then varvar=create_struct(plugin_name,'ADD_SED+ADD_POLSED')
+if (~test) && (~trl) then varvar=create_struct(varvar,plugin_name,'ADD_SED+ADD_POLSED')
 end
  
 'free': begin
  
 plugin_name='FREE_FREE'
-if test && (~trl) then varvar=create_struct(plugin_name,'ADD_SED') else begin
-if (~trl) then varvar=create_struct(varvar,plugin_name,'ADD_SED')
-endelse
+if test && (~trl) then varvar=create_struct(plugin_name,'ADD_SED')
+if (~test) && (~trl) then varvar=create_struct(varvar,plugin_name,'ADD_SED')
 end
  
 'stok': begin
  
 plugin_name='STOKES'
-if test && (~trl) then varvar=create_struct(plugin_name,'MODIFY') else begin ; plugin that allows the computation of Q and U and fitting of PSI_dust and PSI_synch
-if (~trl) then varvar=create_struct(varvar,plugin_name,'MODIFY')
-endelse
+if test && (~trl) then varvar=create_struct(plugin_name,'MODIFY'); plugin that allows the computation of Q_sed and U_sed and fitting of PSI_dust and PSI_synch
+if (~test) && (~trl) then varvar=create_struct(varvar,plugin_name,'MODIFY')
 end
  
 'stex':begin
  
 print,'found'
 plugin_name='STEXT' ;'stokes'+extinction
-if test && (~trl) then varvar=create_struct(plugin_name,'MODIFY') else begin ; plugin that allows the computation of Q_ext = sigext_q and U_ext=sigext_u and fitting of PSI_ext
-if (~trl) then varvar=create_struct(varvar,plugin_name,'MODIFY')
-endelse
+if test && (~trl) then varvar=create_struct(plugin_name,'MODIFY') ; plugin that allows the computation of Q_ext = sigext_q and U_ext=sigext_u and fitting of PSI_ext
+if (~test) && (~trl) then varvar=create_struct(varvar,plugin_name,'MODIFY')
 print, varvar
 end
  
@@ -179,9 +179,9 @@ ENDIF ELSE BEGIN
 	cgoplot,xrange,(st.ext.ext_tot)/norm,psym=psym
  
  
-	; LINES OF CODE TO PLOT THE NEW ADDED PLUGINS TO THE EXTINCTION WILL HAVE TO BE ADDED HERE 
-	
-	cgoplot,xrange,(*(*!dustem_plugin).powerlaw)/norm,psym=psym,color='cadetblue'  
+; 	; LINES OF CODE TO PLOT THE NEW ADDED PLUGINS TO THE EXTINCTION WILL HAVE TO BE ADDED HERE 
+; 	
+; 	cgoplot,xrange,(*(*!dustem_plugin).powerlaw)/norm,psym=psym,color='cadetblue'  
  
  
  
@@ -214,7 +214,7 @@ ENDIF ELSE BEGIN
 	IF keyword_set(ps) AND NOT keyword_set(donotclose) THEN BEGIN
 		device,/close
 		set_plot,previous_device
-		;set_plot,'X'
+		set_plot,'X'
 	ENDIF
  
 ENDELSE
@@ -34,7 +34,6 @@ FOR i=0L,Nparams-1 DO BEGIN
         ;strr= strupcase(strmid(strmid((*(*!dustem_fit).param_descs)(i),0,strlen((*(*!dustem_fit).param_descs)(i))-2),14))
  
         k=where(tag_names(*!dustem_scope) eq strupcase(strr), countt)
-        
  
 ;     stop
  
@@ -52,7 +52,7 @@ IF !dustem_which EQ &#39;VERSTRAETE&#39; THEN BEGIN
 ENDIF
  
 ;IF getenv('DUSTEM_WHICH') EQ 'WEB3p8' THEN BEGIN
-IF !dustem_which EQ 'WEB3p8' THEN BEGIN    ; changed this as well 
+IF !dustem_which EQ 'WEB3p8' THEN BEGIN    
   dustem_write_all_web3p8,st,dir_out
   goto,sortie
 ENDIF
@@ -19,8 +19,7 @@ dustem_write_grain_web3p8,file_out,st.grains
  
 ;== GAS
 file_out=dir_out_dat+'GAS.DAT'
-print, tag_names(st)
-print, n_tags(st)
+
 dustem_write_gas,file_out,st.GAS
  
 ;stop