PRO dustem_plot_polar,st,ps=ps,_Extra=extra,help=help,UV=UV,SED=SED,Pfrac=Pfrac,print_ratio=print_ratio,win=win,cont=cont,xr=xr,yr=yr,donotclose=donotclose,aligned=aligned,noerrbars=noerrbars,multi=multi,almabands=almabands,nsmooth=nsmooth

;+
; NAME:
;     dustem_plot_polar
; PURPOSE:
;    Plots a Dustem polarization curve
; CATEGORY:
;    Dustem
; CALLING SEQUENCE:
;    dustem_plot_curve,st,,/help][,_extra=]
; INPUTS:
;    st        = Dustem model output structure
; OPTIONAL INPUT PARAMETERS:
;    _extra    = extra parameters for the plot routine
; OUTPUTS:
;    None
; OPTIONAL OUTPUT PARAMETERS:
;    None
; ACCEPTED KEY-WORDS:
;    help      = If set, print this help
; COMMON BLOCKS:
;    None
; SIDE EFFECTS:
;    model is plotted
; RESTRICTIONS:
;    The dustem idl wrapper must be installed
; PROCEDURE:
;    None
; EXAMPLES
;    dustem_init,/wrap
;    dir_in=getenv('DUSTEM_SOFT_DIR')+'/src/dustem3.8_web/'
;    dir=getenv('DUSTEM_RES')
;    st=dustem_read_all_res(dir_in)
;    loadct,13
;    dustem_plot_curve,st,yr=[1.e-6,1.e0],/ysty,xr=[1.,400],/xsty
; MODIFICATION HISTORY:
;
;    Adapted by VG from dustem_plot_ext.pro
;
;    Written by J.-Ph. Bernard
;    see evolution details on the dustem cvs maintained at CESR
;    Contact J.-Ph. Bernard (Jean-Philippe.Bernard@cesr.fr) in case of problems.
;-

IF keyword_set(help) THEN BEGIN
  doc_library,'dustem_plot_polar'
  goto,the_end
ENDIF

if not keyword_set(nsmooth) then nsmooth = 1

nsmooth = 4
almabands = 1.

if not keyword_set(multi) then multi=0

if multi ne 2 then begin
	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
	endelse
endif

;stop

use_col_sed_filt = 70
use_col_data_filt = 250

;Ngrains=n_tags(st.sed)-2
Ngrains=(*!dustem_params).Ngrains

colors=dustem_grains_colors(Ngrains,ls,ps=ps)

xtit=textoidl('\lambda (\mum)')

;ext_tot=st.ext.abs_grain(0)*0.
;polext_tot=st.polext.abs_grain(0)*0.
;FOR i=0L,Ngrains-1 DO BEGIN
  ;polext_tot = polext_tot + (st.polext.abs_grain(i) + st.polext.sca_grain(i))
  ;ext_tot    = ext_tot    + (st.ext.abs_grain(i)    + st.ext.sca_grain(i))
;ENDFOR

; Are observational data defined ?
defsysv,'!dustem_data',exists=data_exists

