FUNCTION read_draine,plot_it=plot_it

;+
; NAME:
;       read_draine
; CALLING SEQUENCE:
;       st=read_draine([/plot_it])
; PURPOSE:
;	Reads Draine's optical properties for bulk material and
;	spheriacl grains
; INPUTS:
;       None
; OPTIONAL INPUT:
;       None
; OUTPUTS:
;       st with the following structure
;       st.gra:              Graphite optical properties vs wavelengths
;       st.sil:              Astronomical Silicates optical properties vs wavelengths
;       st.Sil_spheres:      Astronomical Silicates spherical grains Qabs,Qsca,g vs wavelengths
;       st.Sil_spheres_smth: Same as above, but smoothed in UV vs wavelengths
;       st.Gra_spheres:      Graphite spherical grains Qabs,Qsca,g vs wavelengths
;       st.SiC_spheres:      SiC spherical grains Qabs,Qsca,g vs wavelengths
;       st.PAH_neutral:      Neutral PAH Qabs,Qsca,g vs wavelengths
;       st.PAH_ion:          Ionized PAH Qabs,Qsca,g vs wavelengths
; PROCEDURE AND SUBROUTINE USED
;       read_draine_qabs.pro
; SIDE EFFECTS:
;       Requires the DUSTEM_DRAINE_DATA_DIR to be set
; MODIFICATION HISTORY:
;       See cvs logs
;       Written by JPB Mai-2009
;-

silent=1

;==== Read Graphite dielectric function
one_st={wav_mu:0.,e1_1:0.,e2:0.,Ren_1:0.,Imn:0.}

; ICOMP=18: Draine 2003 graphite E parallel to c    
;  0.0100 = grain radius (micron)
;   20.00 = grain temperature (K)
;     386 = number of wavelengths
;  wave(um)      eps_1-1    eps_2      Re(n)-1     Im(n)

dir=getenv('DUSTEM_DRAINE_DATA_DIR')+'/Gra_diel_func/'
file=dir+'callindex.out_CpaD03_0.01.txt'
readcol,file,wav_mu,e1_1,e2,Ren_1,Imn,silent=silent
Ncol=n_elements(wav_mu)
st=replicate(one_st,Ncol)
st.wav_mu=wav_mu
st.e1_1=e1_1
st.e2=e2
st.Ren_1=Ren_1
st.Imn=Imn

st_gra_CpaD03_0p01=st

dir=getenv('DUSTEM_DRAINE_DATA_DIR')+'/Gra_diel_func/'
file=dir+'callindex.out_CpaD03_0.10.txt'
readcol,file,wav_mu,e1_1,e2,Ren_1,Imn,silent=silent
Ncol=n_elements(wav_mu)
st=replicate(one_st,Ncol)
st.wav_mu=wav_mu
st.e1_1=e1_1
st.e2=e2
st.Ren_1=Ren_1
st.Imn=Imn

st_gra_CpaD03_0p1=st

dir=getenv('DUSTEM_DRAINE_DATA_DIR')+'/Gra_diel_func/'
file=dir+'callindex.out_CpeD03_0.01.txt'
readcol,file,wav_mu,e1_1,e2,Ren_1,Imn,silent=silent
Ncol=n_elements(wav_mu)
st=replicate(one_st,Ncol)
st.wav_mu=wav_mu
st.e1_1=e1_1
st.e2=e2
st.Ren_1=Ren_1
st.Imn=Imn

st_gra_CpeD03_0p01=st

dir=getenv('DUSTEM_DRAINE_DATA_DIR')+'/Gra_diel_func/'
file=dir+'callindex.out_CpeD03_0.10.txt'
readcol,file,wav_mu,e1_1,e2,Ren_1,Imn,silent=silent
Ncol=n_elements(wav_mu)
st=replicate(one_st,Ncol)
st.wav_mu=wav_mu
st.e1_1=e1_1
st.e2=e2
st.Ren_1=Ren_1
st.Imn=Imn

st_gra_CpeD03_0p1=st

