FUNCTION dustem_compute_polext ,p_dim,sti,POLEXT_spec,out_st=out_st,dustem_qext,dustem_uext,Q_ext,U_ext,SPEXT_spec,dustem_fpolext,dustem_ext=dustem_ext,_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=dustem_run(p_dim)

ENDIF

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

out_st = sti  

EXT_spec = st.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

    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]
            SPEXT_spec=(*(*!dustem_plugin).(i).spec)[*,3]
            
        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]
            SPEXT_spec+=(*(*!dustem_plugin).(i).spec)[*,3]
            
        ENDIF
      
    ENDFOR

ENDIF    

POLEXT_spec=sqrt(QEXT_spec^2+UEXT_spec^2)

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,st.polext.wav,(((*!dustem_data.qext).wav)(ind_spec)))

out_st=st

;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,st,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