IF ptr_valid(!dustem_data.polext) THEN BEGIN
   if keyword_set(print_ratio) then begin

		; Correction de couleur
		cc = 1.11

	        Bband = 0.44
	        Vband = 0.55
		Rband = 0.66
		Iband = 0.8
		Jband = 1.25
		Hband = 1.6
		Kband = 2.2
		; pV and tauV are In units of 1e-21 cm2/H 
	        pV = interpol(st.polext.ext_tot,st.polext.wav,Vband)
	        pR = interpol(st.polext.ext_tot,st.polext.wav,Rband)
	        pI = interpol(st.polext.ext_tot,st.polext.wav,Iband)
	        pJ = interpol(st.polext.ext_tot,st.polext.wav,Jband)
	        pH = interpol(st.polext.ext_tot,st.polext.wav,Hband)
	        pK = interpol(st.polext.ext_tot,st.polext.wav,Kband)
		j=where(st.polext.wav gt 5 and st.polext.wav lt 20)
		pSil=max(st.polext(j).ext_tot,jmax)
		Silband = st.polext(j(jmax)).wav
		j=where(st.polext.wav gt 0.4 and st.polext.wav lt 0.75)
		pmeanoptical = mean(st.polext(j).ext_tot)

		pV_data = interpol((*!dustem_data.polext).values,(*!dustem_data.polext).wav,Vband)
	        tauV = interpol(st.ext.ext_tot,st.ext.wav,Vband)
		tauV_data = interpol((*!dustem_data.ext).values,(*!dustem_data.ext).wav,Vband)
	        tauB = interpol(st.ext.ext_tot,st.ext.wav,Bband)
		tauB_data = interpol((*!dustem_data.ext).values,(*!dustem_data.ext).wav,Bband)
	    	print


		; I353 and P353 are in fact nu_Pnu, in unit of erg/s/H, not erg/s/sr/H
	                                ; 1 W/m2/sr/1d20H = 1e-17 erg/s/sr/H
	endif
	IF ptr_valid(!dustem_data.polsed) THEN BEGIN

	    dustem_polsed = dustem_compute_polsed(p_dim,st.polsed,cont=cont)
		x=min(((*!dustem_data.polsed).wav-850)^2,jmin)
		P353 = dustem_polsed(jmin)
		P353_data = (*!dustem_data.polsed).values(jmin)

	        dustem_sed = dustem_compute_sed(p_dim,st,cont=cont)
		x=min(((*!dustem_data.sed).wav-850)^2,jmin)
		I353 = dustem_sed(jmin)
		I353_data = (*!dustem_data.sed).values(jmin)

		; Data is Inu, not nInu, in MJy/sr/1e20
		; Convert data, which is nu_Inu SED in MJy/sr/1d20H into nu_Inu SED model in erg/s/H 
		fact = (4*!pi) * 1e-17 / 1e20 * 353e9

		print
		print,"Observational  constraint"
		print,"-------------------------"
		print,"I353/NH data        ",I353_data," MJy/sr for NH=1e20"
		print,"P353/NH data        ",P353_data," MJy/sr for NH=1e20"
	        print,"P353/I353 data      ",P353_data/I353_data*100," %"
	        print
	        print,"Model results"
		print,"-------------"
	        print,"I353/NH model       ",I353," MJy/sr for NH=1e20"
		print,"P353/NH model       ",P353," MJy/sr for NH=1e20"
		print,"P353/I353 model     ",P353/I353*100," %"

	    IF keyword_set(print_ratio) THEN BEGIN
		    print
			print,"Observational  constraint"
			print,"-------------------------"

	        print,"tau_V data      ",tauV_data
	        print,"(p/tau)_V data      ",pV_data/tauV_data*100," %"
	        print,"p_V/A_V data        ",pV_data/(tauV_data*1.086)*100," %"
			print,"with p_V data       ",pV_data," %"
	        print,"Rs/v data           ",(P353_data/I353_data)/(pv_data/tauv_data)
	        print,"RP/p data           ",(P353_data/1e20) / (pv_data*1e-21)," MJy/sr"
	        print,"I353/Av data        ",(I353_data/1e20) / (1.086*tauv_data*1e-21)," MJy/sr"
			print
	        print,"(p/tau)_V model     ",pV/tauV*100," %"
	        print,"p_V/A_V model       ",pV/(tauV*1.086)*100," %"
	        print,"p_V/E(B-V) model       ",pV/((tauB-tauV)*1.086)*100," %"
			print,"with p_V model      ",pV," %"
	        print,"Rs/v model          ",P353/I353/(pV/tauV)
	        print,"Rs/EBV model          ",P353/I353/(pV/(tauB-tauV))
	        print,"RP/p model          ",(P353/1e20) / (pv*1e-21)," MJy/sr"
	        print,"- R band = 0.66 um  ",(P353/1e20) / (pR*1e-21)," MJy/sr"
	        print,"- I band = 0.8 um   ",(P353/1e20) / (pI*1e-21)," MJy/sr"
	        print,"- J band = 1.25 um  ",(P353/1e20) / (pJ*1e-21)," MJy/sr"
	        print,"- H band = 1.6 um   ",(P353/1e20) / (pH*1e-21)," MJy/sr"
	        print,"- K band = 2.2 um   ",(P353/1e20) / (pK*1e-21)," MJy/sr"
	        print,"- "+string(Silband,format='(F5.2)')+" um band       ",(P353/1e20) / (pSil*1e-21)," MJy/sr"
	        print,"- mean optical      ",(P353/1e20) / (pmeanoptical*1e-21)," MJy/sr"
			print,"I353/Av model       ",(I353/1e20) / (1.086*tauv*1e-21)," MJy/sr"
			print
			;print,"Mixt results"
			;print,"-------------"
		        ;print,"(pV_model/tauV_data) 	",pV/tauV_data*100," %"
			;print,"I353_model/AV_data	",(I353/1e20) / (1.086*tauv_data*1e-21)," MJy/sr"
		        ;print,"P353model/Idata	        ",P353/I353_data*100," %"
		        ;print,"Rs/v P353model/I353data/(pVmodel/tauVdata)	",P353/I353_data/(pV/tauV_data)
		        ;print,"RP/p P353model/pVdata	",(P353/1e20)/(pV_data*1e-21)," MJy/sr"
		        ;print
     	ENDIF
    endif    