st_gra={st_gra_CpaD03_0p01:st_gra_CpaD03_0p01, $  ;C parallel for a=0.01 mic
        st_gra_CpaD03_0p1:st_gra_CpaD03_0p1, $    ;C parallel for a=0.1 mic
        st_gra_CpeD03_0p01:st_gra_CpeD03_0p01, $   ;C orth for a=0.01 mic
        st_gra_CpeD03_0p1:st_gra_CpeD03_0p1}     ;C orth for a=0.1 mic

IF keyword_set(plot_it) THEN BEGIN
  yr=[1e-6,1e4]
  plot,st_gra.ST_GRA_CPAD03_0P01.wav_mu,st_gra.ST_GRA_CPAD03_0P01.e1_1,/xlog,/ylog,yr=yr,/ysty
  oplot,st_gra.ST_GRA_CPAD03_0P1.wav_mu,st_gra.ST_GRA_CPAD03_0P1.e1_1,linestyle=2
  oplot,st_gra.ST_GRA_CPAD03_0P01.wav_mu,st_gra.ST_GRA_CPAD03_0P01.e2,color=100
  oplot,st_gra.ST_GRA_CPAD03_0P1.wav_mu,st_gra.ST_GRA_CPAD03_0P1.e2,color=100,linestyle=2
ENDIF

;stop

;==== Read Astronomical Silicate dielectric function
one_st={wav_mu:0.,e1_1:0.,e2:0.,Ren_1:0.,Imn:0.}

; ICOMP=17: Draine 2003 astrosilicate               
;  0.1000 = grain radius (micron)
;   20.00 = grain temperature (K)
;     837 = number of wavelengths
;  wave(um)      eps_1-1    eps_2      Re(n)-1     Im(n)

dir=getenv('DUSTEM_DRAINE_DATA_DIR')+'/Sil_diel_func/'
file=dir+'callindex.out_silD03.txt'
readcol,file,wav_mu,e1_1,e2,Ren_1,Imn,silent=silent
Ncol=n_elements(wav_mu)
st=replicate(one_st,Ncol)
st.wav_mu=wav_mu
st.e1_1=e1_1
st.e2=e2
st.Ren_1=Ren_1
st.Imn=Imn

st_sil=st

IF keyword_set(plot_it) THEN BEGIN
  tit='Astronomical Silicates'
  xtit='wav (mic)' & ytit='Re(n)-1, Im(n)'
  plot,st_sil.wav_mu,st_sil.Ren_1,/xlog,/ylog,yr=yr,/ysty,tit=tit,xtit=xtit,ytit=ytit
  oplot,st_sil.wav_mu,st_sil.Imn,color=100
ENDIF

;==== Read Spherical particles Qabs

;==== Silicates
dir=getenv('DUSTEM_DRAINE_DATA_DIR')+'/Particles/Sil/'
file=dir+'Sil_21'

st=read_draine_qabs(file)
Nsize=n_elements(st.size)

IF keyword_set(plot_it) THEN BEGIN
  ytit='Qabs/a' & tit='Silicate Spheres'
  FOR i=0L,Nsize-1 DO BEGIN
    IF i EQ 0 THEN BEGIN
      plot,st(i).wav_mu,st(i).qabs/st(i).size,/xlog,/ylog,yr=yr,tit=tit,xtit=xtit,ytit=ytit
    ENDIF ELSE BEGIN
      oplot,st(i).wav_mu,st(i).qabs/st(i).size
    ENDELSE
  ENDFOR
ENDIF

st_part_Si=st

;==== Graphite
;Graphite, 1/3-2/3 approximation, T=25K  

dir=getenv('DUSTEM_DRAINE_DATA_DIR')+'/Particles/Gra/'
file=dir+'Gra_21'

st=read_draine_qabs(file)
Nsize=n_elements(st.size)

IF keyword_set(plot_it) THEN BEGIN
  ytit='Qabs/a' & tit='Graphite Spheres'
  FOR i=0L,Nsize-1 DO BEGIN
    IF i EQ 0 THEN BEGIN
      plot,st(i).wav_mu,st(i).qabs/st(i).size,/xlog,/ylog,yr=yr,tit=tit,xtit=xtit,ytit=ytit
    ENDIF ELSE BEGIN
      oplot,st(i).wav_mu,st(i).qabs/st(i).size
    ENDELSE
  ENDFOR
