FUNCTION dustem_read_dustem_lv,file,silent=silent,help=help

;+
; NAME:
;    dustem_read_dustem_lv
;
; PURPOSE:
;    Reads all Dustem SED output
;
; CATEGORY:
;    DustEMWrap, Distributed, High-Level, User-Convenience
;
; CALLING SEQUENCE:
;    st=dustem_read_dustem_lv,file,[/silent],[/help])
;
; INPUTS:
;    file: file to be read
;
; OPTIONAL INPUT PARAMETERS:
;    None
;
; OUTPUTS:
;    st: Dustem output structure
;      structure:
;      IDL> help,st,/str
;      ** Structure <1733e38>, 7 tags, length=28, data length=28, refs=2:
;      WAV             FLOAT         0.0400000
;      EM_GRAIN_1      FLOAT           0.00000
;      EM_GRAIN_2      FLOAT           0.00000
;      EM_GRAIN_3      FLOAT           0.00000
;      EM_GRAIN_4      FLOAT           0.00000
;      EM_GRAIN_5      FLOAT           0.00000
;      EM_TOT          FLOAT           0.00000
;
; OPTIONAL OUTPUT PARAMETERS:
;    None
;
; ACCEPTED KEY-WORDS:
;    help      = If set, print this help
;    silent    = If set, function is silent
;
; COMMON BLOCKS:
;    None
;
; SIDE EFFECTS:
;    None
;
; RESTRICTIONS:
;    The DustEMWrap IDL code must be installed
;
; PROCEDURE:
;    None
;
; EXAMPLES
;    
; MODIFICATION HISTORY:
;    Written by J.-Ph. Bernard
;    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_read_dustem_lv'
  st=0.
  goto,the_end
ENDIF

message,'Reading file '+file,/info
;=== determine the number of grains
str='                               '
OPENR,unit,file,/get_lun
;==read comments
str='' & first_char='#'
WHILE first_char EQ '#' DO BEGIN
  readf,unit,str
  first_char=strmid(str,0,1)
ENDWHILE
;readf,unit,str
CLOSE,unit
free_lun,unit

str=strcompress(str)
str=strtrim(str,2)
vv=str_sep(str,' ')
Ngrains=vv[0]
;Ngrains=!dustem_params.Ngrains

;Read model output
instruc='readcol,file,wav,'
fstr=strarr(Ngrains) & fstr[*]='D,'
format_str='D,'+strjoin(fstr)+'D'
FOR i=0L,Ngrains-1 DO BEGIN
  instruc=instruc+'em_grain_'+strtrim(i+1,2)+','
ENDFOR
instruc=instruc+"em_tot,silent=silent,format='"+format_str+"'"
message,'Executing '+instruc,/info
toto=execute(instruc)

Nlines=n_elements(wav)

instruc='one_st={wav:0.,'
FOR i=0L,Ngrains-1 DO BEGIN
  instruc=instruc+'em_grain_'+strtrim(i+1,2)+':0.D0,'
ENDFOR
instruc=instruc+'em_tot:0.D0}'
message,'Executing '+instruc,/info
toto=execute(instruc)
st=replicate(one_st,Nlines)

st.wav=wav
FOR i=0L,Ngrains-1 DO BEGIN
  instruc='st.(i+1)='+'em_grain_'+strtrim(i+1,2)
  message,'Executing '+instruc,/info
  toto=execute(instruc)
  ;stop
ENDFOR
st.em_tot=em_tot

IF !dustem_redshift NE 0. THEN BEGIN
   st.wav=st.wav*(1+!dustem_redshift)
ENDIF

;stop

the_end:

RETURN,st

END