endif else if keyword_set(UV) then begin

	if not keyword_set(ps) then tit='Dustem polarization cross-section'
	ytit=textoidl('\sigma_{pol}/N_H (cm^2/H)')  
	if ptr_valid(!dustem_data.polext) then begin
	if !dustem_show_plot EQ 2 or multi eq 2 then begin
		;normpol = polext_tot 
		normpol = st.polext.ext_tot 
		dustem_polext = interpol(st.polext.ext_tot,st.polext.wav,(*!dustem_data.polext).wav)  
		yr=[0.,2.]
		ylog = 0
	endif else begin
		dustem_polext = (*!dustem_data.polext).wav * 0. + 1.d21
		normpol = st.polext.ext_tot * 0. + 1.d21
		ylog = 1
	endelse

	if multi eq 0 then plot_oo,st.polext.wav,st.polext.ext_tot/normpol,/nodata,_Extra=extra,xtit=xtit,ytit=ytit,tit=tit,xr=xr,yr=yr,ylog=ylog
	if multi eq 1 then plot_oo,st.polext.wav,st.polext.ext_tot/normpol,/nodata,_Extra=extra,xtit='',ytit=ytit,tit=tit,xr=xr,yr=yr,ylog=ylog,position=[0.17,0.35,0.95,0.95],xtickformat='(A1)' 
    if multi eq 2 then plot_oo,st.polext.wav,st.polext.ext_tot/normpol,/nodata,_Extra=extra,xtit=xtit,ytit='Normalized',tit='',xr=xr,yr=[-1,2],ylog=ylog,position=[0.17,0.14,0.95,0.35],/noerase,yticks=3,ymino=5,xticklen=0.1
    

	FOR i=0L,Ngrains-1 DO if aligned(i) then oplot,st.polext.wav,(st.polext.abs_grain(i)+st.polext.sca_grain(i))/normpol,color=colors(i+1),psym=psym,line=ls(i+1)

	if total(aligned gt 0) then oplot,st.polext.wav,st.polext.ext_tot/normpol,psym=psym

	if data_exists then if ptr_valid(!dustem_data.polext) then begin
        plotsym,0,/FILL
 		oplot,  (*!dustem_data.polext).wav,(*!dustem_data.polext).values/dustem_polext,ps=4
 		if not keyword_set(noerrbars) then err_bar,(*!dustem_data.polext).wav,(*!dustem_data.polext).values/dustem_polext,yrms=(*!dustem_data.polext).sigma/dustem_polext
	endif
	endif
                                ;	; Fit in (pmax,lmax,K)
