PRO readheader,datafile,stars=stars


ON_IOERROR,errorhandler 


IF n_params() EQ 0 THEN BEGIN  ; Skip if file name was passed.
    print,'Select file containing data to read'
    datafile=dialog_pickfile(/read,/must_exist)
    IF datafile EQ '' THEN GOTO,done  ; Cancel.
ENDIF


;--------------------
; Declare parameters.
;--------------------

; Parameters need to be declared before they can be read. Their values
; will be read from file, so they are only used here to set the
; correct data type (e.g. float, longword integer, etc.).
; Note that FORTRAN integers should be declared as IDL longword integers.

; The following parameters are needed by COSMOPLOTTER:
ncdm=0l       ; Number of CDM particles.
nsph=0l       ; Number of SPH particles.
nstar=0l      ; Number of star particles.
nobj=0l       ; Total number of particles (nobj=ncdm+nsph+nstar).
ngrid=0l      ; Grid points per dimension in N-body code.
simboxsize=0. ; Comoving box size, units of h^{-1} Mpc.
redshift=0.   ; Redshift of simulation snapshot.
time=0.       ; Time (Gyr) of simulation snapshot.
omega0=0.     ; Present density parameter in matter (CDM+baryons).
omegab0=0.    ; Present density parameter in baryons.
lambda0=0.    ; Present density parameter in vacuum energy (cosm. constant).
hubble0=0.    ; Present value of Hubble parameter, units of 100 km/(s Mpc).
mcdm=0.       ; Mass of CDM particle, unit: solar mass.
msph=0.       ; Mass of gas/star particle, unit: solar mass.

; These parameters need not all be in your simulation file
; header. However, the ones that are not in your header should be
; computed from the ones that are. At the end of this procedure they
; should all be defined.

;************************************************************
; Declare any other parameters appearing in your header here:

; Part of ibuf1:
;itime=0l     ; Step number.
;itstop=0l    ; Stopping step.
;itdump=0l    ; Dumping interval for backups.
;itout=0l     ; Not used.
;time=0.      ; Time (in code units).
;atime=0.     ; Expansion factor (in code units).
;htime=0.     ; Hubble parameter (in code units).
;dtime=0.     ; Timestep (in code units).
;Est=0.       ; Staring energy (in code units).
;T=0.         ; Kinetic energy (in code units).
;Th=0.        ; Thermal energy (in code units).
;U=0.         ; Potential energy (in code units).
;Radiation=0. ; Radiated energy.
;Esum=0.      ; Energy integral.
;Rsum=0.      ; Total radiated energy.
;cpu=0.       ; CPU Time counter.
;tstop=0.     ; Not used.
;tout=0.      ; Next dump time.
;icdump=0l    ; Position in list of dump times.
;padding=0l   ; Padding for isolated boundary conditions.
;Tlost=0.     ; Lost kinetic energy (isolated boundary conditions only).
;Qlost=0.     ; Lost Thermal energy (isolated boundary conditions only).
;Ulost=0.     ; Lost potential energy (isolated boundary conditions only).

;; Part of ibuf2:
;irun=0l      ; Run number.
;; nobj=0l    ; Total number of particles.
;ngas=0l      ; Number of gas particles.
;ndark=0l     ; Number of dark particles.
;L=0l         ; Grid size (top level).
;intl=0l      ; Interlacing switch.
;nlmx=0l      ; Maximum number of refinement levels.
;perr=0.      ; Maximum 2-body percentage force error.
;dtnorm=0.    ; Timestep mulitiplier.
;sft0=0.      ; Current-day softening (grid units).
;sftmin=0.    ; Minimum softening (grid units).
;sftmax=0.    ; Maximum softening (grid units).
;h100=0.      ; Hubbel parameter.
;box100=0.    ; Boxsize (1/h Mpc).
;zmet=0.      ; Metallicity of gas.
;spc0=0.      ; Average particle spacing.
;lcool=0l     ; Cooling switch.
;rmgas=0.     ; Mass of gas particle.
;rmdark=0.    ; Mass of dark matter particle.
;rmnorm=0.    ; Force normalisation.
;tstart=0.    ; Start time.
;; omega0=0.  ; Omega_0
;xlambda0=0.  ; Lambda_0.
;h0t0=0.      ; H_0t_0.

;; Part of ibuf:
;tout1=fltarr(25)  ; List of desired output times.
;************************************************************


;-----------------------------
; Read and compute parameters.
;-----------------------------

; Open file.
readfile:  ; On I/O-error program execution resumes here.
openr,unit,datafile,/get_lun,/f77_unformatted,/swap_if_little_endian
 
print,'Reading header from file ' + datafile 

; Note that the input buffer arrays contain elements of different types 
; (integers and floats). We therefore need to read the header twice, once
; for both types.
ibuf=lonarr(200)   ; Input buffer.
ibuf1=lonarr(100)
ibuf2=lonarr(100)
readu,unit,ibuf,ibuf1,ibuf2
free_lun,unit

nobj=ibuf2[1]
ngas=ibuf2[2]
ndark=ibuf2[3]
L=ibuf2[4]

openr,unit,datafile,/get_lun,/f77_unformatted,/swap_if_little_endian
ibuf=fltarr(200)   ; Input buffer.
ibuf1=fltarr(100)
ibuf2=fltarr(100)
readu,unit,ibuf,ibuf1,ibuf2

