FUNCTION dustem_read_grain,file,silent=silent,key_str=key_str,G0=G0,help=help

;+
; NAME:
;    dustem_read_grain
;
; PURPOSE:
;    reads the file GRAIN.dat and returns the corresponding grain structure
;
; CATEGORY:
;    Dustem, Distributed, Mid-Level, Fortran
;
; CALLING SEQUENCE:
;    st=dustem_read_grain(file,[/silent][,key_str=][,G0=][,/help])
;
; INPUTS:
;    file = file name to be read
;
; OPTIONAL INPUT PARAMETERS:
;    None
;
; OUTPUTS:
;    st: output structure
;
; OPTIONAL OUTPUT PARAMETERS:
;    key_str : a string containing grain keywords
;    G0      : The G0 value (scaling factor for radiation field)
;
; ACCEPTED KEY-WORDS:
;    help      = If set, print this help
;    silent      = If set, keeps quiet
;
; COMMON BLOCKS:
;    None
;
; SIDE EFFECTS:
;    None
;
; RESTRICTIONS:
;    The DustEM fortran code must be installed
;    The DustEMWrap idl code must be installed
;
; PROCEDURE:
;    None
;
; EXAMPLES
;    
; MODIFICATION HISTORY:
;    Written by JPB,NF,DP Jan-2007
;    Evolution details on the DustEMWrap gitlab.
;    See http://dustemwrap.irap.omp.eu/ for FAQ and help.  
;-

;=== Format for DUSTEM_WHICH=VERSTRAETE
;; # DUSTEM: definition of grain populations
;; # for each grain TYPE make sure you have the following files
;; # Q_TYPE.DAT in code_f90/les_QABS
;; # C_TYPE.DAT in code_f90/les_CAPA
;; # if option(MIX) MIX_TYPE.DAT in code_f90/les_DAT
;; # if option(POL) POL_TYPE.DAT in code_f90/les_DAT
;; # if keyword(SIZE) SIZE_TYPE.DAT in code_f90/les_DAT
;; # 
;; # run keywords
;; # G0 scaling factor for radiation field
;; # nr of grain types
;; # grain type -  Mdust/MH - rho(g/cm3) - amin(cm) - amax(cm) - alpha - nsize
;; #
;; SIZE
;; 1.00D+00
;; 3
;; PAH1    4.50D-04  2.25D+00  3.50D-08  7.22D-08  -3.50D+00   10  NONE
;; Gra     2.30D-03  2.25D+00  1.00D-07  5.00D-07  -3.50D+00   20  NONE
;; aSil    7.00D-03  3.30D+00  5.00D-07  1.10D-05  -3.50D+00   20  NONE
;; #

;=== Format for DUSTEM_WHICH=WEB3p8
;; # DUSTEM: definition of grain populations
;; # 
;; # run keywords 
;; # G0 scaling factor for radiation field
;; # grain type, nsize, type keywords, Mdust/MH, rho, amin, amax, alpha/a0 [, at, ac, gamma (ED)] [, au, zeta, eta (CV)]
;; # cgs units
;; #
;; sdist
;; 1.000000
;; PAH0              10 mix-logn           7.8000E-04  2.2400E+00  3.5000E-08  1.2000E-07  6.4000E-08  1.0000E-01
;; PAH1              10 mix-logn           7.8000E-04  2.2400E+00  3.5000E-08  1.2000E-07  6.4000E-08  1.0000E-01
;; amCBEx            15 logn             1.6500E-04  1.8100E+00  6.0000E-08  2.0000E-06  2.0000E-07  3.5000E-01
;; amCBEx            25 plaw-ed      1.4500E-03  1.8100E+00  4.0000E-07  2.0000E-04 -2.8000E+00  1.5000E-05  1.5000E-05  2.0000E+00
;; aSil              25 plaw-ed        7.8000E-03  3.5000E+00  4.0000E-07  2.0000E-04 -3.4000E+00  2.0000E-05  2.0000E-05  2.0000E+00

IF keyword_set(help) THEN BEGIN
  doc_library,'dustem_read_grain'
  full_st=0.
  goto,the_end
ENDIF

;!run_pol now defined in dustem_init.pro
;defsysv, '!run_pol', 0.         ; Global flag to know if polarisation is being run
defsysv,'!run_tls',0.           ; Global flag to know if TLS is being run

