PRO dustem_plot_polsed, st, p_dim, dustem_polsed, aligned=aligned, win=win

IF keyword_set(ps) THEN BEGIN
    set_plot, 'PS'
    device, filename=ps, /color, /encapsulated
ENDIF ELSE BEGIN
  set_plot,'X'
    IF keyword_set(win) then window,win,title='DUSTEM WRAP (POLSED)'
ENDELSE


;I will try to use the refreshed version of the polsed model spectrum
;This is because the default procedures for dustem

titstq='Polarization Intensity'

ytitstq=textoidl('P_\nu (W/m^2/Hz/sr) for N_H=10^{20} H/cm^2')

xtit=textoidl('\lambda (\mum)')
xr=[10,2e4]

yr=[1e-30,1e-19]

Ngrains=(*!dustem_params).Ngrains
Nwaves=(size(st.polsed.em_tot))[1]

if not keyword_set(st) then begin
    ;activate the plugins here 
    st=dustem_run(p_dim)
    dustem_polsed = dustem_compute_polsed(p_dim,st)
ENDIF

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
I = st.sed.em_tot*fact;This is Total intensity I

out=0.

frac=P/I
tes=where(finite(frac) eq 0)
frac(tes)=0.

colors=dustem_grains_colors(Ngrains,ls,ps=ps)
scopes=tag_names((*!dustem_scope))
IF scopes[0] NE 'NONE' THEN BEGIN

FOR i=0L,n_tags(*!dustem_scope)-1 DO BEGIN   
    
    IF total(strsplit((*(*!dustem_scope).(i)),'+',/extract) EQ 'REPLACE_QSED') THEN BEGIN
        
        Q_spec=(*(*!dustem_plugin).(i))[*,1]
    ENDIF 
    
    IF total(strsplit((*(*!dustem_scope).(i)),'+',/extract) EQ 'REPLACE_USED') THEN BEGIN
        U_spec=(*(*!dustem_plugin).(i))[*,2]
    ENDIF
    
ENDFOR    

    
IF ~isa(Q_spec) && ~isa(U_spec) THEN polar_ippsi2iqu,I,Q_spec,U_spec,frac,replicate(0.,Nwaves) 

        
FOR i=0L,n_tags(*!dustem_scope)-1 DO BEGIN
    
    IF total(strsplit((*(*!dustem_scope).(i)),'+',/extract) EQ 'ADD_QSED') THEN BEGIN
        Q_spec+=(*(*!dustem_plugin).(i))[*,1]
    ENDIF
    
    IF total(strsplit((*(*!dustem_scope).(i)),'+',/extract) EQ 'ADD_USED') THEN BEGIN
        U_spec+=(*(*!dustem_plugin).(i))[*,2]
    ENDIF              
    
ENDFOR


P=sqrt(Q_spec^2+U_spec^2)

ENDIF

normpol1=dustem_polsed*fact1
normpol=P

;normpol=P*0.+1. ;and then use this to plot the normalized data

; for i = 0, Ngrains-1 do begin
;  if aligned(i) then cgoplot,st.polsed.wav,st.polsed.(i+1)*fact/normpol,color=colors(i+1),psym=psym,line=ls(i+1) ;investigate this line of the plotting of the different components in polarization
; endfor


ylog=1.

cgplot,st.polsed.wav,P,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)',charsize=1.3

cgoplot,st.polsed.wav,P,color='black'

cgoplot,(*!dustem_data.polsed).wav,(*!dustem_data.polsed).values*fact1,color='Tomato',psym=16,symsize=1,thick=2
err_bar,(*!dustem_data.polsed).wav,(*!dustem_data.polsed).values*fact1,yrms=3.*((*!dustem_data.polsed).sigma)/2.*fact1,color=cgColor('Tomato')

cgoplot,(*!dustem_data.polsed).wav,dustem_polsed*fact1,psym=6,color='red',symsize=2

cgplot,st.polsed.wav,P/P,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,ytickformat='(F6.2)',charsize=1.3    

cgoplot,st.polsed.wav,P/P,color='black'

cgoplot,(*!dustem_data.polsed).wav,(*!dustem_data.polsed).values*fact1/normpol1,color='Tomato',psym=16,symsize=1,thick=2
err_bar,(*!dustem_data.polsed).wav,(*!dustem_data.polsed).values*fact1/normpol1,yrms=3.*((*!dustem_data.polsed).sigma)/2./normpol1*1E-20,color=cgColor('Tomato')

;stop

IF keyword_set(ps) THEN BEGIN
  device,/close
  set_plot,'X'
ENDIF




END