/* Reads a merger tree from an HDF file in the official mt format.
 * You may edit this routine to read in only those fields that you are interested in. */

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

// HDF5
#include "hdf5.h"

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

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

int mt_read_mergertree(void) {

  int ier;
  hid_t gid;

  // First check to see whether the MergerTree group exists and either open it or make it
  gid=mt_gopen(mt_fid,"MergerTree");

  //----------------------------------------------------------------------------------
#ifdef MTTABLE

  hsize_t nfields;

  // Can't think of any simple way of reading the table in one go
  // Read the table length
  NHalo=mt_treadinfo(gid,"Halo",&nfields);
  if (nfields < mttable_nfields) {
    printf("***mt_read_mergertree: nfields < mttable_nfields***\n");
    return -1;
  }
  if (NHalo <= 0) {
    printf("***mt_read_mergertree: NHalo <= 0***\n");
    return -2;
  }

  // Allocate space for the table
  mtTable=(struct mttable_t *)calloc(NHalo,sizeof(*mtTable));

  // Read in the table
  ier=mt_tread(gid,"Halo",
	       mttable_nfields,mttable_fieldnames,mttable_field_types,
	       NHalo,mtTable);

  //----------------------------------------------------------------------------------
#else // Merger tree stored in arrays

  ier=mt_aget_int(gid,".","NHalo",&NHalo);
  //printf("NHalo=%d\n",(int)NHalo);
  // Mass
  mtMass=(float *)calloc(NHalo,sizeof(*mtMass));
  ier=mt_dread(gid,"Mass","float",mtMass);
  // Density
  mtDensity=(float *)calloc(NHalo,sizeof(*mtDensity));
  ier=mt_dread(gid,"Density","float",mtDensity);
  // HaloID
  mtHaloID=(int64_t *)calloc(NHalo,sizeof(*mtHaloID));
  ier=mt_dread(gid,"HaloID","int64",mtHaloID);
  // Snapshot
  mtSnapshot=(int32_t *)calloc(NHalo,sizeof(*mtSnapshot));
  ier=mt_dread(gid,"Snapshot","int32",mtSnapshot);
  // FirstSubhaloIndex
  mtFirstSubhaloIndex=(int32_t *)calloc(NHalo,sizeof(*mtFirstSubhaloIndex));
  ier=mt_dread(gid,"FirstSubhaloIndex","int32",mtFirstSubhaloIndex);
  // NeighbourIndex
  mtNeighbourIndex=(int32_t *)calloc(NHalo,sizeof(*mtNeighbourIndex));
  ier=mt_dread(gid,"NeighbourIndex","int32",mtNeighbourIndex);
  // HostHaloIndex
  mtHostHaloIndex=(int32_t *)calloc(NHalo,sizeof(*mtHostHaloIndex));
  ier=mt_dread(gid,"HostHaloIndex","int32",mtHostHaloIndex);
  // FirstProgenitorIndex
  mtFirstProgenitorIndex=(int32_t *)calloc(NHalo,sizeof(*mtFirstProgenitorIndex));
  ier=mt_dread(gid,"FirstProgenitorIndex","int32",mtFirstProgenitorIndex);
  // NextSiblingIndex
  mtNextSiblingIndex=(int32_t *)calloc(NHalo,sizeof(*mtNextSiblingIndex));
  ier=mt_dread(gid,"NextSiblingIndex","int32",mtNextSiblingIndex);
  // DescendantIndex
  mtDescendantIndex=(int32_t *)calloc(NHalo,sizeof(*mtDescendantIndex));
  ier=mt_dread(gid,"DescendantIndex","int32",mtDescendantIndex);

#endif

  // Now close the open group
  ier=mt_gclose(gid);

  return ier;

}