;	x=st.polext.wav
;	logx=alog(x)
;
;	logy=alog(polext_tot)
;	j=where(polext_tot gt 0 and x gt 0.2 and x lt 1)
;	res = poly_fit(logx(j),logy(j),2)
;
;	print
;	print,"Serkowski Fit with free K"
;	print,"-------------------------"
;	K    = -res(2)
;	lmax = exp(res(1)/(2*K))
;	pmax = exp(res(0)+K*alog(lmax)^2)
;	print,"* lmax = ",lmax," micron"
;	print,"* K    = ",K
;	print,"* pmax = ",pmax," for NH=1e21"
;	print
;	; overplot the fit
;	;oplot,x(j),exp(res(0)+logx(j)*res(1)-K*logx(j)^2),col=250
;
;	K = 1.15
;	print,"Serkowski Fit with fixed K = "+string(K,format='(F4.2)')
;	print,"-------------------------"
;	res = poly_fit(logx(j),logy(j)+K*logx(j)^2,1)
;	lmax = exp(res(1)/(2*K))
;	pmax = exp(res(0)+K*alog(lmax)^2)
;	print,"* lmax = ",lmax," micron"
;	print,"* pmax = ",pmax," for NH=1e21"
;	print
;
;	; overplot the fit
;	;oplot,x(j),exp(res(0)+logx(j)*res(1)-K*logx(j)^2),col=150

endif else if  ptr_valid(!dustem_data.polsed) then begin
;stop
    ;device, filename = 'plot008.eps', /color, /encapsulated
    if not keyword_set(ps) then tit='Polarization intensity in emission'
    ytit='Polarization intensity'
	ytit=textoidl('\nuP_\nu/N_H (W/m^2/sr/H)') 

	; Convert from erg/s/H into W/m2/sr/H
	fact=1./(4.*!pi)*1.e17/1d20

        ;dustem_polsed = dustem_compute_polsed(p_dim,st.polsed,cont=cont)*fact
	if not keyword_set(xr) then xr=[5e2,5e3]
	if !dustem_show_plot EQ 2 or multi eq 2 then begin
		normpol = st.polsed.(Ngrains+1)*fact
		yr=[0.,2.5]
		ylog = 0
	endif else begin
		normpol = st.polsed.(Ngrains+1)*0. + 1.
        dustem_polsed = ((*!dustem_data.polsed).values)*0. + 1
		ylog=1
		if not keyword_set(yr) then yr=[1e-34,1e-30]
	endelse

    ;if multi eq 0 then plot_oo,st.polsed.wav,/nodata,st.polsed.(Ngrains+1)*fact/normpol,xtit=xtit,ytit=ytit,tit=tit,ylog=ylog,/xlog,xr=xr,yr=yr,/ys,/xs
    ;if multi eq 1 then plot_oo,st.polsed.wav,/nodata,st.polsed.(Ngrains+1)*fact/normpol,xtit='',ytit=ytit,tit=tit,ylog=ylog,/xlog,xr=xr,yr=yr,/ys,/xs,position=[0.17,0.35,0.95,0.95],xtickformat='(A1)' 
    ;if multi eq 2 then plot_oo,st.polsed.wav,/nodata,st.polsed.(Ngrains+1)*fact/normpol,xtit=xtit,ytit='Normalized',tit='',/xlog,xr=xr,/ys,/xs,yr=[0,2],ylog=0,position=[0.17,0.14,0.95,0.35],/noerase,yticks=2,ymino=5,xticklen=0.1
    if multi eq 0 then cgplot,st.polsed.wav,/nodata,st.polsed.(Ngrains+1)*fact/normpol,xtit=xtit,ytit=ytit,tit=tit,ylog=ylog,/xlog,xr=xr,yr=yr,/ys,/xs
    if multi eq 1 then cgplot,st.polsed.wav,/nodata,st.polsed.(Ngrains+1)*fact/normpol,xtit='',ytit=ytit,tit=tit,ylog=ylog,/xlog,xr=xr,yr=yr,/ys,/xs,position=[0.17,0.35,0.95,0.95],xtickformat='(A1)' 
    if multi eq 2 then cgplot,st.polsed.wav,/nodata,st.polsed.(Ngrains+1)*fact/normpol,xtit=xtit,ytit='Normalized',tit='',/xlog,xr=xr,/ys,/xs,yr=[0,2],ylog=0,position=[0.17,0.14,0.95,0.35],/noerase,yticks=2,ymino=5,xticklen=0.1

    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)
	endfor
	;stop
    if (total(aligned) gt 0) then cgoplot,st.polsed.wav,st.polsed.(Ngrains+1)*fact/normpol

    if data_exists then if ptr_valid(!dustem_data.polsed) then begin

		ind_filt=where((*!dustem_data.polsed).filt_names NE 'SPECTRUM',count_filt)
        	;=== Plot the data
        	IF count_filt NE 0 THEN BEGIN

				if !dustem_show_plot EQ 2 or multi eq 2 then begin
				endif else begin
				; Convert data which Pnu in MJy/sr/1d20 into nu*Pnu in W/m2/sr/1d20
     				fact = (3.e10/(1.e-4*(*!dustem_data.polsed).wav(ind_filt)))/1.d40
				endelse
        		dustem_polsed = dustem_compute_polsed(p_dim,st,cont=cont)*fact
				if !dustem_show_plot EQ 2 or multi eq 2 then begin
        			normpol = dustem_polsed
				endif else begin
        			normpol = dustem_polsed * 0. + 1
				endelse
            	plotsym,0,/FILL
            	cgoplot,((*!dustem_data.polsed).wav)(ind_filt),((*!dustem_data.polsed).values)(ind_filt)*fact/normpol(ind_filt),psym=8,_extra=extra,color=use_col_sed_filt,symsize=0.7
            	if not keyword_set(noerrbars) then err_bar,((*!dustem_data.polsed).wav)(ind_filt),((*!dustem_data.polsed).values)(ind_filt)*fact/normpol(ind_filt),yrms=3.*((*!dustem_data.polsed).sigma)(ind_filt)/2.*fact/normpol(ind_filt),color=use_col_sed_filt

        		; PLot color corrected POLSED
            	if keyword_set(ps) then plotsym,8,thick=8 else plotsym,8
        		if multi eq 0 then cgoplot,((*!dustem_data.polsed).wav)(ind_filt),dustem_polsed(ind_filt)/normpol(ind_filt),psym=8,syms=2,_extra=extra,col=use_col_data_filt

        endif
	endif