ENDIF

st_part_Gra=st

;==== SiC

dir=getenv('DUSTEM_DRAINE_DATA_DIR')+'/Particles/SiC/'
file=dir+'SiC_21'

st=read_draine_qabs(file)
Nsize=n_elements(st.size)

IF keyword_set(plot_it) THEN BEGIN
  ytit='Qabs/a' & tit='SiC Spheres'
  FOR i=0L,Nsize-1 DO BEGIN
    IF i EQ 0 THEN BEGIN
      plot,st(i).wav_mu,st(i).qabs/st(i).size,/xlog,/ylog,yr=yr,tit=tit,xtit=xtit,ytit=ytit
    ENDIF ELSE BEGIN
      oplot,st(i).wav_mu,st(i).qabs/st(i).size
    ENDELSE
  ENDFOR
ENDIF

st_part_SiC=st

;==== Smoothed UV Silicates

dir=getenv('DUSTEM_DRAINE_DATA_DIR')+'/Particles/SmthUV_Sil/'
file=dir+'suvSil_21'

st=read_draine_qabs(file)
Nsize=n_elements(st.size)

IF keyword_set(plot_it) THEN BEGIN
  ytit='Qabs/a' & tit='Smoothed UV Silicates Spheres'
  FOR i=0L,Nsize-1 DO BEGIN
    IF i EQ 0 THEN BEGIN
      plot,st(i).wav_mu,st(i).qabs/st(i).size,/xlog,/ylog,yr=yr,tit=tit,xtit=xtit,ytit=ytit
    ENDIF ELSE BEGIN
      oplot,st(i).wav_mu,st(i).qabs/st(i).size
    ENDELSE
  ENDFOR
ENDIF

st_part_SmthUV_Sil=st

;==== Neutral PAH
dir=getenv('DUSTEM_DRAINE_DATA_DIR')+'/Particles/PAH/'
file=dir+'PAHneu_30'

st=read_draine_qabs(file)
Nsize=n_elements(st.size)

IF keyword_set(plot_it) THEN BEGIN
  ytit='Qabs/a' & tit='Neutral PAH'
  FOR i=0L,Nsize-1 DO BEGIN
    IF i EQ 0 THEN BEGIN
      plot,st(i).wav_mu,st(i).qabs/st(i).size,/xlog,/ylog,yr=yr,tit=tit,xtit=xtit,ytit=ytit
    ENDIF ELSE BEGIN
      oplot,st(i).wav_mu,st(i).qabs/st(i).size
    ENDELSE
  ENDFOR
ENDIF

st_part_PAHneu=st

;==== Ionized PAH
dir=getenv('DUSTEM_DRAINE_DATA_DIR')+'/Particles/PAH/'
file=dir+'PAHion_30'

st=read_draine_qabs(file)
Nsize=n_elements(st.size)

IF keyword_set(plot_it) THEN BEGIN
ytit='Qabs/a' & tit='Ionized PAH'
  FOR i=0L,Nsize-1 DO BEGIN
    IF i EQ 0 THEN BEGIN
      plot,st(i).wav_mu,st(i).qabs/st(i).size,/xlog,/ylog,yr=yr,tit=tit,xtit=xtit,ytit=ytit
    ENDIF ELSE BEGIN
      oplot,st(i).wav_mu,st(i).qabs/st(i).size
    ENDELSE
  ENDFOR
ENDIF

st_part_PAHion=st

;=== create final structure
st={gra:st_gra, $
    sil:st_sil, $
    Sil_spheres:st_part_Si, $
    Gra_spheres:st_part_Gra, $
    SiC_spheres:st_part_SiC, $
    Sil_spheres_smth:st_part_SmthUV_Sil, $
    PAH_neutral:st_part_PAHneu, $
    PAH_ion:st_part_PAHion}

RETURN,st

END