CASE !dustem_which OF
  
  'VERSTRAETE':BEGIN
    frmt='(A,D,D,D,D,D,I,A)'
    openr,unit,file,/get_lun
    str=''
    REPEAT BEGIN
      readf,unit,str
      first_char=strmid(str,0,1)
    ENDREP UNTIL first_char NE '#'
    key_str=str
    readf,unit,G0
    readf,unit,nspecies
    one_st={type:'',propmass:0.D0,densgr:0.D0,tmin:0.D0,tmax:0.D0,alpha:0.D0,ndiscr:0,flag:''}
    st=replicate(one_st,nspecies)
    FOR i=0L,nspecies-1 DO BEGIN
      readf,unit,str
      strv=str_sep(strcompress(str),' ')
      st(i).type=strv(0)
      st(i).propmass=double(strv(1))
      st(i).densgr=double(strv(2))
      st(i).tmin=double(strv(3))
      st(i).tmax=double(strv(4))
      st(i).alpha=double(strv(5))
      st(i).ndiscr=fix(strv(6))
      st(i).flag=strv(7)
    ENDFOR
    full_st=st
    close,unit
    free_lun,unit
  END

;   ELSE:BEGIN
;     readcol,file,type,albmax,densgr,ngem,hydro,ioni,silent=silent
;     nlines=n_elements(type)
; 
;     one_st={type:0,albmax:0.,densgr:0.,ngem:0,hydro:0.,ioni:0.}
;     st=replicate(one_st,Nlines)
;     st.type=type
;     st.albmax=albmax
;     st.densgr=densgr
;     st.ngem=ngem
;     st.hydro=hydro
;     st.ioni=ioni
;     full_st=st
;   END
    
  
  
  ELSE: BEGIN 
    frmt='(A,D,D,D,D,D,I,A)'
    ;count # of lines
    openr,unit,file,/get_lun
    Nlines=0L
    str=''
    WHILE not eof(unit) DO BEGIN
      readf,unit,str
      Nlines=Nlines+1
    ENDWHILE
    close,unit
    free_lun,unit
    ;== now read the file
    openr,unit,file,/get_lun
    ncurrent=0L
    REPEAT BEGIN
      readf,unit,str
      ncurrent=ncurrent+1
      first_char=strmid(str,0,1)
    ENDREP UNTIL first_char NE '#'
    key_str=str
    readf,unit,G0 & ncurrent=ncurrent+1
    nspecies=nlines-ncurrent
    one_st={grain_type:'',nsize:0L,type_keywords:'',Mdust_o_MH:0.D0,rho:0.D0,amin:0.D0,amax:0.D0,alpha_o_a0:0.D0, $
            at:0.D0,ac:0.D0,gamma:0.D0, $
            au:0.D0,zeta:0.D0,eta:0.D0}
    st=replicate(one_st,nspecies)
    
    FOR i=0L,nspecies-1 DO BEGIN
      readf,unit,str
      strv=str_sep(strcompress(strtrim(str,2)),' ')
      Nstrv=n_elements(strv)
      ii=0L
      st(i).grain_type=strv(ii) & ii=ii+1
      st(i).nsize=long(strv(ii)) & ii=ii+1
      st(i).type_keywords=strv(ii) & ii=ii+1
      if isa(!dustem_kwords) then st(i).type_keywords=(*!dustem_kwords)(i)
      IF stregex(st(i).type_keywords, 'pol', /bool) THEN !run_pol = 1
      IF stregex(st(i).type_keywords, 'dtls', /bool) THEN !run_tls = 1
      st(i).Mdust_o_MH=double(strv(ii)) & ii=ii+1
      st(i).rho=double(strv(ii))  & ii=ii+1
      st(i).amin=double(strv(ii)) & ii=ii+1
      st(i).amax=double(strv(ii)) & ii=ii+1
      st(i).alpha_o_a0=double(strv(ii)) & ii=ii+1
      IF ii LE Nstrv-1 THEN BEGIN
        st(i).at=double(strv(ii)) & ii=ii+1
      ENDIF
      IF ii LE Nstrv-1 THEN BEGIN
        st(i).ac=double(strv(ii)) & ii=ii+1
      ENDIF
      IF ii LE Nstrv-1 THEN BEGIN
        st(i).gamma=double(strv(ii)) & ii=ii+1
      ENDIF
      IF ii LE Nstrv-1 THEN BEGIN
        st(i).au=double(strv(ii)) & ii=ii+1
      ENDIF
      IF ii LE Nstrv-1 THEN BEGIN
        st(i).zeta=double(strv(ii)) & ii=ii+1
      ENDIF
      IF ii LE Nstrv-1 THEN BEGIN
        st(i).eta=double(strv(ii)) & ii=ii+1
      ENDIF
    ENDFOR
    full_st=st
    close,unit
    free_lun,unit
  END  
    


ENDCASE


the_end:

;stop 

RETURN,full_st

END