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