FUNCTION dustem_read_gas,file,silent=silent,help=help

;+
; NAME:
;    dustem_read_gas
; PURPOSE:
;    reads the file GAS.dat and returns the corresponding grain structure
; CATEGORY:
;    Dustem
; CALLING SEQUENCE:
;    st=dustem_read_gas(file,[/silent][,/help])
; INPUTS:
;    file = file name to be read
; OPTIONAL INPUT PARAMETERS:
;    None
; OUTPUTS:
;    st: output structure
; OPTIONAL OUTPUT PARAMETERS:
;    None
; ACCEPTED KEY-WORDS:
;    silent      = If set, keeps quiet
;    help      = If set, print this help
; COMMON BLOCKS:
;    None
; SIDE EFFECTS:
;    None
; RESTRICTIONS:
;    The dustem idl wrapper must be installed
; PROCEDURE:
;    None
; EXAMPLES
;    
; MODIFICATION HISTORY:
;    Written by J.-Ph. Bernard, N. Flagey, D. Paradis Jan 2007
;    see evolution details on the dustem cvs maintained at CESR
;    Contact J.-Ph. Bernard (Jean-Philippe.Bernard@cesr.fr) in case of problems.
;-

;=== Format for DUSTEM_WHICH=WEB3p8
;# DUSTEM : Gas quantities from Ysard, Juvela & Verstraete (2011)
;# 
;# Gas temperature (K), hydrogen density, H2 density, CR rate, G0 (if>0 supersedes GRAIN.DAT)
;# number of charge type (electrons, H+, C+, ...)
;# ion density (cm-3), mass (amu), charge, polarizability (A^3)
;# >>>>> 1st line is electron <<<<<<<
;# 
;  3.5550e+03 1e-1 -1e0 5e-17   1.0000e+00 3
;-1e0 5.4858e-04 -1e0 0.00
;0e0  1.00794e0   1e0 0.67
; -1.3000e-05 12.0107e0   1e0 1.54

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

CASE !dustem_which OF
  'WEB3p8': BEGIN
    ;==read and count comments
    openr,unit,file,/get_lun
    str=''
    comments=['']
    first_char='#'
    Ncomments=0L
    WHILE first_char EQ '#' DO BEGIN
       readf,unit,str
       first_char=strmid(str,0,1)
       comments=[comments,str]
       IF first_char EQ '#' THEN Ncomments=Ncomments+1
    ENDWHILE
    comments=comments[1:*]
    IF Ncomments NE 0 THEN comments=comments[0:Ncomments-1]
    close,unit
    free_lun,unit
    st={file:file, $
        Tgas:0.D0,nH:0.D0,nh2:0.D0,CR_rate:0.D0,G0:0.D0,n_charges:0L, $
        charges_props:ptr_new(), $
        comments:comments}
    one_charge_st={density:0.D0,mass:0.D0,charge:0.D0,polarizability:0.D0}
    ;== pass comments
    openr,unit,file,/get_lun
    str=''
    FOR i=0L,Ncomments-1 DO BEGIN
      readf,unit,str
    ENDFOR
    ;== read data
    readf,unit,str
    strv=str_sep(strcompress(strtrim(str,2)),' ') & ii=0L
    st.Tgas=double(strv[ii]) & ii=ii+1
    st.nH=double(strv[ii]) & ii=ii+1
    st.nH2=double(strv[ii]) & ii=ii+1
    st.CR_rate=double(strv[ii]) & ii=ii+1

    ;modified because the file geometry did change !!
    ;st.G0=double(strv[ii]) & ii=ii+1

    st.G0=double(strv[ii]) & ii=ii+1

    ;This for the version distributed by Nathalie
    ;readf,unit,str
    ;strv=str_sep(strcompress(strtrim(str,2)),' ') & ii=0L
    ;st.n_charges=long(strv[ii])
    ;This for the right version, not distributed by Nathalie !

    st.n_charges=long(strv[ii])
    charges_st=replicate(one_charge_st,st.n_charges)
    FOR i=0L,st.n_charges-1 DO BEGIN
      readf,unit,str
      strv=str_sep(strcompress(strtrim(str,2)),' ') & ii=0L
      charges_st[i].density=double(strv[ii]) & ii=ii+1
      charges_st[i].mass=double(strv[ii]) & ii=ii+1
      charges_st[i].charge=double(strv[ii]) & ii=ii+1
      charges_st[i].polarizability=double(strv[ii])
    ENDFOR
    close,unit
    free_lun,unit
    st.charges_props=ptr_new(charges_st)
  END
  'VERSTRAETE':BEGIN
    ;GAS.DAT did not exist in that version. Kept empty
    st=0
  END
  ELSE:BEGIN
    ;GAS.DAT did not exist in that version. Kept empty
    st=0
  END
ENDCASE

;stop 

the_end:

RETURN,st

END