endif else if keyword_set(Pfrac) then begin

        ; Also show prediction for polarized SED
        if not keyword_set(ps) then tit='Polarization fraction in emission'
        ytit=textoidl('P/I')

        plot_oi,st.polsed.wav,st.polsed.(Ngrains+1)/st.sed.(Ngrains+1),/nodata,_extra=extra,xtit=xtit,ytit=ytit,tit=tit,/xlog,xr=xr,yr=yr
        if (total(aligned) gt 0) then oplot,st.polsed.wav,st.polsed.(Ngrains+1)/st.sed.(Ngrains+1)
        for i = 0, Ngrains-1 do if aligned(i) then oplot,st.polsed.wav,st.polsed.(i+1)/st.sed.(Ngrains+1),color=colors(i+1),psym=psym,linestyle=ls(i+1)
        ;for i = 0, Ngrains-1 do if aligned(i) then begin
	;	x = st.polsed.wav
	;	y = st.polsed.(i+1)/st.sed.(i+1)
	;	j=where(y eq y)
	;	;oplot,st.polsed.wav,smooth(st.polsed.(i+1)/st.sed.(i+1),nsmooth),color=colors(i+1),psym=psym,linestyle=ls(i+1)
	;	oplot,x(j),smooth(y(j),nsmooth),color=colors(i+1),psym=psym,linestyle=ls(i+1)
	;endif
	;if yr(0) lt 0 then oplot,xr,xr*0,linestyle=1

	if keyword_Set(almabands) then begin
		oplot,870*[1,1],yr,linestyle=1
		oplot,1300*[1,1],yr,linestyle=1
		oplot,3100*[1,1],yr,linestyle=1
	endif
        if data_exists then begin

		if ptr_valid(!dustem_data.polfrac) then begin

           		;=== Plot the data
           		ind_filt_polfrac=where((*!dustem_data.polfrac).filt_names NE 'SPECTRUM',count_filt_polfrac)
           		IF count_filt_polfrac NE 0 THEN BEGIN
              			plotsym,0,/FILL
              			oplot,((*!dustem_data.polfrac).wav)(ind_filt_polfrac),((*!dustem_data.polfrac).values)(ind_filt_polfrac),psym=8,_extra=extra,color=use_col_sed_filt
              			if not keyword_set(noerrbars) then err_bar,((*!dustem_data.polfrac).wav)(ind_filt_polfrac),((*!dustem_data.polfrac).values)(ind_filt_polfrac),yrms=2*((*!dustem_data.polfrac).sigma)(ind_filt_polfrac),color=use_col_sed_filt
				print,((*!dustem_data.polfrac).wav)(ind_filt_polfrac),((*!dustem_data.polfrac).values)(ind_filt_polfrac)
	   		endif

        		; PLot color corrected POLFRAC
			; VG : P & I have almost identical color correction => no effect on P/I
			ind_filt=where((*!dustem_data.polsed).filt_names NE 'SPECTRUM',count_filt)
                        if keyword_set(ps) then plotsym,8,thick=8 else plotsym,8

        		dustem_polsed 	= dustem_compute_polsed(p_dim,st.polsed)
        		dustem_sed_tmp 	= dustem_compute_sed(p_dim,st,cont=cont)
			dustem_sed = dustem_polsed
			for i = 0, count_filt-1 do begin
				j=where((*!dustem_data.polsed).filt_names(ind_filt(i)) eq (*!dustem_data.sed).filt_names,count)
				if count eq 1 then dustem_sed(i) = dustem_sed_tmp(j(0))
			endfor 

           	endif

        endif

