Merge branch 'master' of https://gitlab.irap.omp.eu/OV-GSO-DC/dustem-wrapper_idl

Jean-Philippe Bernard
2 parents bec7b3e8 92e2e66f
Showing 15 changed files with 272 additions and 192 deletions Show diff stats
src/idl/dustem_check_data.pro
src/idl/dustem_compute_polsed.pro
src/idl/dustem_init.pro
src/idl/dustem_init_parinfo.pro
src/idl/dustem_mpfit_run.pro
src/idl/dustem_plot_fit_sed.pro
src/idl/dustem_plot_polsed.pro
src/idl/dustem_plugin_freefree.pro
src/idl/dustem_plugin_mdfy_spinning_pol.pro
src/idl/dustem_plugin_modify_dust_polarization.pro
src/idl/dustem_plugin_stellar_population.pro
src/idl/dustem_read_all_web3p8.pro
src/idl/dustem_sed_plot.pro
src/idl/dustem_set_data.pro
src/idl_misc/JPBLib_for_Dustemwrap/General/convert_mk_mjy.pro
-FUNCTION dustem_check_data, sed=sed,ext=ext,ist=ist,pst=pst,qst=qst,ust=ust,exst=exst,pexst=pexst,qexst=qexst,uexst=uexst
+FUNCTION dustem_check_data, sed=sed,ext=ext,fpol=fpol,ist=ist,pst=pst,fpst=fpst,qst=qst,ust=ust,exst=exst,pexst=pexst,qexst=qexst,uexst=uexst
 ;THIS PROCEDURE CHECKS THE DATA PROVIDED BY THE USER AND RETURNS DIFFERENT STRUCTURES TO DEFINE (SET) THE COMPOSITE DATA FITTING STRUCTURE (!dustem_data)
@@ -23,6 +23,7 @@ IF KEYWORD_SET(sed) THEN BEGIN
     ;INITIALIZING OUTPUT STRUCTURES
     ist = obs_st
     pst = ist
+    fpst = ist
     qst = ist
     ust = ist
     ;This allows the user to fit polsed, qsed and used individually.
@@ -39,8 +40,7 @@ IF KEYWORD_SET(sed) THEN BEGIN
         read, ans, prompt='Would you like automatic format filtering? (Y/N):'
         IF strupcase(ans) EQ 'N' THEN BEGIN
             message,'You have chosen not to use automatic filtering. Please adjust the SED format.',/continue
-            message, 'Aborting ...',/continue
-            stop
+            message, 'Aborting ...'
         ENDIF
     ENDIF 
     IF strupcase(ans) EQ 'Y' THEN message, 'You have chosen automatic format filtering. Removal of undefined data points...',/continue $
@@ -59,11 +59,11 @@ IF KEYWORD_SET(sed) THEN BEGIN
     ans=''
     pst.values=sed.largeP
     pst.sigma=sqrt(abs(sed.sigma_largeP))
