// Example of how to read in a merger tree in HDF5 format

#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <unistd.h>

// HDF5
#include "hdf5.h"

// Load function prototypes
#include "mt.h"

// Define the structure that we want to read in.
#include "mt_data.h"

int main(int argc, char **argv){

  char *file;
  char description[CSIZE] = {""} ;
  int i;

  // Complain if file name not supplied
  if (argc != 2) {
    fprintf(stderr, "Usage: readtree <treefile>\n");
    exit(1);
  }
  file=argv[1];

  /* Open data file
     ier=mt_open(file,ioflags,description)
     Inputs: 
       file - file name
       ioflags - read/write string;
          "r" file must exist and will be opened readonly
	  "rw" file must exist and will be opened for writing
	  "wo" file will be created, over-writing any existing file
	  "w" file will be created, will fail if already exists
       description - comment to be read from or written/appended to file */
  mt_open(file,"r",description);
  printf("Description = %s\n",description);

  // Read attributes
  // There is probably a way to iterate over (unknown) attributes, but not easily
  mt_aget_float(mt_fid,".","BoxsizeMpc",&box);
  printf("BoxsizeMpc=%f'\n",box);
  mt_aget_float(mt_fid,".","H100",&h100);
  printf("H100=%f'\n",h100);
  mt_aget_float(mt_fid,".","OmegaBaryon",&omega_baryon);
  printf("OmegaBaryon=%f'\n",omega_baryon);
  mt_aget_float(mt_fid,".","OmegaCDM",&omega_cdm);
  printf("OmegaCDM=%f'\n",omega_cdm);
  mt_aget_float(mt_fid,".","OmegaLambda",&omega_lambda);
  printf("OmegaLambda=%f'\n",omega_lambda);
  mt_aget_float(mt_fid,".","Sigma8",&sigma8);
  printf("Sigma8=%f'\n",sigma8);
  mt_aget_string(mt_fid,".","Title",title);
  printf("Title=%s'\n",title);

  // Read SnapTable from HDF file
  mt_read_snaptable();
  printf("Snap=\n");
  for(i=0;i<NSnap;i++) printf("%d %f %f %f %f\n",
			      (int)Snap[i].snapnum,
			      (float)Snap[i].a,
			      (float)Snap[i].z,
			      (float)Snap[i].t,
			      (float)Snap[i].t0);
  printf("SnapProp=\n");
  for(i=0;i<NSnapProp;i++) printf("%s %s %s\n",
			      SnapProp[i].Name,
			      SnapProp[i].Description,
			      SnapProp[i].Units);

  // Read MergerTree from HDF file
#ifdef MTTABLE
  printf("***Routines not yet written to read in from MergerTree tables***\n");
#else
  mt_read_mergertree();
  printf("MTHalos=\n");
  for(i=0;i<NHalo;i++) printf("%ld %d %f %f %d %d %d %d %d %d\n",
			      (long)mtHaloID[i],
			      (int)mtSnapshot[i],
			      (float)mtDensity[i],
			      (float)mtMass[i],
			      (int)mtFirstSubhaloIndex[i],
			      (int)mtNeighbourIndex[i],
			      (int)mtHostHaloIndex[i],
			      (int)mtFirstProgenitorIndex[i],
			      (int)mtNextSiblingIndex[i],
			      (int)mtDescendantIndex[i]);
#endif

  // Close file
  mt_close();

  //
  mt_tidy();

  return 0;
}

