FUNCTION dustem_compute_polext,p_dim,$
                              st=st,$
                              POLEXT_spec=POLEXT_spec,$
                              SPEXT_spec=SPEXT_spec,$
                              dustem_fpolext=dustem_fpolext,$
                              dustem_ext=dustem_ext,$
                              help=help

;+
; NAME:
;    dustem_compute_polext
;
; PURPOSE:
;    Computes Polarized Extinction in units of optical depth for a given Dustem spectrum
;
; CATEGORY:
;    Dustem, Distributed, 
;
; CALLING SEQUENCE:
;    polext=dustem_compute_polext(p_dim[,st=][,polext_spec=][,spext_spec=][,dustem_fpolext=][,dustem_ext=][,/help])
;
; INPUTS:
;    p_dim      = parameter values
;
; OPTIONAL INPUT PARAMETERS:
;    st              = Dustem output structure
;    dustem_ext      = Computed polarized extinction for spectrum points in !dustem_data
;
; OUTPUTS:
;    polext       = computed Polarized Extinction for spectrum points in !dustem_data 
;
; OPTIONAL OUTPUT PARAMETERS:
;    st          = Dustem output structure
;    polext_spec = Dustem polarized extiction output 
;    spext_spec  = Dustem polarized exinction fraction
;    dustem_fpolext  = Computed polarized extinction fraction for spectrum points in !dustem_data 
;
; ACCEPTED KEY-WORDS:
;    help      = If set, print this help
;
; COMMON BLOCKS:
;    None
;
; SIDE EFFECTS:
;    None
;
; RESTRICTIONS:
;    The DustEMWrap IDL code must be installed
;
; MODIFICATION HISTORY:
;    Written by JPB,NF,DP Jan-2007
;    Evolution details on the DustEMWrap gitlab.
;    See http://dustemwrap.irap.omp.eu/ for FAQ and help.  
;-



IF keyword_set(help) THEN BEGIN
  doc_library,'dustem_compute_polext'
  dustem_polext=0.
  goto,the_end
ENDIF

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


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

;Hard-coded test to set to potential zero 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
    
    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 = 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
;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,st.polext.wav,(((*(*!dustem_data).qext).wav)(ind_spec)))

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

the_end:

return, dustem_polext

END