FUNCTION dustem_compute_polext ,p_dim,sti,POLEXT_spec,SPEXT_spec,dustem_fpolext,dustem_ext=dustem_ext,out_st=out_st,_extra=extra

;CREATED AS A CONSEQUENCE OF THE INCLUSION OF THE PLUGINS


;written to homogonize DustEM

IF not keyword_set(sti) THEN BEGIN
  dustem_activate_plugins,p_dim/(*(*!dustem_fit).param_init_values),sti
  
ENDIF


POLEXT_spec = sti.polext.ext_tot * (*!dustem_HCD)/1.0e21  ; 

;Hard-coded test to set to zero potential negative values in the polarized extinction arrays:
;so far this has been seen in absorption arrays  
;This test needs to include all grain species that polarize

IF !run_pol THEN BEGIN
    
    ;THIS FALSE CORRECT  IT USING POLEXT_spec
    ;POLEXT_spec.abs_grain[2,*]
    ind_ngtv_plxt = where (POLEXT_spec LT 0, ct_ngtv_plxt)
   
    IF ct_ngtv_plxt NE 0 THEN (POLEXT_spec)[ind_ngtv_plxt] = 0.
        
        
ENDIF 

EXT_spec = sti.ext.ext_tot * (*!dustem_HCD)/1.0e21   ;This is Total intensity I

out=0.
Nwaves=(size(POLEXT_spec))[1]

frac=POLEXT_spec/EXT_spec
tes=where(finite(frac) eq 0)
frac(tes)=0.


scopes=tag_names((*!dustem_plugin))
IF scopes[0] NE 'NONE' THEN BEGIN
;IF ptr_valid(!dustem_plugin) THEN BEGIN
  for i=0L,n_tags(*!dustem_plugin)-1 do begin
    if total(strsplit((*(*!dustem_plugin).(i).scope),'+',/extract) eq 'ADD_EXT') then EXT_spec+=(*(*!dustem_plugin).(i).spec)[*,0]
  endfor
ENDIF


IF scopes[0] NE 'NONE' THEN BEGIN

    FOR i=0L,n_tags(*!dustem_plugin)-1 DO BEGIN   
        
        IF total(strsplit((*(*!dustem_plugin).(i).scope),'+',/extract) EQ 'REPLACE_POLEXT') THEN BEGIN
            
            QEXT_spec=(*(*!dustem_plugin).(i).spec)[*,1]
            UEXT_spec=(*(*!dustem_plugin).(i).spec)[*,2]
            
        ENDIF 
        
    ENDFOR   
ENDIF 
    
    IF ~isa(QEXT_spec) && ~isa(UEXT_spec) THEN BEGIN
        polar_ippsi2iqu,EXT_spec,QEXT_spec,UEXT_spec,frac,replicate(!dustem_psi_ext,Nwaves)
        SPEXT_spec = frac 
    ENDIF

IF scopes[0] NE 'NONE' THEN BEGIN

    FOR i=0L,n_tags(*!dustem_plugin)-1 DO BEGIN
        
        IF total(strsplit((*(*!dustem_plugin).(i).scope),'+',/extract) EQ 'ADD_POLEXT') THEN BEGIN
        
            QEXT_spec+=(*(*!dustem_plugin).(i).spec)[*,1]
            UEXT_spec+=(*(*!dustem_plugin).(i).spec)[*,2]
            
        ENDIF
      
    ENDFOR

ENDIF    

POLEXT_spec=sqrt(QEXT_spec^2+UEXT_spec^2)
SPEXT_spec = POLEXT_SPEC/EXT_spec


dustem_polext = (*!dustem_data.qext).values * 0.

    
if not isarray(POLEXT_spec) THEN stop ;I don't understand this test. The only thing it can indicate is a problem with the dust parameters or the fortran executable. But in my opinion these tests would be  


;NO COLOR CORRECTION FOR EXTINCTION. FOR NOW ...

; IF !dustem_do_cc NE 0 AND !dustem_never_do_cc EQ 0 THEN BEGIN
;   message,'DOING color correction calculations',/info
; ENDIF ELSE BEGIN
;   message,'SKIPPING color correction calculations',/info
; ENDELSE
; 
; ind_polsed=where((*!dustem_data.qsed).filt_names NE 'SPECTRUM',count_polsed)  
;   
; IF count_polsed NE 0 THEN BEGIN
;   filter_names=((*!dustem_data.qsed).filt_names)(ind_polsed)
;   spolsed=dustem_cc(st.polsed.wav,P_spec,filter_names,cc=cc)
;   dustem_polsed[ind_polsed]=spolsed
; ENDIF


;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).
;I believe this to be automatic because the provided dustem spectra are already sampled on a log-log grid.

ind_spec=where((*!dustem_data.qext).filt_names EQ 'SPECTRUM',count_spec) 
IF count_spec NE 0 THEN dustem_polext(ind_spec)=interpol(POLEXT_spec,sti.polext.wav,(((*!dustem_data.qext).wav)(ind_spec)))

out_st=sti

;GENERATING THE INTERPOLATES FOR POLFRAC (dustem_polfrac)
if !run_lin then begin
      		
    
    if not keyword_set(dustem_ext) then dustem_ext = dustem_compute_ext(p_dim,sti,EXT_spec) 
    
    If n_elements(dustem_ext) ne n_elements(dustem_polext) then begin                            
        if n_elements(dustem_ext) gt n_elements(dustem_polext) then begin 
      
            dustem_polext_x = dustem_polext
            dustem_ext_x = dustem_polext ;meaning dustem_sed needs to be modified
      
            nwaves = n_elements(dustem_polext)
      
      
              for i=0L,nwaves-1 do begin
      
                  j=where((*!dustem_data.ext).wav EQ (*!dustem_data.fpolext).wav(i),testwav)
                  if testwav ne 0 then dustem_ext_x(i) = dustem_ext(j(0))                                    
      
              endfor
      
        endif ELSE begin
        
            dustem_polext_x = dustem_ext
            dustem_ext_x = dustem_ext ;meaning dustem_polsed needs to be modified
        
            nwaves = n_elements(dustem_ext)
        
            for i=0L,nwaves-1 do begin
                j=where((*!dustem_data.ext).wav EQ (*!dustem_data.fpolext).wav(i),testwav)
                if testwav ne 0 then dustem_polext_x(i) = dustem_polext(j(0))                                    
            endfor
          
        
        ENDELSE
    
    endif ELSE BEGIN
        
        dustem_polext_x = dustem_polext 
        dustem_ext_x = dustem_ext
    
    ENDELSE
    
    dustem_fpolext = dustem_polext_x/dustem_ext_x 
	       
ENDIF

return, dustem_polext



END