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

;+
; NAME:
;     dustem_plot_ext
; PURPOSE:
;    Plots a Dustem extinction
; CATEGORY:
;    Dustem
; CALLING SEQUENCE:

;    dustem_plot_ext,st[,/help]

; INPUTS:
;    st        = Dustem model output structure
; OPTIONAL INPUT PARAMETERS:

;    _extra    = extra parameters for the plot routine (currently not passed)

; 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
		previous_device=!d.name

   		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')

	;stop

	cgplot,xrange,/nodata,xtit=xtit,ytit=ytit,tit=tit,xlog=xlog,ylog=ylog,yr=yr,ys=ys,xr=xr,xsty=xsty
	cgoplot,xrange,sca_tot/(sca_tot+abs_tot)

ENDIF ELSE BEGIN

	IF multi NE 2 THEN BEGIN
		IF keyword_set(ps) THEN BEGIN
			previous_device=!d.name
   			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)
	Ngrains=(*!dustem_params).Ngrains

	colors=['dark',dustem_grains_colors(Ngrains,ls,/cgplot,ps=ps)]   ;dark added because first color taken for total extinction

	use_lines=replicate(0,Ngrains)

;	stop

	ytit=textoidl('\sigma_{ext}/N_H (cm^2/H)') 

	IF NOT keyword_set(ps) THEN tit='DUSTEMWrapper 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 cgplot,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 cgplot,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 cgplot,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) 

  			IF i EQ 0 THEN cgoplot,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)
  			cgoplot,xrange,(st.ext.abs_grain(i)+st.ext.sca_grain(i))/norm,color=colors(i+1),psym=psym,line=ls(i+1)

		ENDELSE

	ENDFOR

	;cgoplot,xrange,(abs_tot+sca_tot)/norm,psym=psym
	
	;cgoplot,xrange,(abs_tot+sca_tot+powerlaw)/norm,psym=psym
	
	cgoplot,xrange,(st.ext.ext_tot)/norm,psym=psym
	
	

; 	; LINES OF CODE TO PLOT THE NEW ADDED PLUGINS TO THE EXTINCTION WILL HAVE TO BE ADDED HERE 
; 	
; 	cgoplot,xrange,(*(*!dustem_plugin).powerlaw)/norm,psym=psym,color='cadetblue'  

	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 $

      			IF keyword_set(noerrbars) THEN BEGIN
      				cgoplot,xdata,(*!dustem_data.ext).values/norm,psym=4
      			ENDIF ELSE BEGIN
        		    ;else oploterror,xdata,(*!dustem_data.ext).values/norm,xdata*0.,(*!dustem_data.ext).sigma/norm,psym=4
        			;oploterror,xdata,(*!dustem_data.ext).values/norm,xdata*0.,(*!dustem_data.ext).sigma/norm,psym=4
        			cgoplot,xdata,(*!dustem_data.ext).values/norm,psym=4
        			yy=(*!dustem_data.ext).values/norm
        			rms=3.*(*!dustem_data.ext).sigma/norm/2.
        			cgerrplot,xdata,yy-rms,yy+rms,color='black'
        	    ENDELSE

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

		set_plot,'X'

	ENDIF

ENDELSE

the_end:

END