General update

Ilyes Choubani
1 parent 18e4331f
Showing 11 changed files with 355 additions and 305 deletions Show diff stats
src/idl/dustem_activate_plugins.pro
src/idl/dustem_fit_sed_readme.pro
src/idl/dustem_init.pro
src/idl/dustem_init_plugins.pro
src/idl/dustem_mpfit_data.pro
src/idl/dustem_mpfit_run.pro
src/idl/dustem_plot_fit_sed.pro
src/idl/dustem_read_all_web3p8.pro
src/idl/dustem_set_data.pro
src/idl/dustemcgwin_dataset.pro
src/idl_misc/JPBLib_for_Dustemwrap/General/modify_extra.pro
@@ -82,15 +82,6 @@ FOR i=0L,n_elements(p_min)-1 DO BEGIN
  
         ;==============Assigning data to the initialized plugin data arrays==============
  
-        
-        print, 'ftn is:'
-        print, ftn
-        print, 'key is:'
-        print, index
-        print, 'val is:'
-        print, value
-        
-        
         if strmid(ftn,0,13) eq 'dustem_plugin' then begin
  
  
@@ -343,7 +343,7 @@ CASE use_model OF
       END
       'THEMIS':BEGIN
          pd = [ $
-             '(*!dustem_params).G0', $      ;G0
+             '(*!dustem_params).G0', $      h;G0
              '(*!dustem_params).grains(0).mdust_o_mh',$  
              '(*!dustem_params).grains(1).mdust_o_mh',$     
              '(*!dustem_params).grains(2).mdust_o_mh', $  
@@ -358,7 +358,7 @@ CASE use_model OF
  
 ENDCASE
  
-dustem_init,model=use_model,pol=polarization
+dustem_init,model=use_model;,pol=polarization
 !dustem_nocatch=1
  
 !dustem_verbose=1
@@ -382,8 +382,8 @@ ind=where(spec.instru EQ &#39;FIRAS&#39;,count)
 IF count NE 0 THEN spec(ind).sigmaII=(0.2*spec(ind).StokesI)^2
  
 ;=== SET THE OBSERVATION STRUCTURE
-st=dustem_set_data(sed=spec)
-
+st=dustem_set_data(st_sed_fit=spec);sed=spec)
+;stop
 ;== SET THE FITTED PARAMETERS
 dustem_init_parinfo,pd,iv,up_limited=ulimed,lo_limited=llimed,up_limits=ulims,lo_limits=llims
  
@@ -401,8 +401,8 @@ use_Nitermax=5 ;maximum number of iteration. This is the criterium which will st
 IF keyword_set(itermax) THEN use_Nitermax=itermax
  
 t1=systime(0,/sec)
-yr=[1e-4,100]
-xr=[1,6e4]
+yr=[1.00e-4,1.00E2]
+xr=[1.00E0,6.00e4]
 ;legend_xpos=0.6
 ;legend_ypos=0.8
 tit='Spectral Energy Distribution (Running)'
@@ -205,7 +205,7 @@ defsysv, &#39;!dustemcgwin_ncmds&#39;, {    $ ;Keeping track of number of commands run.
                             ext:     {pl:0, nrm:0},       $                         
                             polext:  {pl:0, nrm:0},      $                          
                             polsed:  {pl:0, nrm:0},      $                                      
-                            polfrac: {pl:0, nrm:0},     $                       
+                            polfrac: {pl:0},     $                       
                             psi_em:  {pl:0},      $                      
                             psi_ext: {pl:0},     $                           
                             qsed:    {pl:0, nrm:0},        $                    
 PRO dustem_init_plugins, pd, fpd
  
+;NB: 
+;#Fixed parameters still haven't been taken care of
+;#There is use of pointers when the need does not present itself. This implies modifying the rest of the code.
+
  
 ;INITIALIZING THE SCOPES DATA ARRAYS--------------------
  
@@ -79,7 +83,7 @@ IF Nplugins NE 0 THEN BEGIN
         message, 'varvar is:',/continue
         print, varvar
  
-        defsysv, '!dustem_scope', ptr_new(varvar)
+        defsysv, '!dustem_scope', ptr_new(varvar) 
  
    ENDIF ELSE BEGIN ;case of one plugin
  
@@ -58,7 +58,23 @@ IF keyword_set(Nitermax) THEN use_Nitermax=Nitermax
  
 ; Build array concatenating observational datas
 ; Order : that of !dustem_data (sed, ext, polext, polsed)
-ind_data_tag=dustem_data_tag(count=count_data_tag)
+
+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'
+    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.
+        ;I also only need the points that match.
+        ;testind1 = where(ind1 ne -1,ctind1)
+        ;if ctind1 ne 0 then ind_data_tag = ind_data_tag[testind1]
+        ind_data_tag = ind_data_tag[ind1]
+        count_data_tag = n_elements(ind_data_tag) 
+    endif
+endif else ind_data_tag=dustem_data_tag(count=count_data_tag)
+
  
 fields=['wav','values','sigma']
 FOR j=0,n_elements(fields)-1 DO BEGIN
@@ -104,6 +120,7 @@ p_min = dustem_mpfitfun(func_name,wav,values,sigma, $
                  parinfo=(*!dustem_parinfo), perror=perror, yfit=yfit, $
                  bestnorm=bestnorm, dof=dof,quiet=quiet, status=status, xtol=xtol, niter=niter_completed, ftol=use_tol, $
                  maxiter=use_Nitermax,gtol=gtol,errmsg=errmsg,functargs=_extra,nocatch=nocatch)
+                 
 message,errmsg,/info
 message,'mpfitfun Status='+strtrim(status,2),/info
 p_dim=p_min*(*(*!dustem_fit).param_init_values)
@@ -112,26 +112,12 @@ FOR i_tag=0,count_data_tag-1 DO BEGIN
 	        print,"Best SED chi2 obtained if f_HI = ",alpha, " or if we multiply model by ",beta
 	        (*!dustem_fit).chi2=chi2_sed
 	        (*!dustem_fit).rchi2=rchi2_sed
-			; Plot if needed
-	        IF !dustem_show_plot NE 0 THEN BEGIN
-	          
-	            ;window,win,title='DUSTEM WRAP (SED)' 
-	            
-	            ;cgwindow,WMULTI=[2,2,1],/CURRENT
-                
-                
-                
-                ;stop
-                ;dustemwrap_plott,p_dim,st,dustem_sed,dustem_polsed,dustem_qsed,dustem_polfrac,_extra=_extra
-	            ;dustem_plot_fit_sed,st,dustem_sed,_extra=_extra,res=p_min*(*(*!dustem_fit).param_init_values),chi2=(*!dustem_fit).chi2,rchi2=(*!dustem_fit).rchi2,fpol=fpol
-	       
-		 	ENDIF
-
-	        dustem_out = [ dustem_out, dustem_sed ]
+              
+	        if isa(!dustem_data.sed) then dustem_out = [ dustem_out, dustem_sed ]
  
 		END
  
-		'EXT' 	: BEGIN 
+		'EXT' 	: BEGIN ;STILL HAVEN'T DONE THE EXTINCTION PART
  
     		;ADDING PLUGIN TO SPECTRUM----------------
     		IF ptr_valid(*!dustem_plugin) THEN BEGIN
@@ -266,67 +252,65 @@ FOR i_tag=0,count_data_tag-1 DO BEGIN
  
 			END
  
- ;		'POLFRAC': BEGIN 
-;      		if !run_lin then begin
-;                  dustem_polsed = dustem_compute_polsed(p_dim,st) ; this quantity hasn't been computed yet that is why this line is present here.      
-;                  ind_filt=where((*!dustem_data.polfrac).filt_names NE 'SPECTRUM',count_filt) ;locating the filters of polfrac=smallp in SED xcat        
-;                  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    
-; ;                 
-; ;                 ;dustem_sed and dustem_polsed might not have the same length thus this block.
-; ;                 
-;                  If n_elements(dustem_sed) ne n_elements(dustem_polsed) then begin                            
-;                      if n_elements(dustem_sed) gt n_elements(dustem_polsed) then begin 
-;                          dustem_polsed_x = dustem_polsed
-;                          dustem_sed_x = dustem_polsed ;meaning dustem_sed needs to be modified
-;                          ;matching the filter data points
-;                          IF count_filt ne 0 then begin
-;                              for i=0L,count_filt-1 do begin    
-;                                  j=where((*!dustem_data.sed).filt_names EQ (*!dustem_data.polfrac).filt_names(ind_filt(i)),testfilt)
-;                                  if testfilt ne 0 then dustem_sed_x(ind_filt(i)) = dustem_sed(j(0))                                                 
-;                              endfor                                        
-;                          endif
-;                          ;matching the spectrum data points    
-;                          IF count_fspec ne 0 then begin
-;                              for i=0L,count_fspec-1 do begin
-;                                  j=where((*!dustem_data.sed).wav EQ (*!dustem_data.polfrac).wav(ind_fspec(i)),testspec)
-;                                  if testspec ne 0 then dustem_sed_x((ind_fspec(i))) = dustem_sed(j(0))                                    
-; ;                             endfor
-;                         endif     
-;                     endif ELSE begin
-;                         
-;                         dustem_polsed_x = dustem_sed
-;                         dustem_sed_x = dustem_sed ;meaning dustem_polsed needs to be modified
-;                         ;matching the filter data points
-;                         IF count_filt ne 0 then begin
-;                             for i=0L,count_filt-1 do begin    
-;                                 j=where((*!dustem_data.sed).filt_names EQ (*!dustem_data.polfrac).filt_names(ind_filt(i)),testfilt)
-;                                 if testfilt ne 0 then dustem_polsed_x(ind_filt(i)) = dustem_polsed(j(0))                                                 
-;                             endfor                                        
-;                         endif
-;                         ;matching the spectrum data points    
-;                         IF count_fspec ne 0 then begin
-;                             for i=0L,count_fspec-1 do begin
-;                                 j=where((*!dustem_data.sed).wav EQ (*!dustem_data.polfrac).wav(ind_fspec(i)),testspec)
-;                                 if testspec ne 0 then dustem_polsed_x((ind_fspec(i))) = dustem_polsed(j(0))                                    
-;                             endfor
-;                         endif
-;                     
-;                     ENDELSE
-;                 endif ELSE BEGIN
-;                     dustem_polsed_x = dustem_polsed 
-;                     dustem_sed_x = dustem_sed
-;                 ENDELSE
+;!!!!!!!!!!!!!YOU NEED TO RUN 'COMPUTE_SED' AGAIN because you have no idea if the sed block will be run.
+
+ 		'POLFRAC': BEGIN 
+      		if !run_lin then begin
+                  if ~isa(dustem_sed) then dustem_sed = dustem_compute_sed(p_dim,st) 
+                  dustem_polsed = dustem_compute_polsed(p_dim,st) ; this quantity hasn't been computed yet that is why this line is present here.      
+                  ind_filt=where((*!dustem_data.polfrac).filt_names NE 'SPECTRUM',count_filt) ;locating the filters of polfrac=smallp in SED xcat        
+                  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    
 ;                 
-;                 dustem_polfrac = dustem_polsed_x/dustem_sed_x 
-; 	        	chi2_polfrac   = total(((*!dustem_data.polfrac).values-dustem_polfrac)^2/(*!dustem_data.polfrac).sigma^2)
-; 				n_polfrac = n_elements((*!dustem_data.polfrac).values)
-; 	        	rchi2_polfrac  = chi2_polfrac / n_polfrac
-; 	        	wrchi2_polfrac  = rchi2_polfrac * !fit_rchi2_weight.polfrac
-; 	        	message,"chi2 POLFRAC = "+strtrim(chi2_polfrac,2)+"  rchi2 POLFRAC = "+strtrim(rchi2_polfrac,2),/continue    ;," weighted rchi2 POLFRAC=",wrchi2_polfrac
-;                 dustem_out = [ dustem_out, dustem_polfrac ]
-;                 ;fpol=dustem_polfrac ; problem because extra structure does not get updated so you actually need the 'extra' procedure you haven't finihed coding.
-;             endif 
- ;			END
+;                 ;dustem_sed and dustem_polsed might not have the same length thus this block.
+;                 
+                  If n_elements(dustem_sed) ne n_elements(dustem_polsed) then begin                            
+                      if n_elements(dustem_sed) gt n_elements(dustem_polsed) then begin 
+                          dustem_polsed_x = dustem_polsed
+                          dustem_sed_x = dustem_polsed ;meaning dustem_sed needs to be modified
+                          ;matching the filter data points
+                          IF count_filt ne 0 then begin
+                              for i=0L,count_filt-1 do begin    
+                                  j=where((*!dustem_data.sed).filt_names EQ (*!dustem_data.polfrac).filt_names(ind_filt(i)),testfilt)
+                                  if testfilt ne 0 then dustem_sed_x(ind_filt(i)) = dustem_sed(j(0))                                                 
+                              endfor                                        
+                          endif
+                          ;matching the spectrum data points    
+                          IF count_fspec ne 0 then begin
+                              for i=0L,count_fspec-1 do begin
+                                  j=where((*!dustem_data.sed).wav EQ (*!dustem_data.polfrac).wav(ind_fspec(i)),testspec)
+                                  if testspec ne 0 then dustem_sed_x((ind_fspec(i))) = dustem_sed(j(0))                                    
+                             endfor
+                        endif     
+                    endif ELSE begin
+                        
+                        dustem_polsed_x = dustem_sed
+                        dustem_sed_x = dustem_sed ;meaning dustem_polsed needs to be modified
+                        ;matching the filter data points
+                        IF count_filt ne 0 then begin
+                            for i=0L,count_filt-1 do begin    
+                                j=where((*!dustem_data.sed).filt_names EQ (*!dustem_data.polfrac).filt_names(ind_filt(i)),testfilt)
+                                if testfilt ne 0 then dustem_polsed_x(ind_filt(i)) = dustem_polsed(j(0))                                                 
+                            endfor                                        
+                        endif
+                        ;matching the spectrum data points    
+                        IF count_fspec ne 0 then begin
+                            for i=0L,count_fspec-1 do begin
+                                j=where((*!dustem_data.sed).wav EQ (*!dustem_data.polfrac).wav(ind_fspec(i)),testspec)
+                                if testspec ne 0 then dustem_polsed_x((ind_fspec(i))) = dustem_polsed(j(0))                                    
+                            endfor
+                        endif
+                    
+                    ENDELSE
+                endif ELSE BEGIN
+                    dustem_polsed_x = dustem_polsed 
+                    dustem_sed_x = dustem_sed
+                ENDELSE
+                
+                dustem_polfrac = dustem_polsed_x/dustem_sed_x 
+                	       
+            ENDIF
+            ;stop
+ 			END
  
  
 ;-----------------------------------------------------------------
@@ -335,40 +319,16 @@ FOR i_tag=0,count_data_tag-1 DO BEGIN
 		    ;PLUGIN IS ADDED INSIDE DUSTEM_COMPUTE_SED
  
  			IF !run_lin THEN BEGIN
-                ;make polsed only used for when polsed is fitted. Otherwise it's the fitting of Q_SED AND USED that fits the polarization data
+
 	        	dustem_polsed = dustem_compute_polsed(p_dim,st);dustem_compute_polsed(p_dim,st,cont=cont,freefree=freefree,synchrotron=synchrotron)
  
-	        	;chi2_polsed   = total(((*!dustem_data.polsed).values-dustem_polsed)^2/(*!dustem_data.polsed).sigma^2)
-				;n_polsed = n_elements((*!dustem_data.polsed).values)
-	        	;rchi2_polsed  = chi2_polsed / n_polsed
-	        	;wrchi2_polsed  = rchi2_polsed * !fit_rchi2_weight.polsed
-	        	;message,"chi2 POLSED = "+strtrim(chi2_polsed,2)+"  rchi2 POLSED = "+strtrim(rchi2_polsed,2),/continue;," weighted rchi2 POLSED=",wrchi2_polsed
-	        	;print,"chi2 per wl: ",((*!dustem_data.polsed).values-dustem_polsed)^2/(*!dustem_data.polsed).sigma^2
-	        	;print,"SNR per wl: ",(*!dustem_data.polsed).values/(*!dustem_data.polsed).sigma
-              
- 	            ;dustem_out = [ dustem_out, dustem_polsed ] ; classic command
-		        ;-----------------------------------------------
-	         
-				; Plot if needed
-; 				IF !dustem_show_plot NE 0 THEN BEGIN
-;  					win+=1
-; ; 					
-;  					xr=[10,8e3]
-;  					yr=[1e-34,1e-28]
-;  					
-;  					dustem_plot_polsed,st,p_dim,dustem_polsed,aligned=st_model.align.grains.aligned,win=win
-; ; 					dustem_plot_polar,st,p_dim,aligned=st_model.align.grains.aligned,/SED,ps=filename,win=win,cont=cont,freefree=freefree,synchrotron=synchrotron,xr=xr,yr=yr,/xsty,/donotclose,multi=1
-; ; 					dustem_plot_polar,st,p_dim,aligned=st_model.align.grains.aligned,/SED,ps=filename,win=win,cont=cont,freefree=freefree,synchrotron=synchrotron,xr=xr,yr=yr,/xsty,multi=2
-; 
-;                      					
-;  					
-; 				ENDIF
  			ENDIF
 		END
  
         'QSED': BEGIN 
  
 		        toto = dustem_compute_stokes(p_dim,st,dustem_qsed,dustem_used,Q_spec,U_spec) 
+		       
 		        ;taking care of the color correction for spass ; SHOULD NOT LEAVE THIS ON
 		        ;dustem_qsed(-1)=dustem_qsed(-2)
 ; 		        testspass1 = ((*!dustem_data.qsed).filt_names)(-1) eq 'SPASS1' 
@@ -378,118 +338,24 @@ 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 
-	            dustem_out = [ dustem_out, dustem_qsed ] ; classic command
+	            if isa(!dustem_data.qsed) then 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
  
-		        ;stop
-    		      
-		        ;stop
-; 		        if !dustem_show_plot NE 0 then begin
-; 		        fact = 1.00E-20;(3.e10/(1.e-4*(st.polsed.wav)))/1.d40
-;     	        fact1 = 1.00E-20;(3.e10/(1.e-4*((*!dustem_data.qsed).wav)))/1.d40 ; for data plotting   
-; 	            
-; ; 	            fact=1/(4.*!pi)/(3.e8/1.e-6/(st.polsed).wav)*1.00e17
-; ;               fact1=1.E-20
-; ; 	            
-; 		        win+=1
-;                 window  ,win ,title='DUSTEM WRAP (QSED)'
-; 		        titstq='Stokes Q Polarization Intensity'
-;                 ytitstq=textoidl('Q_\nu (W/m^2/Hz/sr) for N_H=10^{20} H/cm^2')
-;                 xtit=textoidl('\lambda (\mum)')
-;                 xr=[10,1e6]
-;                
-;                 ;normpol1=(*!dustem_data.qsed).values*fact1
-;                 yr=[1e-28,1e-17]
-;                 ;yr=[1e-34,1e-23]
-;                 indq=where(Q_spec lt 0, countq)
-;                 ;indq=where((*!dustem_data.qsed).values lt 0, countq)
-;                 
-;                 ylog=1
-;                 if countq ne 0 then begin
-; ;                 ylog=0
-; ;                 yr=[-1E-23,-1e-19];e-31 
-;                     Q_spec=-Q_spec
-;                     dustem_qsed =-dustem_qsed
-;                     ((*!dustem_data.qsed).values)=-((*!dustem_data.qsed).values)
-;                 endif
-;                 normpol=Q_spec*fact
-;                 normpol1=dustem_qsed*fact1
-;                 
-;                 
-;                 cgplot,st.polsed.wav,Q_spec*fact,xtit='',ytit=ytitstq,tit=titstq,ylog=ylog,/xlog,xr=xr,yr=yr,/ys,/xs,position=[0.17,0.35,0.93,0.93],xtickformat='(A1)'
-;                 cgoplot,st.polsed.wav,Q_spec*fact,color='black'
-;                 cgoplot,(*!dustem_data.qsed).wav,(*!dustem_data.qsed).values*fact1,color='Tomato',psym=16,symsize=1,thick=2
-;                 err_bar,(*!dustem_data.qsed).wav,(*!dustem_data.qsed).values*fact1,yrms=3.*((*!dustem_data.qsed).sigma)/2.*fact1,color=cgColor('Tomato')
-;                 cgoplot,(*!dustem_data.qsed).wav,dustem_qsed*fact1,psym=6,color='red',symsize=2
-;                 
-;                 cgplot,st.polsed.wav,Q_spec*fact/normpol,xtit=xtit,ytit='Normalized',tit='',/xlog,xr=xr,/ys,/xs,yr=[0,2],ylog=0,position=[0.17,0.14,0.93,0.35],/noerase,yticks=2,ymino=5,xticklen=0.1    
-;                 cgoplot,st.polsed.wav,Q_spec*fact/normpol,color='black'
-;                 cgoplot,(*!dustem_data.qsed).wav,(*!dustem_data.qsed).values*fact1/normpol1,color='Tomato',psym=16,symsize=1,thick=2
-;                 err_bar,(*!dustem_data.qsed).wav,(*!dustem_data.qsed).values*fact1/normpol1,yrms=3.*((*!dustem_data.qsed).sigma)/2.*fact1/normpol1,color=cgColor('Tomato')
-;                         
-;                 endif
 		    END   
  
         'USED': BEGIN
-		    
-                ;toto = dustem_compute_stokes(p_dim,st,dustem_qsed,dustem_used,Q_spec,U_spec)
-                ;stop
+        		    ;THIS IMPLIES THAT QSED HAS BEEN RAN BEFORE 
  			    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
- 	            dustem_out = [ dustem_out, dustem_used] ; classic command
+ 	            if isa(!dustem_data.used) then 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
-                dustemwrap_plot,p_dim,st,dustem_sed,dustem_qsed,dustem_used,dustem_polsed,_extra=_extra  
- 	            
- 	            ;if !dustem_show_plot ne 0 then begin
- 	            ;fact = 1.00E-20;(3.e10/(1.e-4*(st.polsed.wav)))/1.d40
-     	        ;fact1 =1.00E-20;(3.e10/(1.e-4*((*!dustem_data.used).wav)))/1.d40 ; for data plotting   
-	     
- 	            ;win+=1
-                ;window  ,win ,title='DUSTEM WRAP (USED)'                 
-                ;titstu='Stokes U Polarization Intensity'
-                ;ytitstu=textoidl('U_\nu (W/m^2/Hz/sr) for N_H=10^{20} H/cm^2')
-                ;xtit=textoidl('\lambda (\mum)')
-     	        ;xr=[10,1e6]
-     	        ;normpol2=U_spec*fact
-                ;normpol3=dustem_used*fact1
-                ;yr=[min(normpol2)-min(normpol2)/10,max(normpol2)+max(normpol2)*10]	
-                  
-                  ;conditions here are on computed spectra and they should be on the data itself.
-                  ;This even impairs the visualization of the data. Correct this!!!!!
-                  
-                  ;yr=[1e-28,1e-17];[1e-34,1e-23]
-                  ;indu=where(U_spec lt 0, countu)
-                  ;indu=where((*!dustem_data.used).values lt 0, countu)
-                  ;ylog=1
-                  ;if countu ne 0 then begin
-                  ;ylog=0
-                  ;yr=[-2e-31,5e-31]  
-                  ;endif 
-                  
-                ;cgplot,st.polsed.wav,U_spec*fact,xtit='',ytit=ytitstu,tit=titstu,ylog=ylog,/xlog,xr=xr,yr=yr,/ys,/xs,position=[0.17,0.35,0.93,0.93],xtickformat='(A1)',color='black'
-                ;cgoplot,st.polsed.wav,U_spec*fact,color='black'
-                ;cgoplot,(*!dustem_data.used).wav,(*!dustem_data.used).values*fact1,color='Teal',psym=16,symsize=1,thick=2
-                ;err_bar,(*!dustem_data.used).wav,(*!dustem_data.used).values*fact1,yrms=3.*((*!dustem_data.used).sigma)/2.*fact1,color=cgColor('Teal')
-                ;cgoplot,(*!dustem_data.used).wav,dustem_used*fact1,psym=6,color='red',symsize=2
-            
-                
-                
-                ;cgplot,st.polsed.wav,U_spec*fact/normpol2,xtit=xtit,ytit='Normalized',tit='',/xlog,xr=xr,/ys,/xs,yr=[0,2],ylog=0,position=[0.17,0.14,0.93,0.35],/noerase,yticks=2,ymino=5,xticklen=0.1    
-                ;cgoplot,st.polsed.wav,U_spec*fact/normpol2,color='black'
-                ;cgoplot,(*!dustem_data.used).wav,(*!dustem_data.used).values*fact1/normpol3,color='Teal',psym=16,symsize=1,thick=2
-                ;err_bar,(*!dustem_data.used).wav,(*!dustem_data.used).values*fact1/normpol3,yrms=3.*((*!dustem_data.used).sigma)/2.*fact1/normpol3,color=cgColor('Teal')
-                ;                
-                ;cgwindow, 'plot', 
-                	 
- 	            ;endif
+                dustemwrap_plot,p_dim,st,dustem_sed,dustem_qsed,dustem_used,dustem_polsed,dustem_polfrac,_extra=_extra  
+ 	
   		END
-;  		
-;  		
- 		
  
   		'QEXT': BEGIN ; Q - EXTINCTION CROSS SECTION
  
@@ -603,41 +469,27 @@ if !fit_rchi2_weight.sed ne 0 then begin
 	total_rchi2 += rchi2_sed
 endif
 if !run_pol then begin
-; if !fit_rchi2_weight.polext ne 0 then begin
-; 	total_n += n_polext
-; 	total_chi2 += chi2_polext
-; 	total_rchi2 += rchi2_polext
-; endif
-; if !fit_rchi2_weight.polsed ne 0 then begin
-; 	total_n += n_polsed 
-; 	total_chi2 += chi2_polsed 
-; 	total_rchi2 += rchi2_polsed 
-; endif
-; if !fit_rchi2_weight.polfrac ne 0 then begin
-; 	total_n += n_polfrac
-; 	total_chi2 += chi2_polfrac
-; 	total_rchi2 += rchi2_polfrac
-; endif
-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
-	total_n += n_used
-	total_chi2 += chi2_used
-	total_rchi2 += rchi2_used
-endif
-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
-	total_n += n_uext
-	total_chi2 += chi2_uext
-	total_rchi2 += rchi2_uext
-endif
+
+    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
+    	total_n += n_used
+    	total_chi2 += chi2_used
+    	total_rchi2 += rchi2_used
+    endif
+    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
+    	total_n += n_uext
+    	total_chi2 += chi2_uext
+    	total_rchi2 += rchi2_uext
+    endif
  
 endif
  
@@ -74,7 +74,7 @@ fact = 1.e4*(*!dustem_HCD)/(4.*!pi)/(3.e8/1.e-6/st.sed.wav)*1.e20/1.e7
 use_col_data_filt='blue'
  
 ;use_col_sed_spec='red'
-use_col_sed_spec='grey'
+use_col_sed_spec='Powder Blue'
  
 IF not keyword_set(col_sed) THEN BEGIN
   ;use_col_sed_filt=250     ;red
@@ -150,7 +150,8 @@ if keyword_set(fpol) then begin
       IF !d.name NE 'PS' THEN cgDisplay, 600, 500
       cgplot,(*!dustem_data.sed).wav,(*!dustem_data.sed).values/norm,/nodata,xtit='',ytit=ytit,tit='',/ylog,/xlog,/ys,/xs,position=[0.12,0.25,0.96,0.76],xtickformat='(A1)',_extra=_extra,charsize=1.3,/noerase    
     endif
-endif ELSE cgplot,(*!dustem_data.sed).wav,(*!dustem_data.sed).values/norm,/nodata,xtit='',ytit=ytit,tit=title,/ylog,/xlog,/ys,/xs,position=[0.12,0.35,0.96,0.90],xtickformat='(A1)',_extra=_extra,charsize=1.3
+endif ELSE cgplot,(*!dustem_data.sed).wav,(*!dustem_data.sed).values/norm,/nodata,ytit=ytit,tit=title,/ylog,/xlog,/ys,/xs,position=[0.12,0.35,0.96,0.90],charsize=1.3,_extra=_extra,xtickformat='(A1)'
+;stop  
 ;###############################
  
 ;this is where you will add the normalized sed
@@ -167,7 +168,7 @@ IF count_spec NE 0 THEN BEGIN
   xx=((*!dustem_data.sed).wav)[ind_spec]
   yy=((*!dustem_data.sed).values)[ind_spec]/norm[ind_spec]
   rms=3.*((*!dustem_data.sed).sigma)[ind_spec]/2./norm[ind_spec]
-  cgoplot,xx,yy,psym=8,syms=0.5,color=use_col_sed_spec
+  cgoplot,xx,yy,psym=8,color=use_col_sed_spec
   IF not keyword_set(no_spec_error) THEN BEGIN
     cgerrplot,xx,yy-rms,yy+rms,color=use_col_sed_spec
     ;err_bar,((*!dustem_data.sed).wav)(ind_spec),((*!dustem_data.sed).values)(ind_spec)/norm(ind_spec),yrms=3.*((*!dustem_data.sed).sigma)(ind_spec)/2./norm(ind_spec)
@@ -183,7 +184,7 @@ IF count_filt NE 0 THEN BEGIN
   ;err_bar,((*!dustem_data.sed).wav)(ind_filt),((*!dustem_data.sed).values)(ind_filt)/norm(ind_filt),yrms=3.*((*!dustem_data.sed).sigma)(ind_filt)/2./norm(ind_filt),color=use_col_data_filt
   cgerrplot,xx,yy-rms,yy+rms,color='Dodger Blue';use_col_data_filt
 ENDIF
-;=== Plot the computed SED - color corrected interpolates
+;=== Plot the computed SED
 IF count_filt NE 0 THEN BEGIN
   plotsym,8
   xx=((*!dustem_data.sed).wav)[ind_filt]
@@ -194,7 +195,7 @@ IF count_spec NE 0 THEN BEGIN
   plotsym,0
   xx=((*!dustem_data.sed).wav)[ind_spec]
   yy=sed[ind_spec]/norm[ind_spec]
-  cgoplot,xx,yy,color=use_col_sed_filt,psym=6,syms=2
+  cgoplot,xx,yy,color='Indian Red',psym=7,syms=1
 ENDIF
  
  
@@ -50,7 +50,7 @@ file=dir_in_dat+(*!dustem_inputs).grain
 ;st_grains=dustem_read_grain(file,silent=silent)
 st_grains=dustem_read_grain(file,silent=silent,key_str=key_str,G0=G0)
 Ngrains=n_elements(st_grains)
-
+    
 file=dir_in_dat+'ISRF.DAT'
 st_isrf=dustem_read_isrf(file,silent=silent,Nisrf=Nisrf)
  
@@ -63,7 +63,7 @@ for i=0L,n_tags(!dustem_data)-1 do begin
 endfor
  
 xx=dustem_check_data(data_sed=st_sed_fit,data_ext=st_ext_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)    
-
+;stop
 ;outputs whatever the input structure has, into organized seperate structures formatted for the fitting process.
  
 ;If SED/EXT data is present in the structure it is fitted automatically. 
@@ -91,42 +91,63 @@ if !run_pol then begin
  
         teststks = isa(qsed) and isa(used)
  
-        IF (not teststks)  then message, 'Stokes parameters are not set correctly. Please check the input structure. Aborting...'
+        IF ~teststks  then message, 'Stokes parameters are not set correctly. Please check the input structure. Aborting...'
  
         ;EMISSION 
-        if teststks then nowfitting = 'QSED and USED' 
  
-        message, 'Polarization component(s) to fit: '+nowfitting ,/continue
+        message, 'Polarization component(s) to fit: QSED and USED' ,/continue
  
         IF isa(qsed) THEN BEGIN        
             !dustem_data.qsed = ptr_new(qsed)
             if keyword_set(f_HI) then if f_HI gt 0 then $
           	(*!dustem_data.qsed).values *= f_HI
         ENDIF    
+        
         if isa(used) then BEGIN
             !dustem_data.used = ptr_new(used)
             if keyword_set(f_HI) then if f_HI gt 0 then $
           	(*!dustem_data.used).values *= f_HI
         ENDIF
+        
+        if isa(polsed) then BEGIN
+            
+            !dustem_data.polsed = ptr_new(polsed)
+            if keyword_set(f_HI) then if f_HI gt 0 then $
+          	(*!dustem_data.polsed).values *= f_HI
+            
+        ENDIF 
+        
+        if isa(polfrac) then BEGIN
+            
+            !dustem_data.polfrac = ptr_new(polfrac)
+            if keyword_set(f_HI) then if f_HI gt 0 then $
+          	(*!dustem_data.polfrac).values *= f_HI
+        
+        ENDIF 
+        
+        
+        if isa(psi_em) then BEGIN
+            
+            !dustem_data.psi_em = ptr_new(psi_em)
+            if keyword_set(f_HI) then if f_HI gt 0 then $
+          	(*!dustem_data.psi_em).values *= f_HI
+        
+        ENDIF 
+        
+        
     endif
  
     if keyword_set(st_ext_fit) then begin
         ;Two tests that determine what the user will eventually fit
  
         testexstks = isa(qext) and isa(uext)
-        testpexstks = (isa(qext) and isa(uext) and isa(polext))
-        
-        IF (not testexstks and not isa(polext)) then message, 'Extinction stokes parameters/Polarization data is not set correctly. Please check the input structure. Aborting...'
-        
-        IF testpexstks THEN $
-        message, 'Presence of both extinction polarization and stokes parameter values. Fitting structure is not correctly filtered. Aborting...'
  
+        IF ~testexstks then message, 'Extinction stokes parameters/Polarization data is not set correctly. Please check the input structure. Aborting...'
+         
         ;EXTINCTION - allowed for the default fitting of both because I do not know if the same formalism should be applied to the extinction data because it also implies deviding extinction data into two components. 
-        if testexstks then nowfitting = 'QEXT and UEXT' else begin ;Serkowski is default
-            if isa(polext) then nowfitting = 'POLEXT'      
-        endelse
-        
-        message, 'Extinction polarization component(s) to fit: '+nowfitting ,/continue
+       
+        message, 'Extinction polarization component(s) to fit: QEXT and UEXT',/continue
+       
         if isa(polext) then BEGIN        
             !dustem_data.polext = ptr_new(polext)
             if keyword_set(f_HI) then if f_HI gt 0 then $
@@ -142,6 +163,7 @@ if !run_pol then begin
             if keyword_set(f_HI) then if f_HI gt 0 then $
            	(*!dustem_data.uext).values *= f_HI
         ENDIF
+        
     endif
 endif
  
@@ -149,8 +171,8 @@ endif
 if not keyword_set(st_sed_show) then st_sed_show=st_sed_fit
  
 ;Setting of the structure that will be displayed
-yy=dustem_check_data(data_sed=st_sed_show,data_ext=st_ext_show,sed=sed,ext=ext,polext=polext,polsed=polsed,polfrac=polfrac,qsed=qsed,used=used,qext=qext,uext=uext)
-
+yy=dustem_check_data(data_sed=st_sed_show,data_ext=st_ext_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)
  
@@ -173,16 +195,16 @@ if !run_pol then begin
  
         teststks = isa(qsed) and isa(used)
  
-        IF (not teststks and not isa(polsed)) or (not teststks and not isa(polfrac)) then message, 'Stokes parameters/Polarization data is not set correctly. Please check the input structure. Aborting...'
+        IF (~teststks and ~isa(polsed)) or (~teststks and ~isa(polfrac)) then message, 'Stokes parameters/Polarization showing data is not set correctly. Please check the input structure. Aborting...'
  
-        if teststks and isa(polsed) and not (isa(polfrac)) and not(isa(psi_em)) then nowshowing = 'QSED, USED and POLSED'
-        if teststks and isa(polsed) and not(isa(polfrac)) and isa(psi_em) then nowshowing = 'QSED, USED, POLSED and PSI_EM'
-        if teststks and not(isa(polsed)) and isa(polfrac) and not(isa(psi_em)) then nowshowing = 'QSED, USED and POLFRAC'
-        if teststks and not(isa(polsed)) and isa(polfrac) and isa(psi_em) then nowshowing = 'QSED, USED, POLFRAC and PSI_EM'
-        if teststks and not(isa(polsed)) and not(isa(polfrac)) and not(isa(psi_em)) then nowshowing = 'QSED and USED'
-        if teststks and isa(polsed) and isa(polfrac) and not(isa(psi_em)) then nowshowing = 'QSED, USED, POLSED and POLFRAC'
+        if teststks and isa(polsed) and ~(isa(polfrac)) and ~(isa(psi_em)) then nowshowing = 'QSED, USED and POLSED'
+        if teststks and isa(polsed) and ~(isa(polfrac)) and isa(psi_em) then nowshowing = 'QSED, USED, POLSED and PSI_EM'
+        if teststks and ~(isa(polsed)) and isa(polfrac) and ~(isa(psi_em)) then nowshowing = 'QSED, USED and POLFRAC'
+        if teststks and ~(isa(polsed)) and isa(polfrac) and isa(psi_em) then nowshowing = 'QSED, USED, POLFRAC and PSI_EM'
+        if teststks and ~(isa(polsed)) and ~(isa(polfrac)) and ~(isa(psi_em)) then nowshowing = 'QSED and USED'
+        if teststks and isa(polsed) and isa(polfrac) and ~(isa(psi_em)) then nowshowing = 'QSED, USED, POLSED and POLFRAC'
         if teststks and isa(polsed) and isa(polfrac) and isa(psi_em) then nowshowing = 'QSED, USED, POLFRAC, POLSED and PSI_EM'
-        if teststks and not(isa(polsed)) and not(isa(polfrac)) and isa(psi_em) then nowshowing = 'QSED, USED and PSI_EM'
+        if teststks and ~(isa(polsed)) and ~(isa(polfrac)) and isa(psi_em) then nowshowing = 'QSED, USED and PSI_EM'
  
  
         message, 'Polarization component(s) to show: '+nowshowing ,/continue
@@ -217,12 +239,12 @@ if !run_pol then begin
     if isa(st_ext_show) then begin
         teststks = isa(qext) and isa(uext)
  
-        IF (not teststks and not isa(polext)) then message, 'Extinction stokes parameters/Polarization data is not set correctly. Please check the input structure. Aborting...'
+        IF (~teststks) then message, 'Extinction stokes parameters/Polarization data is not set correctly. Please check the input structure. Aborting...'
  
-        if teststks and isa(polext) and not(isa(psi_ext)) then nowshowing = 'QEXT, UEXT and POLEXT'
+        if teststks and isa(polext) and ~(isa(psi_ext)) then nowshowing = 'QEXT, UEXT and POLEXT'
         if teststks and isa(polext) and isa(psi_ext) then nowshowing = 'QEXT, UEXT, POLEXT and PSI_EXT'
-        if teststks and not(isa(polext)) and isa(psi_ext) then nowshowing = 'QEXT,UEXT and PSI_EXT' 
-        if teststks and not(isa(polext)) and not(isa(psi_ext)) then nowshowing = 'QEXT and UEXT'
+        if teststks and ~(isa(polext)) and isa(psi_ext) then nowshowing = 'QEXT,UEXT and PSI_EXT' 
+        if teststks and ~(isa(polext)) and ~(isa(psi_ext)) then nowshowing = 'QEXT and UEXT'
  
         message, 'Polarization component(s) to show: '+nowshowing ,/continue
  
@@ -255,12 +277,10 @@ IF keyword_set(rchi2_weight) THEN BEGIN
 	!fit_rchi2_weight.ext     = rchi2_weight.ext
  
 	IF !run_pol THEN BEGIN
-    	;IF isa(polsed) THEN !fit_rchi2_weight.polsed  = rchi2_weight.polsed
-    	;IF isa(polfrac) THEN !fit_rchi2_weight.polfrac = rchi2_weight.polfrac
+    	
     	IF isa(qsed) THEN !fit_rchi2_weight.qsed  = rchi2_weight.qsed
     	IF isa(used) THEN !fit_rchi2_weight.used  = rchi2_weight.used
  
-    	IF isa(polext) THEN !fit_rchi2_weight.polext  = rchi2_weight.polext
     	IF isa(qext) THEN !fit_rchi2_weight.qext  = rchi2_weight.qext
     	IF isa(uext) THEN !fit_rchi2_weight.uext  = rchi2_weight.uext
     ENDIF
@@ -271,13 +291,13 @@ ENDIF ELSE BEGIN
     !fit_rchi2_weight.ext=1.
  
     If !run_pol THEN BEGIN
-        ;IF isa(polsed) THEN !fit_rchi2_weight.polsed=1.
-        ;IF isa(polfrac) THEN !fit_rchi2_weight.polfrac=1.
+        
         IF isa(qsed) THEN !fit_rchi2_weight.qsed=1.
         IF isa(used) THEN !fit_rchi2_weight.used=1.
-        IF isa(polext) THEN !fit_rchi2_weight.polext=1.
+        
         IF isa(qext) THEN !fit_rchi2_weight.qext=1.
         IF isa(uext) THEN !fit_rchi2_weight.uext=1.
+    
     ENDIF
  
 ENDELSE
@@ -44,10 +44,14 @@ if !run_pol then begin
  
     specpol = st.polsed.em_tot * fact
  
+;    stop
     ;Instead of matching the non-zero values in spec and specpol (dustem outputs), 
     ;I will only test on specpol as spec does not contain any null values FOR NOW. 
     ;If this happens to change, these lines will have to be adjusted.
  
+    ;CHECKING IF THIS IS CORRECT FOR EACH OF THE SPECIES SPECTRA = ? 
+    
+    
     specpolfrac = specpol *0.
     psi_em = specpol*0.
     indpolsed = where(specpol ne 0, ct_nullpolsed) ;I am not using the counter for now as it serves me no direct purpose.
@@ -72,8 +76,12 @@ IF scopes[0] NE &#39;NONE&#39; THEN BEGIN
  
                 specq = (*(*!dustem_plugin).(i))[*,1]
                 specu = (*(*!dustem_plugin).(i))[*,2]
-                ;specpol = sqrt(1.0E10*specq^2+1.0E10*specu^2) ;sqrt(((*(*!dustem_plugin).(i))[*,1])^2+(*(*!dustem_plugin).(i))[*,2]^2)
+                
+                ;STILL LACKING A WORKAROUND. 
                 specpol = (specq*1.0E10)^2+(specu*1.0E10)^2 & specpol=1.0E-10*sqrt(specpol) ;for some reason cgwindow couldn't handle these mathematical operations. It seems this is caused by values that are too close to zero
+                
+                ;specpol = (specq^2 + specu^2) & specpol = sqrt(specpol)
+                
                 specpolfrac[indpolsed] = specpol[indpolsed]/spec[indpolsed]
                 psi_em[indpolsed] = 0.5*atan(specu[indpolsed],specq[indpolsed])/degtorad
                 ;stop
@@ -293,7 +301,9 @@ if keyword_set(dataset) then begin
  
         end
  
-        'POLSED': begin
+        'POLSED': begin ;#NB: There seems to be an issue with the display of the dust species. 'Carbonaceous grains' 's emission exceeds the total emission.
+                        ;MABE AN ISSUE WITH HOW 'SPECPOL' refreshes?
+                        ;PROCEEDING TO POLFRAC  
  
             idx_filt=where((*!dustem_data.polsed).filt_names NE 'SPECTRUM',ct_filt)
             idx_spec=where((*!dustem_data.polsed).filt_names EQ 'SPECTRUM',ct_spec)
@@ -426,7 +436,7 @@ if keyword_set(dataset) then begin
                         ;Locating all the hidden data points (spectrum+filter)
                         match2,((*!dustem_data.polsed).wav),((*!dustem_show.polsed).wav),fit_polsedpts,show_polsedpts ;only show_polsedpts is needed
                         idx_rmv_polsed=where(show_polsedpts eq -1, ct_hdnpts) ; indices of the points to hide
-                        
+                        ;stop
                         if ct_hdnpts ne 0 then begin ;Hidden data points are present 
                             ;stop
                             ;Locating the hidden spectrum and filter data points                                        
@@ -442,7 +452,7 @@ if keyword_set(dataset) then begin
                                 cgerrplot,(((*!dustem_show.polsed).wav)[idx_rmv_polsed])(idx_spec_hdn),(((*!dustem_show.polsed).values)[idx_rmv_polsed])(idx_spec_hdn)-rms,(((*!dustem_show.polsed).values)[idx_rmv_polsed])(idx_spec_hdn)+rms,color='Black'
                             endif                                            
  
-                            if ct_filt_hdn then begin
+                            if ct_filt_hdn ne 0 then begin
                             ;Plotting of hidden filter data points 
                                 cgplot,(((*!dustem_show.polsed).wav)[idx_rmv_polsed])(idx_filt_hdn),(((*!dustem_show.polsed).values)[idx_rmv_polsed])(idx_filt_hdn),pos=position,/ylog,/xlog,/ys,xs=9,noerase=1,charsize=1.15,xtickformat='(A1)',color='Black',xr=xr,yr=yr,psym=8,syms=0.8                                           
  
@@ -466,7 +476,161 @@ if keyword_set(dataset) then begin
  
         end
  
-        'POLFRAC': begin
+        'POLFRAC': begin ; YOU NEED TO CODE THIS PART - TOU STILL HAVE'NT FINISHED. 
+            idx_filt=where((*!dustem_data.polfrac).filt_names NE 'SPECTRUM',ct_filt)
+            idx_spec=where((*!dustem_data.polfrac).filt_names EQ 'SPECTRUM',ct_spec)
+            
+            ;#1) get the plotting keywords (pertaining to each data set) @here when the _extra structure is ready
+            
+            if keyword_set(nodata) then begin ;when the data is not present
+            
+                    xtit = ''
+                    cgplot,wavs,wavs*0.,/nodata,/xlog,/ys,xs=1,pos=position,noerase=1,charsize=1.15,xtickformat='(A1)',color='Powder Blue',xr=xr,xtit='',yr=[1.00E-03,50.00],/ylog,psym=8,syms=0.8  
+                    xyouts,pospltxt[0],pospltxt[1],textoidl('p_{\nu} (%)'),color=0,/normal,charsize=1.1
+                 
+            endif else begin ;when the data is present ; normal plot     
+                        
+                    if keyword_set(refresh) then begin ;The data points in the plot are being refreshed
+                        
+                        ;##################TAKEN FROM SED
+                        
+                        ;Plotting of the spectra of the dust species
+                        
+                        ;I THINK I WILL HAVE AN ERROR HERE BECAUSE I WILL HAVE TO DIVIDE TWO ARRAYS FOR EACH DUST SPECIES. 
+                        ; (st.polsed.(i+1)/st.sed.(i+1))*fact
+                        ;There is an issue with the fact that there are grain species that are not polarized
+                        ;test and set a counter? YES FOR NOW 
+                        ;
+                        
+                        FOR i=0L,Ngrains-1 DO BEGIN
+                            
+                            testneq = where(st.polsed.(i+1) ne 0, countneq)
+                            
+                            vecfin = st.polsed.(i+1)*0.0
+                            ;SET the condition later
+                            if countneq ne 0 then vecfin[testneq] = (st.polsed.(i+1))[testneq]/(st.sed.(i+1))[testneq]
+                             
+                            cgoplot,st.polsed.wav,vecfin*100,pos=position,noerase=1,color=use_cols[i],xr=xr,/xlog,/ylog,ytickformat='(A1)',yr=[1.00E-03,50.00]
+                       
+                        ENDFOR
+                       
+                        ;Plotting of the plugins. 
+                        
+                        ;NB: YOU HAVE YET TO TEST THIS BLOCK. DON'T FORGET TO ADD A POLARIZATION FRACTION TO FF for instance. 
+                        ; The thing is you've already had errors of that kind. 
+                        
+                        ;PLUGINS SEEM TO HAVE ZERO NULL VALUES. THIS IS GOOD NEWS. LESS TESTING.
+                        ;BUT if you get errors you know it has to be the division of the two arrays.
+                        ;THE CHANGES SEEM LIKE THEY'RE GONNA WORK.
+                        
+                        
+                        ;MISSING KEYWORDS? 
+                        for i=0L,n_plgns-1 do begin
+                            IF total(strsplit((*(*!dustem_scope).(i)),'+',/extract) EQ 'ADD_POLSED') THEN begin
+                                
+                                cgoplot,st.polsed.wav,sqrt(((*(*!dustem_plugin).(i))[*,1])^2+((*(*!dustem_plugin).(i))[*,2])^2)/((*(*!dustem_plugin).(i))[*,0])*100,color=clrs_plgns[i],pos=position,noerase=1,linestyle=2               
+                            
+                            ENDIF
+                        endfor
+                       
+                        
+                         ;PLotting of the interpolates corresponding to spectrum and filter points
+                        
+                         IF ct_spec NE 0 THEN BEGIN
+                         
+                             xx=((*!dustem_data.polfrac).wav)[idx_spec]
+                             yy=prediction[idx_spec]
+                             cgoplot,xx,yy*100,color='red',pos=position,psym=7,syms=2,noerase=1   
+                         ENDIF
+                         
+                         
+                         IF ct_filt NE 0 THEN BEGIN
+                             xx=((*!dustem_data.polfrac).wav)[idx_filt]
+                             yy=prediction[idx_filt]
+                             cgoplot,xx,yy*100,color='red',pos=position,psym=6,syms=2,noerase=1
+                         ENDIF
+                         
+                         
+                         ;Plotting of the total dust emission spectrum
+                         cgoplot,st.polsed.wav,specpolfrac*100,pos=position,noerase=1,/xlog,/ys,/xs
+                                        
+                                
+                    endif else begin ;The data points in the plot that remain unchanged.
+                        
+                        
+                        ;Plotting of frequency axis
+                        ;cgaxis, xaxis=1, xlog=1, xs=1, xminor=10, xticklen=0.05, xrange=((!const.c*1E6/(_extra.(ind_xr)))*1E-9),charsize=1.15,title=textoidl('\nu (GHz)')
+                        
+                        
+                        ;Spectrum points are treated before for plotting reasons. 
+                        if ct_spec ne 0 then begin 
+                            
+                            ;Plotting of spectrum data points (to be fitted)
+                            cgplot,((*!dustem_data.polfrac).wav)(idx_spec),((*!dustem_data.polfrac).values)(idx_spec)*100,/xlog,/ylog,/ys,xs=1,pos=position,noerase=1,xtit=' ',charsize=1.15,xtickformat='(A1)',color='Powder Blue',xr=xr,yr=[1.00E-03,50.00],psym=16,syms=0.8;,ytickformat='dstmwrp_exp'                
+                            
+                            ;Plotting of the spectrum error points
+                            rms=3.*((*!dustem_data.polfrac).sigma)(idx_spec)/2.
+                            cgerrplot,((*!dustem_data.polfrac).wav)(idx_spec),(((*!dustem_data.polfrac).values)(idx_spec)-rms)*100,(((*!dustem_data.polfrac).values)(idx_spec)+rms)*100,color='Powder Blue'
+
+                        endif
+                        
+                        if ct_filt ne 0 then begin
+                        
+                            ;Plotting of filter data points (to be fitted) + testing for prior plotting                         
+                            if ct_spec ne 0 then cgoplot,((*!dustem_data.polfrac).wav)(idx_filt),((*!dustem_data.polfrac).values)(idx_filt)*100,charsize=1.15,color='Dodger Blue',psym=16,syms=0.8,pos=position,/ys,xs=1,noerase=1,xtickformat='(A1)',xr=xr,yr=[1.00E-03,50.00],/xlog,/ylog else $;,ytickformat='dstmwrp_exp'
+                                cgplot,((*!dustem_data.polfrac).wav)(idx_filt),((*!dustem_data.polfrac).values)(idx_filt)*100,charsize=1.15,color='Dodger Blue',psym=16,syms=0.8,pos=position,/ys,xs=1,noerase=1,xtickformat='(A1)',xtit=' ',xr=xr,yr=[1.00E-03,50.00],/xlog,/ylog;,ytickformat='dstmwrp_exp'
+                            
+                            
+                            ;Plotting of the filter error points
+                            rms=3.*((*!dustem_data.polfrac).sigma)(idx_filt)/2.;/dustem_polfrac(idx_filt)
+                            cgerrplot,((*!dustem_data.polfrac).wav)(idx_filt),(((*!dustem_data.polfrac).values)(idx_filt)-rms)*100,(((*!dustem_data.polfrac).values)(idx_filt)+rms)*100,color='Dodger Blue'
+                            
+                        endif
+                        
+                        
+                        
+                        xyouts,pospltxt[0],pospltxt[1],textoidl('p_{\nu} (%)'),color=0,/normal,charsize=1.1
+                         
+                        ;Locating all the hidden data points (spectrum+filter)
+                        
+                        ;The filtering has been done prior (in the primary routine)
+                        match2,((*!dustem_data.polfrac).wav),((*!dustem_show.polfrac).wav),fit_polfracpts,show_polfracpts ;only show_polfracpts is needed
+                        idx_rmv_polfrac=where(show_polfracpts eq -1, ct_hdnpts) ; indices of the points to hide
+                        ;stop
+                        if ct_hdnpts ne 0 then begin ;Hidden data points are present 
+                            ;stop
+                        
+                            ;Locating the hidden spectrum and filter data points                                        
+                            idx_filt_hdn = where(((*!dustem_show.polfrac).filt_names)(idx_rmv_polfrac) NE 'SPECTRUM',ct_filt_hdn)
+                            idx_spec_hdn = where(((*!dustem_show.polfrac).filt_names)(idx_rmv_polfrac) EQ 'SPECTRUM',ct_spec_hdn)
+                            
+                            if ct_spec_hdn ne 0 then begin
+                                ;Plotting of hidden spectrum data points
+                                cgplot,(((*!dustem_show.polfrac).wav)[idx_rmv_polfrac])(idx_spec_hdn),(((*!dustem_show.polfrac).values)[idx_rmv_polfrac])(idx_spec_hdn)*100,pos=position,/xlog,/ylog,/ys,xs=1,noerase=1,charsize=1.15,xtickformat='(A1)',color='Black',xr=xr,yr=[1.00E-03,50.00],psym=8,syms=0.8
+                                
+                                ;Plotting of hidden spectrum error points
+                                rms=3.*(((*!dustem_show.polfrac).sigma)[idx_rmv_polfrac])(idx_spec_hdn)/2.
+                                cgerrplot,(((*!dustem_show.polfrac).wav)[idx_rmv_polfrac])(idx_spec_hdn),((((*!dustem_show.polfrac).values)[idx_rmv_polfrac])(idx_spec_hdn)-rms)*100,((((*!dustem_show.polfrac).values)[idx_rmv_polfrac])(idx_spec_hdn)+rms)*100,color='Black'
+                            endif                                           
+                            
+                            if ct_filt_hdn ne 0 then begin
+                                ;stop
+                                plotsym,0, /fill
+                                ;Plotting of hidden filter data points 
+                                cgplot,(((*!dustem_show.polfrac).wav)[idx_rmv_polfrac])(idx_filt_hdn),(((*!dustem_show.polfrac).values)[idx_rmv_polfrac])(idx_filt_hdn)*100,pos=position,/xlog,/ylog,/ys,xs=1,noerase=1,charsize=1.15,xtickformat='(A1)',color='Black',xr=xr,yr=[1.00E-03,50.00],psym=8,syms=0.8                                           
+                                ;stop
+                                ;Plotting of hidden filter error point
+                                rms=3.*(((*!dustem_show.polfrac).sigma)[idx_rmv_polfrac])(idx_filt_hdn)/2.
+                                cgerrplot,(((*!dustem_show.polfrac).wav)[idx_rmv_polfrac])(idx_filt_hdn),((((*!dustem_show.polfrac).values)[idx_rmv_polfrac])(idx_filt_hdn)-rms)*100,((((*!dustem_show.polfrac).values)[idx_rmv_polfrac])(idx_filt_hdn)+rms)*100,color='Black'
+                                 
+                            endif 
+                                                           
+                        endif 
+                        
+                    
+                    endelse 
+                
+            endelse 
  
         end
  
@@ -694,8 +858,8 @@ if keyword_set(dataset) then begin
                     endif else begin ;The data points in the plot that remain unchanged
                         ;stop
                         xtit=textoidl('\lambda (\mum)')
-                        if !run_pol then xtit = ''
-                        cgplot,wavs,wavs/wavs,/xlog,/ys,xs=1,pos=position,noerase=1,xtickformat='(A1)',color='Black',xr=xr,xtit=xtit,yr=[0.0,2.0],yticks=2,ymino=2,xticklen=0.1,ytickformat='(F6.2)',charsize=1.0
+                        ;if !run_pol then xtit = ''
+                        cgplot,wavs,wavs/wavs,/xlog,/ys,xs=1,pos=position,noerase=1,color='Black',xr=xr,xtit=xtit,yr=[0.0,2.0],yticks=2,ymino=2,xticklen=0.1,ytickformat='(F6.2)',charsize=1.0
                         xyouts,pospltxt[0],pospltxt[1],textoidl('norm'),color=0,/normal,charsize=1.1
  
                     endelse
@@ -1053,8 +1217,8 @@ if keyword_set(dataset) then begin
                     endif else begin ;The data points in the plot that remain unchanged
                         ;stop
                         xtit=textoidl('\lambda (\mum)')
-                        if !run_pol then xtit = ''
-                        cgplot,wavs,wavs/wavs,/xlog,/ys,xs=1,pos=position,noerase=1,xtickformat='(A1)',color='Black',xr=xr,xtit=xtit,yr=[0.0,2.0],yticks=2,ymino=2,xticklen=0.1,ytickformat='(F6.2)',charsize=1.0
+                        ;if !run_pol then xtit = ''
+                        cgplot,wavs,wavs/wavs,/xlog,/ys,xs=1,pos=position,noerase=1,color='Black',xr=xr,xtit=xtit,yr=[0.0,2.0],yticks=2,ymino=2,xticklen=0.1,ytickformat='(F6.2)',charsize=1.0
                         xyouts,pospltxt[0],pospltxt[1],textoidl('norm'),color=0,/normal,charsize=1.1
  
                     endelse
@@ -75,6 +75,7 @@ IF keyword_set(replace) THEN BEGIN
 	  _extra_out=_extra
 	  tags=tag_names(_extra)
 	  ind=where(tags EQ tag_name,count)
+	  ;stop
 	  IF count NE 0 THEN BEGIN
 	  	_extra_out.(ind[0])=tag_value
 	  ENDIF ELSE BEGIN