-    ;Testig if there is LargeP data. If not move to qsed data
+    ;Testig if there is LargeP data. If not move to polfrac data (fpol)
     ind = where(sed.largeP NE la_undef() and finite(sed.largeP) EQ 1 and sed.sigma_largeP NE la_undef() and finite(sed.sigma_largeP) EQ 1,NP) ;Testing on the nan values is needed to move to testing the other data fields (because if the dataset is full of Nans but not la_undef(), the code will think it's valid dta)
     IF NP EQ 0 THEN BEGIN
         pst = ptr_new()
-        goto, qzone
+        goto, fpzone
     ENDIF ELSE BEGIN
     ind = where(finite(sed.largeP) EQ 0 or finite(sed.sigma_largeP) EQ 0, countnokp) ;Format is not okay if there are nan values but the wrapper recognizes the la_undef() value 
     IF countnokp NE 0 THEN BEGIN 
@@ -72,8 +72,7 @@ IF KEYWORD_SET(sed) THEN BEGIN
         read, ans, prompt='Would you like automatic format filtering? (Y/N):' 
         IF strupcase(ans) EQ 'N' THEN BEGIN
             message,'You have chosen not to use automatic filtering. Please adjust the SED format.',/continue
-            message, 'Aborting ...',/continue
-            stop
+            message, 'Aborting ...'
         ENDIF
     ENDIF
     IF strupcase(ans) EQ 'Y' THEN message, 'You have chosen automatic format filtering. Removal of undefined data points...',/continue $
@@ -86,8 +85,41 @@ IF KEYWORD_SET(sed) THEN BEGIN
         values = pst.values(testp)
         sigma = pst.sigma(testp)
         pst={instru_names:instru_names,filt_names:filt_names,wav:wav,values:values,sigma:sigma}
-    ENDIF 
-    ENDELSE       
+    ENDIF
+    ENDELSE 
+    fpzone:
+    ans=''
+    fpst.values=sed.smallp
+    fpst.sigma=sqrt(abs(sed.sigma_smallp))
+    ;Testig if there is smallp data. If not move to qsed data
+    ind = where(sed.smallP NE la_undef() and finite(sed.smallP) EQ 1 and sed.sigma_smallP NE la_undef() and finite(sed.sigma_smallP) EQ 1,NSP) ;Testing on the nan values is needed to move to testing the other data fields (because if the dataset is full of Nans but not la_undef(), the code will think it's valid dta)
+    IF NSP EQ 0 THEN BEGIN
+        fpst = ptr_new()
+        ;goto, qzone
+        goto, qzone
+    ENDIF ELSE BEGIN
+    ind = where(finite(sed.smallP) EQ 0 or finite(sed.sigma_smallP) EQ 0, countnoksp) ;Format is not okay if there are nan values but the wrapper recognizes the la_undef() value 
+    IF countnoksp NE 0 THEN BEGIN 
+        message, "SED format isn't met. NaN value(s) detected in SmallP data.",/continue
+        message, 'Without automatic format fitering you cannot proceed.',/continue
+        read, ans, prompt='Would you like automatic format filtering? (Y/N):' 
+        IF strupcase(ans) EQ 'N' THEN BEGIN
+            message,'You have chosen not to use automatic filtering. Please adjust the SED format.',/continue
+            message, 'Aborting ...'
+        ENDIF
+    ENDIF
+    IF strupcase(ans) EQ 'Y' THEN message, 'You have chosen automatic format filtering. Removal of undefined data points...',/continue $
+    ELSE message, 'Data format is met. Removal of undefined data points...',/continue
+    testsp=where(sed.sigma_smallP NE 0. and finite(sed.sigma_smallP) EQ 1 and sed.sigma_smallP NE la_undef() and sed.smallP NE la_undef() and finite(sed.smallP) EQ 1,ctestsp)                
+    IF ctestsp NE 0 THEN BEGIN 
+        instru_names = fpst.instru_names(testsp)
+        filt_names = fpst.filt_names(testsp)
+        wav = fpst.wav(testsp)
+        values = fpst.values(testsp)
+        sigma = fpst.sigma(testsp)
+        fpst={instru_names:instru_names,filt_names:filt_names,wav:wav,values:values,sigma:sigma}
+    ENDIF
+    ENDELSE
     qzone:
     ans=''
     qst.values=sed.StokesQ
@@ -105,8 +137,7 @@ IF KEYWORD_SET(sed) THEN BEGIN
         read, ans, prompt='Would you like automatic format filtering? (Y/N):' 
         IF strupcase(ans) EQ 'N' THEN BEGIN
             message,'You have chosen not to use automatic filtering. Please adjust the SED format.',/continue
-            message, 'Aborting ...',/continue
-            stop
+            message, 'Aborting ...'
         ENDIF
     ENDIF
     IF strupcase(ans) EQ 'Y' THEN message, 'You have chosen automatic format filtering. Removal of undefined data points...',/continue $
@@ -138,8 +169,7 @@ IF KEYWORD_SET(sed) THEN BEGIN
         read, ans, prompt='Would you like automatic format filtering? (Y/N):' 
         IF strupcase(ans) EQ 'N' THEN BEGIN
             message,'You have chosen not to use automatic filtering. Please adjust the SED format.',/continue
-            message, 'Aborting ...',/continue
-            stop
+            message, 'Aborting ...'
         ENDIF
     ENDIF  
     IF strupcase(ans) EQ 'Y' THEN message, 'You have chosen automatic format filtering. Removal of undefined data points...',/continue $
@@ -191,8 +221,7 @@ IF KEYWORD_SET(ext) THEN begin
         read, ans, prompt='Would you like automatic format filtering? (Y/N):' 
         IF strupcase(ans) EQ 'N' THEN BEGIN
             message,'You have chosen not to use automatic filtering. Please adjust the SED format.',/continue
-            message, 'Aborting ...',/continue
-            stop
+            message, 'Aborting ...'
         ENDIF
     ENDIF      
     IF strupcase(ans) EQ 'Y' THEN message, 'You have chosen automatic format filtering. Removal of undefined data points...',/continue $ 
@@ -225,8 +254,7 @@ IF KEYWORD_SET(ext) THEN begin
         read, ans, prompt='Would you like automatic format filtering? (Y/N):' 
         IF strupcase(ans) EQ 'N' THEN BEGIN
             message,'You have chosen not to use automatic filtering. Please adjust the SED format.',/continue
-            message, 'Aborting ...',/continue
-            stop
+            message, 'Aborting ...'
         ENDIF
     ENDIF
     IF strupcase(ans) EQ 'Y' THEN message, 'You have chosen automatic format filtering. Removal of undefined data points...',/continue $ 
@@ -260,8 +288,7 @@ IF KEYWORD_SET(ext) THEN begin
         read, ans, prompt='Would you like automatic format filtering? (Y/N):' 
         IF strupcase(ans) EQ 'N' THEN BEGIN
             message,'You have chosen not to use automatic filtering. Please adjust the SED format.',/continue
-            message, 'Aborting ...',/continue
-            stop
+            message, 'Aborting ...'
         ENDIF
     ENDIF
     IF strupcase(ans) EQ 'Y' THEN message, 'You have chosen automatic format filtering. Removal of undefined data points...',/continue $
@@ -294,8 +321,7 @@ IF KEYWORD_SET(ext) THEN begin
         read, ans, prompt='Would you like automatic format filtering? (Y/N):' 
         IF strupcase(ans) EQ 'N' THEN BEGIN
             message,'You have chosen not to use automatic filtering. Please adjust the SED format.',/continue
-            message, 'Aborting ...',/continue
-            stop
+            message, 'Aborting ...'
         ENDIF
     ENDIF
     IF strupcase(ans) EQ 'Y' THEN message, 'You have chosen automatic format filtering. Removal of undefined data points...',/continue $
@@ -48,14 +48,14 @@ IF not keyword_set(stp) THEN BEGIN
 ENDIF
 If keyword_set(result) THEN result=stp
-fact = 1.e4*1.E20/(4.*!pi)/(3.e8/1.e-6/(stp.polsed).wav)*1.e20/1.e7
+fact = 1.e4*(*!dustem_HCD)/(4.*!pi)/(3.e8/1.e-6/(stp.polsed).wav)*1.e20/1.e7
 P = stp.polsed.em_tot*fact;1.00e17/(4.*!pi) ;This is Polarized intensity P
-I = stp.sed.em_tot*fact;1.00e17/(4.*!pi)    ;This is Total intensity I
+stI = stp.sed.em_tot*fact;1.00e17/(4.*!pi)    ;This is Total intensity I
 out=0.
 Nwaves=(size(P))[1]
-frac=P/I
+frac=P/stI
 tes=where(finite(frac) eq 0)
 frac(tes)=0.
@@ -77,7 +77,7 @@ FOR i=0L,n_tags(*!dustem_scope)-1 DO BEGIN
 ENDFOR    
-IF ~isa(Q_spec) && ~isa(U_spec) THEN polar_ippsi2iqu,I,Q_spec,U_spec,frac,replicate(0.,Nwaves) 
+IF ~isa(Q_spec) && ~isa(U_spec) THEN polar_ippsi2iqu,stI,Q_spec,U_spec,frac,replicate(0.,Nwaves) 
 FOR i=0L,n_tags(*!dustem_scope)-1 DO BEGIN
@@ -92,13 +92,12 @@ FOR i=0L,n_tags(*!dustem_scope)-1 DO BEGIN
 ENDFOR
-
 P=sqrt(Q_spec^2+U_spec^2)
 ENDIF
-dustem_polsed = (*!dustem_data.polsed).values * 0.
-
+IF isa(!dustem_data.polsed) THEN dustem_polsed = (*!dustem_data.polsed).values * 0. ELSE dustem_polsed = (*!dustem_data.polfrac).values * 0. 
+;stop
 if not isarray(P) THEN stop 
@@ -109,22 +108,25 @@ ENDIF ELSE BEGIN
   message,'SKIPPING color correction calculations',/info
 ENDELSE
-ind_polsed=where((*!dustem_data.polsed).filt_names NE 'SPECTRUM',count_polsed)
-  
-IF count_polsed NE 0 THEN BEGIN 
-  filter_names=((*!dustem_data.polsed).filt_names)(ind_polsed)
+IF isa(!dustem_data.polsed) THEN ind_polsed=where((*!dustem_data.polsed).filt_names NE 'SPECTRUM',count_polsed) ELSE ind_polsed=where((*!dustem_data.polfrac).filt_names NE 'SPECTRUM',count_polsed) 
+IF count_polsed NE 0 THEN BEGIN
+  IF isa(!dustem_data.polsed) THEN filter_names=((*!dustem_data.polsed).filt_names)(ind_polsed) ELSE filter_names=((*!dustem_data.polfrac).filt_names)(ind_polsed)
   spolsed=dustem_cc(stp.polsed.wav,P,filter_names,cc=cc)
   dustem_polsed[ind_polsed]=spolsed
-     
 ENDIF
-;For spectrum data points, interpolate in log-log
+
+;stop
+
+
+;For spectrum data points, interpolate in log-log.
 ;Linear interpolation leads to wrong values, in particular where few
 ;wavelengths points exist in the model (long wavelengths).
-ind_spec=where((*!dustem_data.polsed).filt_names EQ 'SPECTRUM',count_spec)
+
+IF isa(!dustem_data.polsed) THEN ind_spec=where((*!dustem_data.polsed).filt_names EQ 'SPECTRUM',count_spec) ELSE ind_spec=where((*!dustem_data.polfrac).filt_names EQ 'SPECTRUM',count_spec)
 IF count_spec NE 0 THEN BEGIN
-  dustem_polsed(ind_spec)=interpol(P,stp.polsed.wav,(((*!dustem_data.polsed).wav)(ind_spec)))
+  IF isa(!dustem_data.polsed) THEN dustem_polsed(ind_spec)=interpol(P,stp.polsed.wav,(((*!dustem_data.polsed).wav)(ind_spec))) ELSE dustem_polsed(ind_spec)=interpol(P,stp.polsed.wav,(((*!dustem_data.polfrac).wav)(ind_spec)))
 ENDIF
 out_st=stp
@@ -200,7 +200,11 @@ defsysv,&#39;!dustem_show_plot&#39;,1           ;0=show no plot
 ;== This is for the default mode
 dustem_inputs={GRAIN:'GRAIN.DAT',ISRF:'ISRF.DAT',ALIGN:'ALIGN.DAT'}
 !dustem_inputs=ptr_new(dustem_inputs)
+
+
 defsysv, '!dustem_HCD', ptr_new(1.00E20) ;This is a dumb thing I added but will shortly remove it 
+
+
 ;===This is the Desert option
 IF keyword_set(model) THEN BEGIN
@@ -247,7 +251,11 @@ defsysv, &#39;!dustem_model&#39;, model
     'NY_MODELA':BEGIN
       (*!dustem_inputs).grain='GRAIN_NY_MODELA.DAT'
       (*!dustem_inputs).align='ALIGN_G17_MODELA.DAT' 
-    END  
+    END
+    ;This is a user-defined model    
+    'USER_MODEL':BEGIN
+      (*!dustem_inputs).grain='GRAIN_USER.DAT' 
+    END
     ELSE:BEGIN
       message,'model '+model+' unknown',/continue
       message,'Known models are COMPIEGNE_ETAL2010,DBP90,DL01,DL07,AJ13,G17_MODELA,G17_MODELB,G17_MODELC,G17_MODELD,THEMIS',/continue
@@ -10,7 +10,7 @@ PRO dustem_init_parinfo,pd,iv, $
 dustem_set_func_ind,pd,iv
-one_parinfo={value:0.D0, fixed:0, limited:[0,0], limits:[0.D0,0], parname:'', relstep:0., step:0.5, mpside:0, tied:''}
+one_parinfo={value:0.D0, fixed:0, limited:[0,0], limits:[0.D0,0], parname:'', relstep:0.5, step:0.0, mpside:0, tied:''}
 Npar=n_elements(pd)
 dustem_parinfo=replicate(one_parinfo,Npar)
 FOR i=0L,Npar-1 DO BEGIN
@@ -60,5 +60,5 @@ ENDIF
 suite2:
 defsysv,'!dustem_parinfo',ptr_new(dustem_parinfo)
-
+;stop
 END
-FUNCTION dustem_mpfit_run,x,p_min,_EXTRA=extra
+FUNCTION dustem_mpfit_run,x,p_min,_EXTRA=_extra
 ;+
 ; NAME:
@@ -55,12 +55,6 @@ dustem_activate_plugins,p_dim/(*(*!dustem_fit).param_init_values)
 st = dustem_run(p_dim)
 st_model=dustem_read_all(!dustem_soft_dir)
-
- 
-;problem that needs fixing - case of the composite ISRF
-; st=st if you want to make the isrf thingy work. This is because dustem_activate_plugins requires the st structure but in the case of the composite ISRF the problem is that the wrapper's first read isrf is a MATHIS one so the first abundances the wrapper gets from dustem are the product of dust heated by the MATHIS ISRF
-
-
 chi2_sed      = 0  
 chi2_ext      = 0  
 chi2_polext   = 0 
@@ -105,7 +99,7 @@ dustem_out = [0]
 ind_data_tag=dustem_data_tag(count=count_data_tag)
 ;if exist('./noplot') then !dustem_show_plot = 0 else if !dustem_show_plot eq 0 then !dustem_show_plot = 1
-
+;stop
 FOR i_tag=0,count_data_tag-1 DO BEGIN
 ;stop
@@ -113,11 +107,9 @@ FOR i_tag=0,count_data_tag-1 DO BEGIN
 ;print,i_tag,ind_data_tag(i_tag)
 ;print,tagnames(ind_data_tag(i_tag))
-
     CASE tagnames(ind_data_tag(i_tag)) OF
-
+        
 		'SED'	: BEGIN
-            
 	        dustem_sed = dustem_compute_sed(p_dim,st);dustem_compute_sed(p_dim,st,cont=cont,freefree=freefree,synchrotron=synchrotron)
 	        chi2_sed   = total(((*!dustem_data.sed).values-dustem_sed)^2/(*!dustem_data.sed).sigma^2)
 			n_sed = n_elements((*!dustem_data.sed).values)
@@ -136,22 +128,96 @@ 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
+	            
 	            ;title='DUSTEM WRAP('+strtrim(string(win),2)+')'
 	            window,win,title='DUSTEM WRAP (SED)' 
 	            ;cgwindow,WMULTI=[2,2,1],/CURRENT
-                
-	            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
-	           
-	            
+                ;I.C
+                ;undeprecated the use of polfrac here
+                ;**********
+                ;I added this block so I can plot polfrac inside the sed window
+                ;I will probably remove it
+                ;**********
+                ;test if polfrac is in tagnames: 
+                ;another problem is that sed and pol data don't have the same filters. 
+                varvar=where(tagnames eq 'POLFRAC', testfpol)                
+                if testfpol ne 0 then 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.
+                        ;I'm not so sure about this...
+                        
+                        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
+                            
+                                ;I think I will have errors here!! The thing is I'm trying to make this work without (*!dustem_data.polsed) which is hard.
+                                
+                                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 
+                        ;stop
+        	        	;stop
+        	        	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
+                endif
+                ;stop
+	            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
+	       
 	            ;stop
 		 	ENDIF
 	        dustem_out = [ dustem_out, dustem_sed ]
-
+            
 		END
 		'EXT' 	: BEGIN 
@@ -394,7 +460,7 @@ FOR i_tag=0,count_data_tag-1 DO BEGIN
 		        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')
+                ytitstq=textoidl('Q_\nu (W/m^2/Hz/sr) for N_H=10^{20} H/cm^2')
                 xtit=textoidl('\lambda (\mum)')
                 xr=[10,2e4]
                 normpol=Q_spec*fact
@@ -579,9 +645,7 @@ FOR i_tag=0,count_data_tag-1 DO BEGIN
   		END
- 		
-	    ELSE : stop
-
+	    ELSE : IF tagnames(ind_data_tag(i_tag)) NE 'POLFRAC' THEN stop 
     ENDCASE
 ENDFOR
@@ -2,7 +2,7 @@
                          res=res,errors=errors,chi2=chi2,rchi2=rchi2, $
                          no_spec_error=no_spec_error,title=title,$;xtit=xtit,ytit=ytit,xr=xr,yr=yr $
                          help=help,win=win,ps=ps, $
-                         legend_xpos=legend_xpos,legend_ypos=legend_ypos,legend_offset=legend_offset, $
+                         legend_xpos=legend_xpos,legend_ypos=legend_ypos,legend_offset=legend_offset,fpol=fpol, $
                          _extra=_extra
 ; NAME:
@@ -59,11 +59,10 @@ IF keyword_set(ps) THEN BEGIN
 ENDIF ELSE BEGIN
   set_plot,'X'
-    IF keyword_set(win) then window,win
+    IF keyword_set(win) then window,win;,xsize=600,ysize=800
 ENDELSE
 fact = 1.e4*(*!dustem_HCD)/(4.*!pi)/(3.e8/1.e-6/st.sed.wav)*1.e20/1.e7
-
 ;use_col_data_filt=70
 ;use_col_sed_spec=170
 use_col_data_filt='blue'
@@ -89,7 +88,7 @@ ENDELSE
 spec = st.sed.em_tot * fact
-
+if keyword_set(fpol) then specpol = st.polsed.em_tot * fact
 ;ADDING PLUGIN(S) TO SPECTRUM----------------
 ;if n_tags(!dustem_data.sed) gt 1 then begin 
 scopes=tag_names((*!dustem_scope))
@@ -97,6 +96,14 @@ IF scopes[0] NE &#39;NONE&#39; THEN BEGIN
 ;IF ptr_valid(!dustem_plugin) THEN BEGIN
   FOR i=0L,n_tags(*!dustem_scope)-1 DO BEGIN
     IF total(strsplit((*(*!dustem_scope).(i)),'+',/extract) EQ 'ADD_SED') THEN spec+=(*(*!dustem_plugin).(i))[*,0]
+    if keyword_set(fpol) then begin
+        IF total(strsplit((*(*!dustem_scope).(i)),'+',/extract) EQ 'REPLACE_QSED') THEN specpol=sqrt(((*(*!dustem_plugin).(i))[*,1])^2+((*(*!dustem_plugin).(i))[*,2])^2)
+    endif
+  ENDFOR
+  FOR i=0L,n_tags(*!dustem_scope)-1 DO BEGIN
+    if keyword_set(fpol) then begin    
+        IF total(strsplit((*(*!dustem_scope).(i)),'+',/extract) EQ 'ADD_QSED') THEN specpol+=sqrt(((*(*!dustem_plugin).(i))[*,1])^2+((*(*!dustem_plugin).(i))[*,2])^2)
+    endif  
   ENDFOR
 ENDIF
 ;endif
@@ -127,12 +134,23 @@ IF keyword_set(xtit) THEN xtit = xtit ELSE xtit = textoidl(&#39;\lambda (\mum)&#39;)
 IF keyword_set(ytit) THEN ytit = ytit ELSE ytit = textoidl('Brightness/N_H (MJy/sr/H)')
+
 ;deffo define a title variable here. The procedures are not communicating with each other. 
+;###############################
+if keyword_set(fpol) then begin
+    if !run_lin then begin
+    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
+;###############################
+
+
 ;cgplot,(*!dustem_data.sed).wav,(*!dustem_data.sed).values/norm,/nodata,xtit='',ytit=textoidl('Brightness/N_H (MJy/sr/H)'),tit=title,ylog=1,/xlog,xr=xr,yr=yr,/ys,/xs,position=[0.12,0.35,0.96,0.90],xtickformat='(A1)'
-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
+;stop
 ;cgplot,st.polsed.wav,Q_sed*fact,xtit='',ytit=ytitstq,tit=titstq,ylog=ylog,/xlog,xr=xr,yr=yr,/ys,/xs,position=[0.17,0.35,0.95,0.95],xtickformat='(A1)'
 ;this is where you will add the normalized sed
@@ -185,8 +203,8 @@ IF !dustem_show_plot EQ 2 THEN BEGIN
 ENDIF ELSE BEGIN
 	norm = spec * 0. + 1
 ENDELSE
+use_cols[1]='Cornflower'
 FOR i=0L,Ngrains-1 DO BEGIN
-  use_cols[1]='Cornflower'
   cgoplot,st.sed.wav,st.sed.(i+1)*fact/norm,color=use_cols[i],linestyle=use_lines[i]
 ENDFOR
@@ -213,6 +231,7 @@ ENDIF
 cgoplot,st.sed.wav,spec/norm,color=use_col_tot,linestyle=use_line_tot
+if keyword_set(fpol) then cgoplot,st.polsed.wav,specpol,color='Teal',linestyle=4
 ;plot the normalized data as well. 
 ;----------------------------------------
@@ -221,7 +240,7 @@ frmt0=&#39;(A36)&#39;
 frmt1='(1E10.2)'
 frmt2='(F7.2)'
 use_legend_xpos=0.50
-use_legend_ypos=0.84
+if keyword_set(fpol) then use_legend_ypos=0.70 ELSE use_legend_ypos=0.84
 use_legend_offset=0.03             ;This is the offset between lines of the legend in normalized units
 legend_charsize=0.86
@@ -236,7 +255,7 @@ IF keyword_set(legend_ypos) THEN use_legend_ypos=legend_ypos
 IF keyword_set(legend_offset) THEN use_legend_offset=legend_offset
 IF !d.name NE 'PS' THEN cleanplot
-;oo=0
+
 tg=0
 prv_str=''
 IF keyword_set(res) THEN BEGIN
@@ -327,25 +346,31 @@ IF keyword_set(res) THEN BEGIN
 ENDIF
 IF keyword_set(chi2) THEN BEGIN
   xxpos=1.29*use_legend_xpos
-  yypos=1.13*use_legend_ypos
+  yypos=1.13*0.84;use_legend_ypos
   xyouts,xxpos,yypos,string('chi2=',format=frmt0)+string(chi2,format=frmt2),color=0,/normal,charsize=legend_charsize
 ENDIF
 IF keyword_set(rchi2) THEN BEGIN
   xxpos=1.29*use_legend_xpos
-  yypos=1.13*use_legend_ypos-use_legend_offset
+  yypos=1.13*0.84-use_legend_offset;use_legend_ypos-use_legend_offset
   xyouts,xxpos,yypos,string('red. chi2=',format=frmt0)+string(rchi2,format=frmt2),color=0,/normal,charsize=legend_charsize
 ENDIF
 xtit=textoidl('\lambda (\mum)')
 ;cgplot,(*!dustem_data.sed).wav,(*!dustem_data.sed).values/sed,/nodata,xtit=xtit,ytit='Normalized',tit='',/xlog,xr=xr,/ys,/xs,yr=[0,2],ylog=0,position=[0.12,0.14,0.96,0.35],/noerase,yticks=2,ymino=2,xticklen=0.1
 ;stop
+
 IF keyword_set(_extra) THEN BEGIN
   extra_kept={XR:[0.,0.]}
   extra_tags=tag_names(_extra)
   ind=where(extra_tags EQ 'XR',count)
   IF count NE 0 THEN extra_kept.XR=_extra.(ind[0]) ;ELSE extra_kept=0
 ENDIF
-cgplot,(*!dustem_data.sed).wav,(*!dustem_data.sed).values/sed,_extra=extra_kept,/nodata,xtit=xtit,ytit='Normalized',tit='',/xlog,/ys,/xs,yr=[0,2],ylog=0,position=[0.12,0.14,0.96,0.35],/noerase,yticks=2,ymino=2,xticklen=0.1,ytickformat='(F6.2)',charsize=1.3
+
+;stop
+;,position=[0.12,0.14,0.96,0.35] old position before your polfrac display modification
+if keyword_set(fpol) then begin
+cgplot,(*!dustem_data.sed).wav,(*!dustem_data.sed).values/sed,_extra=extra_kept,/nodata,xtit=xtit,ytit='Normalized',tit='',/xlog,/ys,/xs,yr=[0,2],ylog=0,position=[0.12,0.1,0.96,0.25],/noerase,yticks=2,ymino=2,xticklen=0.1,ytickformat='(F6.2)',charsize=1.3
+endif else cgplot,(*!dustem_data.sed).wav,(*!dustem_data.sed).values/sed,_extra=extra_kept,/nodata,xtit=xtit,ytit='Normalized',tit='',/xlog,/ys,/xs,yr=[0,2],ylog=0,position=[0.12,0.14,0.96,0.35],/noerase,yticks=2,ymino=2,xticklen=0.1,ytickformat='(F6.2)',charsize=1.3
 ;plot the normalized data as well. 
 IF count_spec NE 0 THEN BEGIN
   xx=((*!dustem_data.sed).wav)[ind_spec]
@@ -361,6 +386,7 @@ IF count_filt NE 0 THEN BEGIN
   trois_sigma=(3.*((*!dustem_data.sed).sigma)/2./sed)[ind_filt]
   cgerrplot,xx,yy-trois_sigma,yy+trois_sigma,color='Dodger Blue'
 ENDIF
+
 ;stop
 ;cgoplot,(*!dustem_data.sed).wav,(*!dustem_data.sed).values/sed,psym=16,symsize=1,thick=2,color='Dodger Blue'
 cgoplot,10^!x.crange,replicate(1.,2),color='black',linestyle=0
@@ -368,6 +394,24 @@ cgoplot,10^!x.crange,replicate(1.,2),color=&#39;black&#39;,linestyle=0
 ;cgerrplot,((*!dustem_data.sed).wav),((*!dustem_data.sed).values)/sed-3.*((*!dustem_data.sed).sigma)/2./sed,((*!dustem_data.sed).values)/sed+3.*((*!dustem_data.sed).sigma)/2./sed,color='Dodger Blue'
+if keyword_set(fpol) then begin
+    plotsym,0,/fill
+    cgplot, (*!dustem_data.polfrac).wav,(*!dustem_data.polfrac).values*100,/nodata,xtit='',tit=title,ytit='',/xlog,/ylog,/ys,/xs,position=[0.12,0.76,0.96,0.91],charsize=1.3,yr=[1.0E-1,30.0],/noerase,xticklen=0.1,xtickformat='(A1)',_extra=extra_kept;,xr=_extra.xr
+    xxpos=use_legend_xpos*0.35
+    yypos=use_legend_ypos*1.25
+    xyouts,xxpos,yypos,'Polarization fraction (%)',color=cgcolor('purple'),charsize=legend_charsize,/normal
+    yypos=use_legend_ypos*1.20
+    stringg='Polarization SED '+'('+textoidl('P_{\nu}')+')'
+    xyouts,xxpos,yypos,stringg,color=cgcolor('teal'),charsize=legend_charsize,/normal
+    pilotfx=fpol(0);((*!dustem_data.polfrac).values)(0);fpol(0);
+    cgoplot, (*!dustem_data.polfrac).wav,(*!dustem_data.polfrac).values*100,charsize=1.3,psym=8,syms=1,thick=2,color='Dodger Blue'   
+    cgoplot, st.polsed.wav,(specpol/spec)*100/pilotfx*((*!dustem_data.polfrac).values)(0),color='purple'
+    cgoplot, (*!dustem_data.polfrac).wav,fpol*100/pilotfx*((*!dustem_data.polfrac).values)(0),color='Red',psym=6,syms=1
+    ;fpol=abs(fpol)
+    ;delp=(*!dustem_data.polfrac).values*100-fpol*100/pilotfx*((*!dustem_data.polfrac).values)(0)
+    ;cgoplot, ((*!dustem_data.polfrac).wav),((*!dustem_data.polfrac).values*100+delp),charsize=1.3,color='black',linestyle=2
+    ;stop
+endif
 IF keyword_set(ps) THEN BEGIN
   device,/close
-PRO dustem_plot_polsed, st, p_dim, dustem_polsed, aligned=aligned, win=win
+        PRO dustem_plot_polsed, st, p_dim, dustem_polsed, aligned=aligned, win=win
 IF keyword_set(ps) THEN BEGIN
@@ -29,7 +29,7 @@ if not keyword_set(st) then begin
 ENDIF
-fact=1/(4.*!pi)/(3.e8/1.e-6/(st.polsed).wav)*1.00e17
+fact=1/(4.*!pi)*(*!dustem_HCD)*1.E-20/(3.e8/1.e-6/(st.polsed).wav)*1.00e17 ;applied correction here to take into account the HCD
 fact1=1.E-20 ;This will convert data (real data) that is in MJy/sr to W
 P = st.polsed.em_tot*fact    ;This is Polarized intensity P
@@ -42,14 +42,20 @@ IF keyword_set(help) THEN BEGIN
 ENDIF
 ;default values of input parameters
-Tgas=10000.                 ;default gas temperature
+Tgas=10000.           ;default gas temperature
 Amplitude=1.          ;Amplitude
-paramtag=['T_gas (K)','Amp']
+smallp=0.0            ;default polarization fraction
+psi=0.                ;default polarization angle
+paramtag=['T_gas (K)','Amp','p','Psi (deg)']
 IF keyword_set(key) THEN BEGIN 
   a=where(key EQ 1,count1)
   b=where(key EQ 2,count2)
+  c=where(key EQ 3,count3) ;default polarization fraction -newly added
+  d=where(key EQ 4,count4) ;default polarization angle -newly added
   IF count1 NE 0 then Tgas=(val(a))(0)
   IF count2 NE 0 then Amplitude=(val(b))(0)
+  IF count3 NE 0 then smallp=val[c[0]]   ;-newly added
+  IF count4 NE 0 then psi=val[d[0]]    ;-newly added
 ENDIF
 ;IF !dustem_which EQ 'DESERT' THEN BEGIN
@@ -59,12 +65,15 @@ ENDIF
 ;ENDELSE
 ;replaced by the following to ease the life of plugins writters
 lambir=dustem_get_wavelengths()
-
+Nwavs=n_elements(lambir)
 cmic=3.e14
 nu=cmic/lambir   ;Hz
 mjy=1           ;output is in MJy/sr
 lambir_ref=10000.
+
+output=fltarr(Nwavs,3) ;newly added
+
 ;stop
 ;use_method='Deschenes2008'
 ;use_method='WallsGabaud1998'
@@ -73,22 +82,22 @@ use_method=&#39;Dickinson2003_norm&#39;
 CASE use_method OF
     'WallsGabaud1998':BEGIN
 		em=1.   ;This is a stupid value for the Emission measure. Result will be rescaled based on I_halpha_R anyway
-		output=intensity_free_free(nu,Tgas,em,I_halpha_R=I_halpha_R,nHe=nHe,ergs=ergs,mjy=mjy)
+		output[*,0]=intensity_free_free(nu,Tgas,em,I_halpha_R=I_halpha_R,nHe=nHe,ergs=ergs,mjy=mjy)
 		fact=Amplitude/I_halpha_R
         ;=== normalize to the requested Halpha value.
-		output=output*fact
+		output[*,0]=output[*,0]*fact
 	END
     'WallsGabaud1998_scaled':BEGIN
 		em=1.   ;This is a stupid value for the Emission measure. Result will be rescaled based on I_halpha_R anyway
-		output=intensity_free_free(nu,Tgas,em,I_halpha_R=I_halpha_R,nHe=nHe,ergs=ergs,mjy=mjy)
-		norm=interpol(output,lambir,lambir_ref)
-		output=output/norm*Amplitude
+		output[*,0]=intensity_free_free(nu,Tgas,em,I_halpha_R=I_halpha_R,nHe=nHe,ergs=ergs,mjy=mjy)
+		norm=interpol(output[*,0],lambir,lambir_ref)
+		output[*,0]=output[*,0]/norm*Amplitude
 	END
 	'Deschenes2008':BEGIN
 		beta_freefree=2.+1./(10.48+1.5*alog(Tgas/8.e3)-alog(nu/1.e9))
-		output=nu^(-1.*beta_freefree)
-		norm=interpol(output,lambir,lambir_ref)
-		output=output/norm*Amplitude
+		output[*,0]=nu^(-1.*beta_freefree)
+		norm=interpol(output[*,0],lambir,lambir_ref)
+		output[*,0]=output[*,0]/norm*Amplitude
 	END
 	'Dickinson2003_norm':BEGIN
 		;stop
@@ -99,12 +108,17 @@ CASE use_method OF
 		Tb=8.369e3*a*(nu_ghz)^(-2.)*T4^0.667*10^(0.029*T4)*(1.+0.08)  ;in mK for 1 Rayleigh
 		convert_mk_mjy, lambir, Tb, I_Mjy, /RJ
 		norm=interpol(I_Mjy,lambir,lambir_ref)
-		output=I_Mjy/norm*Amplitude
+		output[*,0]=I_Mjy/norm*Amplitude
 	END
 ENDCASE
-scope='ADD_SED'
+scope='ADD_SED+ADD_QSED+ADD_USED'
+
+;polarization this actually need to 
+polar_ippsi2iqu,output[*,0],Q,U,replicate(smallp,Nwavs),replicate(psi,Nwavs)
+output[*,1]=Q
+output[*,2]=U
 the_end:
-FUNCTION dustem_plugin_mdfy_spinning_pol, key=key, val=val, scope=scope, help=help
+FUNCTION dustem_plugin_mdfy_spinning_pol, key=key, val=val, scope=scope, paramtags=paramtags, help=help
 ;+
 ; NAME:
@@ -65,6 +65,10 @@ IF keyword_set(key) THEN BEGIN
 ENDIF
+
+
+paramtags=['p','Psi (deg)']
+
 lambir=dustem_get_wavelengths()
 Nwaves=n_elements(lambir)
@@ -1,93 +0,0 @@
-Function dustem_plugin_modify_dust_polarization, key=key, val=val, scope=scope,help=help
-
-;+
-; NAME:
-;    dustem_create_stokes
-; PURPOSE:
-;    Produces stokes emission parameters for the dust and synchrotron components
-; CATEGORY:
-;    DUSTEM Wrapper
-; CALLING SEQUENCE:
-;    dustem_create_stokes(st,key=key,val=val)
-; INPUTS:
-;    st (st = dustem_run(p_dim))
-; OPTIONAL INPUT PARAMETERS:
-;    key  = input parameter number
-;    val  = input parameter value
-; OUTPUTS:
-;    out = array containing the stokes emission parameters associated to the dust/synchrotron component or a concatenated version for both
-; OPTIONAL OUTPUT PARAMETERS:
-;    None
-; ACCEPTED KEY-WORDS:
-;    help                  = if set, print this help
-; COMMON BLOCKS:
-;    None
-; SIDE EFFECTS:
-;    None
-; RESTRICTIONS:
-;    The dustem fortran code must be installed
-;    The dustem idl wrapper must be installed
-; PROCEDURE:
-;    This is a dustem plugin
-;-
-
-IF keyword_set(help) THEN BEGIN
-  doc_library,'dustem_create_stokes'
-  out=0.
-  goto,the_end
-ENDIF
-
-
-
-IF not KEYWORD_SET(st) THEN BEGIN
-IF  ptr_valid((*!dustem_fit).CURRENT_PARAM_VALUES) THEN BEGIN
-    p_di=(*(*!dustem_fit).CURRENT_PARAM_VALUES)
-    
-ENDIF ELSE BEGIN
-          p_di=(*(*!dustem_fit).PARAM_INIT_VALUES)    
-      ENDELSE 
-          
-st=dustem_run(p_di)      
-ENDIF
-
-
-psi=0.
-smallp=1.
-IF keyword_set(key) THEN BEGIN 
-  ind1=where(key EQ 1,count1)
-  ind2=where(key EQ 2,count2)
-  IF count1 NE 0 THEN smallp=val[ind1[0]] ; setting psi from pd
-  IF count2 NE 0 THEN psi=val[ind2[0]] ; setting psi from pd   
-ENDIF
-
-
-fact = 1.e4*1.E20/(4.*!pi)/(3.e8/1.e-6/((st.polsed).wav))*1.e20/1.e7
-P=((st.polsed).em_tot)*fact
-I=((st.sed).em_tot)*fact
-
-
-
-Nwaves=(size(I))[1]
-
-frac=P/I*smallp
-tes=where(finite(frac) eq 0)
-frac(tes)=0.
-
-
-polar_ippsi2iqu,I,Q,U,frac,replicate(psi,Nwaves) 
-
-out=fltarr(Nwaves,3)
-
-out[*,0]=I
-out[*,1]=Q
-out[*,2]=U
-
-scope='REPLACE_QSED+REPLACE_USED'
-
-return, out
-
-the_end:
-
-
-end
-
@@ -3288,8 +3288,6 @@ ENDELSE
 if keyword_set(key) then begin
 print, 'comp_pop is:'
 print, comp_pop
-endif
-
 ;integrating the mathis isrf
 ma_isrf_dat=!dustem_soft_dir+'data/ISRF_MATHIS.DAT'
 ma_isrf=dustem_read_isrf(ma_isrf_dat)
@@ -3303,6 +3301,9 @@ int_isrf=INT_TABULATED((ma_isrf.lambisrf)[un],(st.isrf)[un])
 print, 'G0_stellar_population='
 print, int_isrf/int_mathis
 ;stop
+
+endif
+
 the_end:
@@ -70,7 +70,7 @@ FOR i=0L,Ngrains-1 DO BEGIN
   st=dustem_read_qabs_lv(Qabs_file,silent=silent,read_densities=read_densities)
   st_qabs(i)=ptr_new(st)
 ENDFOR
-
+;stop
 ;=== Read heat capacities
 st_calor=ptrarr(Ngrains)
 FOR i=0L,Ngrains-1 DO BEGIN
-PRO dustem_sed_plot,p_min,_extra=_extra,function_name=function_name,pol=pol,legend_xpos=legend_xpos,legend_ypos=legend_ypos,ps=ps,png=png
+PRO dustem_sed_plot,p_min,_extra=_extra,function_name=function_name,pol=pol,legend_xpos=legend_xpos,legend_ypos=legend_ypos,ps=ps,png=png,fpol=fpol
 IF not keyword_set(function_name) THEN BEGIN
   ;Run dustem with as an interface to mpfitfun; what does this mean? This whole procedure will have to be changed
@@ -15,20 +15,20 @@ IF not keyword_set(function_name) THEN BEGIN
   st=out_st
   ;dustem_plot_fit_sed,st,dustem_sed,cont,freefree,synchrotron,_extra=_extra,legend_xpos=legend_xpos,legend_ypos=legend_ypos
   win=1;win-act_wins+1
-  IF keyword_set(ps) THEN dustem_plot_fit_sed,st,dustem_sed,win=win,_extra=_extra,legend_xpos=legend_xpos,legend_ypos=legend_ypos,ps=ps,png=png,use_model=use_model ELSE dustem_plot_fit_sed,st,dustem_sed,_extra=_extra,win=win,legend_xpos=legend_xpos,legend_ypos=legend_ypos,use_model=use_model
+  IF keyword_set(ps) THEN dustem_plot_fit_sed,st,dustem_sed,win=win,_extra=_extra,legend_xpos=legend_xpos,legend_ypos=legend_ypos,ps=ps,png=png,use_model=use_model,fpol=fpol ELSE dustem_plot_fit_sed,st,dustem_sed,_extra=_extra,win=win,legend_xpos=legend_xpos,legend_ypos=legend_ypos,use_model=use_model,fpol=fpol
   ;I do not understand the reason begind the stop below - either way polarization is handled by other procedures
   IF keyword_set(pol) THEN BEGIN
     ;stop
     IF tag_exist(!dustem_data,'polsed') THEN BEGIN
-        IF isa(*!dustem_data.polsed) THEN BEGIN
+        IF isa(!dustem_data.polsed) THEN BEGIN
             win+=1
             dustem_polsed=dustem_compute_polsed(p_min,st=st)
             dustem_plot_polsed, st, p_dim, dustem_polsed, aligned=aligned, win=win,_Extra=extra
         ENDIF
     ENDIF
     IF tag_exist(!dustem_data,'polext') THEN BEGIN
-        IF isa(*!dustem_data.polext) THEN BEGIN
+        IF isa(!dustem_data.polext) THEN BEGIN
             win+=1
             dustem_polext=dustem_compute_polext(p_min,st=st)
             dustem_plot_polext, st, p_dim, dustem_polext, aligned=aligned, win=win,_Extra=extra
-FUNCTION dustem_set_data,help=help,sed=sed,ext=ext,polfrac=polfrac,rchi2_weight=rchi2_weight,f_HI=f_HI
+FUNCTION dustem_set_data,help=help,sed=sed,ext=ext,fpol=fpol,rchi2_weight=rchi2_weight,f_HI=f_HI
 ;+
 ; NAME:
@@ -46,7 +46,6 @@ FUNCTION dustem_set_data,help=help,sed=sed,ext=ext,polfrac=polfrac,rchi2_weight=
 ;    The function uses two main input structures: SED and EXT. Polarization data is included within each. 
 ;-
-;fully deprecate the use of polfrac?===>YES should do so
 IF keyword_set(help) THEN BEGIN
   doc_library,'dustem_set_data'
@@ -61,8 +60,10 @@ ENDIF
 ;SETTING THE DATA STRUCTURES FOR THE FITTING PROCESS - SED
 IF keyword_set(sed) THEN BEGIN    
-               
-    xx=dustem_check_data(sed=sed,ist=ist,pst=pst,qst=qst,ust=ust)            
+    
+    
+    ;maybe to activate fpol here you need to set it to one? CHECK THIS OUT.           
+    if keyword_set(fpol) then xx=dustem_check_data(sed=sed,ist=ist,pst=pst,fpol=fpol,fpst=fpst,qst=qst,ust=ust) else xx=dustem_check_data(sed=sed,ist=ist,pst=pst,fpst=fpst,qst=qst,ust=ust)          
     ;stop
     ;=== fill !dustem_data
     IF isa(ist) THEN !dustem_data.sed = ptr_new(ist) ;ELSE !dustem_data.sed = ptr_new()   
@@ -79,6 +80,13 @@ IF keyword_set(sed) THEN BEGIN
             if keyword_set(f_HI) then if f_HI gt 0 then $
           	(*!dustem_data.polsed).values *= f_HI
         ENDIF
+        If isa(fpst) then BEGIN        
+            ;=== fill !dustem_data
+            !dustem_data.polfrac = ptr_new(fpst)
+            ; If f_HI is specified, multiply the data by f_HI
+            if keyword_set(f_HI) then if f_HI gt 0 then $
+          	(*!dustem_data.polfrac).values *= f_HI
+        ENDIF
         IF isa(qst) THEN BEGIN        
             ;=== fill !dustem_data
             !dustem_data.qsed = ptr_new(qst)
@@ -139,6 +147,7 @@ ENDIF
 If !run_pol THEN BEGIN
     ;SED
     IF isa(pst) THEN !fit_rchi2_weight.polsed=1.
+    IF isa(fpst) THEN !fit_rchi2_weight.polfrac=1.
     IF isa(qst) THEN !fit_rchi2_weight.qsed=1.
     IF isa(ust) THEN !fit_rchi2_weight.used=1.
@@ -155,6 +164,7 @@ IF keyword_set(rchi2_weight) THEN BEGIN
 	IF !run_pol THEN BEGIN
     	IF isa(pst) THEN !fit_rchi2_weight.polsed  = rchi2_weight.polsed
+    	IF isa(fpst) THEN !fit_rchi2_weight.polfrac = rchi2_weight.polfrac
     	IF isa(qst) THEN !fit_rchi2_weight.qsed  = rchi2_weight.qsed
     	IF isa(ust) THEN !fit_rchi2_weight.used  = rchi2_weight.used
@@ -16,7 +16,7 @@
 ; Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301  USA
 ;+
-;  Convert mK, either themodynamic or Rayleigh-Jeans in MJy/sr
+;  Convert mK, either thermodynamic or Rayleigh-Jeans in MJy/sr
 ; 
 ;
 ;  @param lam {in}{required}{type=float} wavelength in micron