endif else begin

	if not keyword_set(ps) then tit='Polarization fraction in extinction'
	ytit='Polarization fraction'
	plot_oi,st.polext.wav,st.polext.ext_tot/st.ext.ext_tot,/nodata,_Extra=extra,xtit=xtit,ytit=ytit,tit=tit,xr=xr,yr=yr
	FOR i=0L,Ngrains-1 DO BEGIN
  		oplot,st.polext.wav,smooth((st.polext.abs_grain(i)+st.polext.sca_grain(i))/(st.ext.abs_grain(i)+st.ext.sca_grain(i)),nsmooth),color=colors(i+1),psym=psym,line=ls(i+1)
	ENDFOR
	if (total(aligned) gt 0) then oplot,st.polext.wav,st.polext.ext_tot/st.ext.ext_tot,psym=psym
	oplot,st.polsed.wav,st.polsed.(Ngrains+1)/st.sed.(Ngrains+1),psym=psym,linestyle=2

        if data_exists and ptr_valid(!dustem_data.polext) and ptr_valid(!dustem_data.ext) then begin
        	values_ext_for_polext = interpol((*!dustem_data.ext).values,(*!dustem_data.ext).wav,(*!dustem_data.polext).wav)
        	sigma_ext_for_polext  = interpol((*!dustem_data.ext).sigma ,(*!dustem_data.ext).wav,(*!dustem_data.polext).wav)
		if not keyword_set(noerrbars) then err_bar,(*!dustem_data.polext).wav,(*!dustem_data.polext).values/values_ext_for_polext, $
	 	   yrms=1./values_ext_for_polext * (*!dustem_data.polext).sigma + (*!dustem_data.polext).values/values_ext_for_polext^2 * sigma_ext_for_polext
        endif

endelse

if keyword_set(ps) and not keyword_set(donotclose) then begin
        device,/close
        set_plot,'X'
endif

the_end:

END