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 curve
;
; CATEGORY:
;    Dustem, Plotting, Low-level
;  
; 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 DustEMWrap idl code must be installed
;
; PROCEDURE:
;    None
;
; EXAMPLES
;    dustem_init,/wrap
;    dir_in=getenv('DUSTEM_SOFT_DIR')
;    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:
;    Written by J.-Ph. Bernard
;    Modified by VG with flags /UV,/IR
;    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_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