PRO dustem_plot_ext,st,ps=ps,_Extra=extra,help=help,UV=UV,win=win,print_Rv=print_Rv,mergePAH=mergePAH,ylog=ylog,xlog=xlog,donotclose=donotclose,yr=yr,ysty=ysty,xr=xr,xsty=xsty,noerrbars=noerrbars,albedo=albedo,multi=multi

;+
; NAME:
;     dustem_plot_ext
; PURPOSE:
;    Plots a Dustem extinction
; CATEGORY:
;    Dustem
; CALLING SEQUENCE:
;    dustem_plot_ext,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_ext,st,yr=[1.e-6,1.e0],/ysty,xr=[1.,400],/xsty
; MODIFICATION HISTORY:
;
;    Modified by VG with flags /UV,/IR
;
;    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_ext'
  goto,the_end
ENDIF

if not keyword_set(multi) then multi=0

defsysv,'!dustem_data',exists=data_exists

if keyword_set(print_Rv) then begin

	Vband = 0.55
	Bband = 0.44
	tauV = interpol(st.ext.ext_tot,st.ext.wav,Vband)
	tauB = interpol(st.ext.ext_tot,st.ext.wav,Bband)
	print
	print,"E(B-V) model/NH",1.086*(tauB-tauV)/1d21
	print,"Av model",1.086*tauV
	print,"Rv model = ",tauV/(tauB-tauV)
	print,"NH/Av model = ",1e21/(tauv*1.086)
        print
	if data_exists then begin
		tauV = interpol((*!dustem_data.ext).values,(*!dustem_data.ext).wav,Vband)
		tauB = interpol((*!dustem_data.ext).values,(*!dustem_data.ext).wav,Bband)
		print,"E(B-V) data/NH",1.086*(tauB-tauV)/1d21
		print,"Av data",1.086*tauV
		print,"Rv data = ",tauv/(tauB-tauV)
		print,"NH/Av data = ",1e21/(tauv*1.086)
        	print
	endif

endif else if keyword_set(albedo) 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
	
	Ngrains=n_tags(st.sed)-2

        abs_tot=st.ext.abs_grain(0)*0.
        sca_tot=st.ext.sca_grain(0)*0.
        FOR i=0L,Ngrains-1 DO BEGIN
                abs_tot=abs_tot+st.ext.abs_grain(i)
                sca_tot=sca_tot+st.ext.sca_grain(i)
        ENDFOR

        xrange=1./st.ext.wav
	xr=[0,10]
	yr=[0,1]
	xtit=textoidl('1/\lambda (\mum^{-1})')
	ytit=textoidl('albedo')
	plot,xrange,/nodata,xtit=xtit,ytit=ytit,tit=tit,xlog=xlog,ylog=ylog,yr=yr,ys=ys,xr=xr,xsty=xsty
	oplot,xrange,sca_tot/(sca_tot+abs_tot)

	
endif else begin

	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
	
	Ngrains=n_tags(st.sed)-2
	colors=dustem_grains_colors(Ngrains,ls,ps=ps)

;	stop

	ytit=textoidl('\sigma_{ext}/N_H (cm^2/H)') 
	if not keyword_set(ps) then tit='Dustem extinction'
	
	abs_tot=st.ext.abs_grain(0)*0.
	sca_tot=st.ext.sca_grain(0)*0.
	FOR i=0L,Ngrains-1 DO BEGIN
  		abs_tot=abs_tot+st.ext.abs_grain(i)
  		sca_tot=sca_tot+st.ext.sca_grain(i)
	ENDFOR

	if !dustem_show_plot EQ 2 or multi eq 2 then begin
		norm = abs_tot+sca_tot
		yr=[0.,2.]
		ylog = 0
		ys = 1
	endif else begin
		norm = abs_tot * 0. + 1.d21
		if keyword_set(UV) then ylog = 0 else ylog = 1
	endelse
	
	;VG
	if keyword_set(UV) then begin
		xtit=textoidl('1/\lambda (\mum^{-1})')
        	xrange=1./st.ext.wav
		xlog = 0
	endif else begin
		xtit=textoidl('\lambda (\mum)')
        	xrange=st.ext.wav
		xlog = 1
	endelse
	if multi eq 0 then plot,xrange,/nodata,xtit=xtit,ytit=ytit,tit=tit,xlog=xlog,ylog=ylog,yr=yr,ys=ys,xr=xr,xsty=xsty
	if multi eq 1 then plot,xrange,/nodata,xtit='',ytit=ytit,tit=tit,xlog=xlog,ylog=ylog,yr=yr,/ys,xr=xr,xsty=xsty,position=[0.17,0.35,0.95,0.95],xtickformat='(A1)',yticks=3,yminor=10
	if multi eq 2 then plot,xrange,/nodata,xtit=xtit,xlog=xlog,ylog=ylog,ys=ys,xsty=xsty,ytit='Normalized',tit='',xr=xr,yr=[0,2],position=[0.17,0.14,0.95,0.35],/noerase,yticks=2,yminor=5,xticklen=0.1
	;VG
	
	;VG
	; Plot model
	FOR i=0L,Ngrains-1 DO BEGIN
		if i le 1 and keyword_set(mergePAH) then begin
  			if i eq 0 then oplot,xrange,(st.ext.abs_grain(0)+st.ext.sca_grain(0)+st.ext.abs_grain(1)+st.ext.sca_grain(1))/norm,color=colors(1),psym=psym,line=ls(1) 
		endif else begin
  			oplot,xrange,(st.ext.abs_grain(i)+st.ext.sca_grain(i))/norm,color=colors(i+1),psym=psym,line=ls(i+1)
		endelse
	ENDFOR
	oplot,xrange,(abs_tot+sca_tot)/norm,psym=psym
	
	if !dustem_show_plot EQ 2 or multi eq 2 then begin
		norm = interpol(st.ext.ext_tot,st.ext.wav,(*!dustem_data.ext).wav)
	endif else begin
		norm = st.ext.ext_tot * 0. + 1.d21
	endelse
	
	; Plot data
	if data_exists then begin
   		if ptr_valid(!dustem_data.ext) then begin
      			if keyword_set(UV) then xdata = 1./(*!dustem_data.ext).wav else xdata = (*!dustem_data.ext).wav
      			if keyword_set(noerrbars) then oplot,xdata,(*!dustem_data.ext).values/norm,psym=4 $
        		else oploterror,xdata,(*!dustem_data.ext).values/norm,xdata*0.,(*!dustem_data.ext).sigma/norm,psym=4
   		endif
	endif
	
	if keyword_set(ps) and not keyword_set(donotclose) then begin
		device,/close
		set_plot,'X'
	endif

endelse

the_end:

END