#include <stdlib.h>

#define H5_USE_16_API 1
#include "hdf5.h"

#include "cctk.h"
#include "cctk_Parameters.h"
#include "cctk_Arguments.h"
#include "cctk_Functions.h"

// Interface of function for fortran module eostable
void CCTK_FNAME(allocate_eosmodule)
                (CCTK_INT* , CCTK_INT* , CCTK_INT*, CCTK_REAL*,
                 CCTK_REAL*, CCTK_REAL*, CCTK_REAL*, CCTK_REAL*);

// Catch HDF5 errors
#define HDF5_ERROR(fn_call)                                              \
  do {                                                                   \
    int _error_code = fn_call;                                           \
    if (_error_code < 0)                                                 \
      CCTK_VWarn (CCTK_WARN_ABORT, __LINE__, __FILE__, CCTK_THORNSTRING, \
                  "HDF5 call '%s' returned error code %d",               \
                  #fn_call, _error_code);                                \
  } while (0)

static int file_is_readable(const char* filename);
static int file_is_readable(const char* filename)
{
    FILE* fp = NULL;
    fp = fopen(filename, "r");
    if(fp != NULL)
    {
        fclose(fp);
        return 1;
    }
    return 0;
}

// Cactus calls this function. It reads in the table and calls a fortran
// function to setup values for the fortran eos module
void EOS_OMNI_ReadTable(CCTK_ARGUMENTS)
{
  DECLARE_CCTK_PARAMETERS
  DECLARE_CCTK_ARGUMENTS

  CCTK_Info(CCTK_THORNSTRING,"*******************************");
  CCTK_Info(CCTK_THORNSTRING,"Reading nuc_eos table file:");
  CCTK_Info(CCTK_THORNSTRING,nuceos_table_name);
  CCTK_Info(CCTK_THORNSTRING,"*******************************");


  hid_t file;
  if (!file_is_readable(nuceos_table_name))
    CCTK_VWarn(CCTK_WARN_ABORT, __LINE__, __FILE__, CCTK_THORNSTRING,
               "Could not read nuceos_table_name'%s'",
               nuceos_table_name);
  HDF5_ERROR(file = H5Fopen(nuceos_table_name, H5F_ACC_RDONLY, H5P_DEFAULT));

// Use these two defines to easily read in a lot of variables in the same way
// The first reads in one variable of a given type completely
#define READ_EOS_HDF5(NAME,VAR,TYPE,MEM)                                      \
  do {                                                                        \
    hid_t dataset;                                                            \
    HDF5_ERROR(dataset = H5Dopen(file, NAME));                                \
    HDF5_ERROR(H5Dread(dataset, TYPE, MEM, H5S_ALL, H5P_DEFAULT, VAR));       \
    HDF5_ERROR(H5Dclose(dataset));                                            \
  } while (0)
// The second reads a given variable into a hyperslab of the alltables array
#define READ_EOSTABLE_HDF5(NAME,OFF)                                     \
  do {                                                                   \
    hsize_t offset[2]     = {OFF,0};                                     \
    H5Sselect_hyperslab(mem3, H5S_SELECT_SET, offset, NULL, var3, NULL); \
    READ_EOS_HDF5(NAME,alltables,H5T_NATIVE_DOUBLE,mem3);                \
  } while (0)

  // Read size of tables
  CCTK_INT nrho, ntemp, nye;

  READ_EOS_HDF5("pointsrho",  &nrho,  H5T_NATIVE_INT, H5S_ALL);
  READ_EOS_HDF5("pointstemp", &ntemp, H5T_NATIVE_INT, H5S_ALL);
  READ_EOS_HDF5("pointsye",   &nye,   H5T_NATIVE_INT, H5S_ALL);


  // Allocate memory for tables
  CCTK_REAL *alltables, *logrho, *logtemp, *ye, energy_shift;

  if (!(alltables = (CCTK_REAL*)malloc(nrho * ntemp * nye * 19 * sizeof(CCTK_REAL))))
    CCTK_WARN(CCTK_WARN_ABORT, "Cannot allocate memory for EOS table\n");
  if (!(logrho = (CCTK_REAL*)malloc(nrho * sizeof(CCTK_REAL))))
    CCTK_WARN(CCTK_WARN_ABORT, "Cannot allocate memory for EOS table\n");
  if (!(logtemp = (CCTK_REAL*)malloc(ntemp * sizeof(CCTK_REAL))))
    CCTK_WARN(CCTK_WARN_ABORT, "Cannot allocate memory for EOS table\n");
  if (!(ye = (CCTK_REAL*)malloc(nye * sizeof(CCTK_REAL))))
    CCTK_WARN(CCTK_WARN_ABORT, "Cannot allocate memory for EOS table\n");

  // Prepare HDF5 to read hyperslabs into alltables
  hsize_t table_dims[2] = {19, nrho * ntemp * nye};
  hsize_t var3[2]       = { 1, nrho * ntemp * nye};
  hid_t mem3 =  H5Screate_simple(2, table_dims, NULL);

  // Read alltables
  READ_EOSTABLE_HDF5("logpress",  0);
  READ_EOSTABLE_HDF5("logenergy", 1);
  READ_EOSTABLE_HDF5("entropy",   2);
  READ_EOSTABLE_HDF5("munu",      3);
  READ_EOSTABLE_HDF5("cs2",       4);
  READ_EOSTABLE_HDF5("dedt",      5);
  READ_EOSTABLE_HDF5("dpdrhoe",   6);
  READ_EOSTABLE_HDF5("dpderho",   7);
  // chemical potentials
  READ_EOSTABLE_HDF5("muhat",     8);
  READ_EOSTABLE_HDF5("mu_e",      9);
  READ_EOSTABLE_HDF5("mu_p",     10);
  READ_EOSTABLE_HDF5("mu_n",     11);
  // compositions
  READ_EOSTABLE_HDF5("Xa",       12);
  READ_EOSTABLE_HDF5("Xh",       13);
  READ_EOSTABLE_HDF5("Xn",       14);
  READ_EOSTABLE_HDF5("Xp",       15);
  // average nucleus
  READ_EOSTABLE_HDF5("Abar",     16);
  READ_EOSTABLE_HDF5("Zbar",     17);
  // Gamma
  READ_EOSTABLE_HDF5("gamma",    18);


  // Read additional tables and variables
  READ_EOS_HDF5("logrho",       logrho,        H5T_NATIVE_DOUBLE, H5S_ALL);
  READ_EOS_HDF5("logtemp",      logtemp,       H5T_NATIVE_DOUBLE, H5S_ALL);
  READ_EOS_HDF5("ye",           ye,            H5T_NATIVE_DOUBLE, H5S_ALL);
  READ_EOS_HDF5("energy_shift", &energy_shift, H5T_NATIVE_DOUBLE, H5S_ALL);

  HDF5_ERROR(H5Sclose(mem3));
  HDF5_ERROR(H5Fclose(file));

  // Give all values to fortran - which will copy them until I can find out how
  // I can tell Fortran to use those pointers inside the module
  CCTK_FNAME(allocate_eosmodule)
    (&nrho, &ntemp, &nye, alltables, logrho, logtemp, ye, &energy_shift);


  // Free the memory again because fortran copied the whole thing now
  free(ye);
  free(logtemp);
  free(logrho);
  free(alltables);

}