time=ibuf1[4]
atime=ibuf1[5]
h100=ibuf2[12]
box100=ibuf2[13]
msph=ibuf2[17]
mcdm=ibuf2[18]
omega0=ibuf2[21]
xlambda0=ibuf2[22]
h0t0=ibuf2[23]
;*******************************************************************


;--------------------
; Compute parameters.
;--------------------

;*********************************************************************
; If your the parameters in your header are not the ones needed by
; COSMOPLOTTER, or are not in the correct units, you should add some
; statements here to define them.
ngrid=L
simboxsize=box100
redshift=1.0/atime-1.0
; Periodic simulations:
; Code time unit is such that t_0=1.
; Isolated simulations:
; Code time unit is 10^10 years.
Mpc=3.086e24                    ; cm.
year=3.156e7                    ; s.
H0=1.0e7/Mpc                    ; s^{-1}.
tunit=h0t0/H0/h100              ; s.
tunit=tunit/year*1.e-9          ; Gyr.
time=time*tunit
lambda0=xlambda0
hubble0=h100
; Particle masses for periodic b.c. are computed later.
; omegab0 is computed below.
;*********************************************************************


if keyword_set(stars)then begin
;	 check particle types
readu,unit	  ; masses
readu,unit		  ; positions
readu,unit    ; velocities
readu,unit    ; smoothing length
readu,unit    ; temperature
itype = lonarr(nobj,/nozero)
readu,unit,itype
; count types
nn = where(itype eq -2,ngals)
nn = where(itype eq -1,nstars)
nn = where(itype eq 0,ndark)
nn = where(itype eq 1,ngas)
nn = where(itype eq 3,nsn)
print,' # gals  : ',ngals
print,' # stars : ',nstars
print,' # Nsn   : ',nsn
print,' # ndark : ',ndark
print,' # ngas  : ',ngas
endif
free_lun,unit

nsph_sim=nsph
ncdm_sim=ncdm
nstar_sim=nstar
nobj_sim=nobj

aexp0=1.0  ; Expansion factor at present time.
aexp=aexp0/(redshift+1.0)
hubble=hubble(redshift,hubble0,omega0,lambda0)
omega=omega(redshift,omega0,lambda0)
lambda=lambda(redshift,hubble0,omega0,lambda0)
omegab=omegab(redshift,omegab0,omega0,lambda0)


;--------------
; Print header.
;--------------

print
print,'-----------------------------------------------------------------------'
print
print,'Data parameters'
print,'==============='
print
print,format='("Grid points per dimension for N-body code: ",7x,i8)',ngrid
print,format='($,"CDM particles per dimension: ",21x,i8)',ncdm_sim^(1./3.)
IF ncdm EQ 0 THEN print,' (not loaded)' ELSE print,''
IF nstar_sim EQ 0 THEN $
  print,format='($,"Gas particles per dimension: ",21x,i8)',nsph_sim^(1./3.) $
ELSE print,format='($,"Gas particles: ",33x,i10)',nsph_sim
IF nsph EQ 0 AND nsph_sim NE 0 THEN print,' (not loaded)' ELSE print,''
IF nstar_sim NE 0 THEN BEGIN
    print,format='($,"Star particles: ",32x,i10)',nstar_sim
    IF nstar EQ 0 THEN print,' (not loaded)' ELSE print,''
ENDIF
print,format='("Comoving size of the box: ",22x,f10.2," h^{-1} Mpc")', $
  simboxsize
print,format='("Mass per CDM particle: ",25x,e10.2," M_sun")',mcdm
IF nsph_sim NE 0 THEN $
  print,format='("Mass per gas/star particle: ",20x,e10.2," M_sun")',msph


;-------------------------------
; Print cosmological parameters.
;-------------------------------

print
print,format='(a9,17x,a17,a17)','Parameter','Value','Present value'
print,format='(a9,17x,a17,a17)','---------','-----','-------------'
print,format='("Expansion factor:",9x,f17.4,f17.4)',aexp,aexp0
print,format='("Redshift:",17x,f17.4,f17.4)',redshift,0.0
print,format='("Hubble parameter (h):",5x,f17.4,f17.4)',hubble,hubble0
print,format='("Omega_matter:",13x,f17.4,f17.4)',omega,omega0
print,format='("Omega_Lambda:",13x,f17.4,f17.4)',lambda,lambda0
;    print,format='("Omega_b*h^2:",14x,f17.4,f17.4)', $
;      omegab*hubble^2,omegab0*hubble0^2
print,format='("Physical length of box:",3x,f17.4,f17.4," h^{-1} Mpc")', $
  aexp*simboxsize,aexp0*simboxsize
IF nstar NE 0 THEN $
  print,format='("Mass fraction of baryons in stars:",x,f8.4," %")', $
  float(nstar_sim)/float(nsph_sim)*100.0
print


;---------------------
; Some error handling.
;---------------------

GOTO,done  ; Skip error handling if okay.

; When an i/o error is encountered, execution proceeds from here. 
; In this way, the program will continue if the file name contains typos.
errorhandler: 
print,!err_string
print,'Your readdata.pro file may not be correct,'
print,'but I will assume that the filename ' + datafile + ' is wrong...'
print,'Please select the correct file'
datafile=dialog_pickfile(/read,/must_exist)
datafilename=datafile   ; This will go into datapars.
IF datafilename EQ '' THEN return  ; Cancel.
GOTO,readfile

done:

END