/*******************************************************************************
 * This file is part of SWIFT.
 * Copyright (c) 2020 Mladen Ivkovic (mladen.ivkovic@hotmail.com)
 *
 * This program is free software: you can redistribute it and/or modify
 * it under the terms of the GNU Lesser General Public License as published
 * by the Free Software Foundation, either version 3 of the License, or
 * (at your option) any later version.
 *
 * This program is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.
 *
 * You should have received a copy of the GNU Lesser General Public License
 * along with this program.  If not, see <http://www.gnu.org/licenses/>.
 *
 ******************************************************************************/
#ifndef SWIFT_RT_IO_GEAR_H
#define SWIFT_RT_IO_GEAR_H

#define RT_LABELS_SIZE 10

/**
 * @file src/rt/GEAR/rt_io.h
 * @brief Main header file for GEAR M1 Closure radiative transfer
 * scheme IO routines.
 */

/**
 * @brief Specifies which particle fields to read from a dataset
 *
 * @param parts The particle array.
 * @param list The list of i/o properties to read.
 *
 * @return Returns the number of fields to read.
 */
INLINE static int rt_read_particles(const struct part *parts,
                                    struct io_props *list) {

  /* List what we want to read */

  char fieldname[30];
  int count = 0;
  for (int phg = 0; phg < RT_NGROUPS; phg++) {
    sprintf(fieldname, "PhotonEnergiesGroup%d", phg + 1);
    list[count++] =
        io_make_input_field(fieldname, FLOAT, 1, OPTIONAL, UNIT_CONV_ENERGY,
                            parts, rt_data.radiation[phg].energy_density);
    sprintf(fieldname, "PhotonFluxesGroup%d", phg + 1);
    list[count++] = io_make_input_field(fieldname, FLOAT, 3, OPTIONAL,
                                        UNIT_CONV_ENERGY_VELOCITY, parts,
                                        rt_data.radiation[phg].flux);
  }

  list[count++] = io_make_input_field("MassFractionHI", FLOAT, 1, OPTIONAL,
                                      UNIT_CONV_NO_UNITS, parts,
                                      rt_data.tchem.mass_fraction_HI);
  list[count++] = io_make_input_field("MassFractionHII", FLOAT, 1, OPTIONAL,
                                      UNIT_CONV_NO_UNITS, parts,
                                      rt_data.tchem.mass_fraction_HII);
  list[count++] = io_make_input_field("MassFractionHeI", FLOAT, 1, OPTIONAL,
                                      UNIT_CONV_NO_UNITS, parts,
                                      rt_data.tchem.mass_fraction_HeI);
  list[count++] = io_make_input_field("MassFractionHeII", FLOAT, 1, OPTIONAL,
                                      UNIT_CONV_NO_UNITS, parts,
                                      rt_data.tchem.mass_fraction_HeII);
  list[count++] = io_make_input_field("MassFractionHeIII", FLOAT, 1, OPTIONAL,
                                      UNIT_CONV_NO_UNITS, parts,
                                      rt_data.tchem.mass_fraction_HeIII);

  return count;
}

/**
 * @brief Specifies which particle fields to read from a dataset
 *
 * @param sparts The star particle array.
 * @param list The list of i/o properties to read.
 *
 * @return Returns the number of fields to read.
 */
INLINE static int rt_read_stars(const struct spart *sparts,
                                struct io_props *list) {
  return 0;
}

/**
 * @brief Extract radiation energies of radiation struct for all photon groups
 * Note: "allocation" of `float* ret` happens in io_copy_temp_buffer()
 *
 * @param engine the engine
 * @param part the particle to extract data from
 * @param xpart the according xpart to extract data from
 * @param ret (return) the extracted data
 */
INLINE static void rt_convert_radiation_energies(const struct engine *engine,
                                                 const struct part *part,
                                                 const struct xpart *xpart,
                                                 float *ret) {

  for (int g = 0; g < RT_NGROUPS; g++) {
    ret[g] = part->rt_data.radiation[g].energy_density * part->geometry.volume;
  }
}

/**
 * @brief Extract radiation fluxes of radiation struct for all photon groups
 * Note: "allocation" of `float* ret` happens in io_copy_temp_buffer()
 *
 * @param engine the engine
 * @param part the particle to extract data from
 * @param xpart the according xpart to extract data from
 * @param ret (return) the extracted data
 */
INLINE static void rt_convert_radiation_fluxes(const struct engine *engine,
                                               const struct part *part,
                                               const struct xpart *xpart,
                                               float *ret) {

  int i = 0;
  for (int g = 0; g < RT_NGROUPS; g++) {
    ret[i++] = part->rt_data.radiation[g].flux[0];
    ret[i++] = part->rt_data.radiation[g].flux[1];
    ret[i++] = part->rt_data.radiation[g].flux[2];
  }
}

/**
 * @brief Extract mass fractions of constituent species from tchem struct.
 * Note: "allocation" of `float* ret` happens in io_copy_temp_buffer()
 *
 * @param engine the engine
 * @param part the particle to extract data from
 * @param xpart the according xpart to extract data from
 * @param ret (return) the extracted data
 */
INLINE static void rt_convert_mass_fractions(const struct engine *engine,
                                             const struct part *part,
                                             const struct xpart *xpart,
                                             float *ret) {

  ret[0] = part->rt_data.tchem.mass_fraction_HI;
  ret[1] = part->rt_data.tchem.mass_fraction_HII;
  ret[2] = part->rt_data.tchem.mass_fraction_HeI;
  ret[3] = part->rt_data.tchem.mass_fraction_HeII;
  ret[4] = part->rt_data.tchem.mass_fraction_HeIII;
}

/**
 * @brief Creates additional output fields for the radiative
 * transfer data of hydro particles.
 *
 * @param parts The particle array.
 * @param list The list of i/o properties to write.
 *
 * @return Returns the number of fields to write.
 */
INLINE static int rt_write_particles(const struct part *parts,
                                     struct io_props *list) {

  int num_elements = 3;

  list[0] = io_make_physical_output_field_convert_part(
      "PhotonEnergies", FLOAT, RT_NGROUPS, UNIT_CONV_ENERGY, 0, parts,
      /*xparts=*/NULL,
      /*convertible to comoving=*/1, rt_convert_radiation_energies,
      "Photon Energies (all groups)");

  list[1] = io_make_physical_output_field_convert_part(
      "PhotonFluxes", FLOAT, 3 * RT_NGROUPS, UNIT_CONV_RADIATION_FLUX, 0, parts,
      /*xparts=*/NULL,
      /*convertible to comoving=*/1, rt_convert_radiation_fluxes,
      "Photon Fluxes (all groups; x, y, and z coordinates)");
  list[2] = io_make_output_field_convert_part(
      "IonMassFractions", FLOAT, 5, UNIT_CONV_NO_UNITS, 0, parts,
      /*xparts=*/NULL, rt_convert_mass_fractions,
      "Mass fractions of all constituent species");

#ifdef SWIFT_RT_DEBUG_CHECKS
  num_elements += 8;
  list[3] =
      io_make_output_field("RTDebugInjectionDone", INT, 1, UNIT_CONV_NO_UNITS,
                           0, parts, rt_data.debug_injection_done,
                           "How many times rt_injection_update_photon_density "
                           "has been called");
  list[4] = io_make_output_field(
      "RTDebugCallsIactGradientInteractions", INT, 1, UNIT_CONV_NO_UNITS, 0,
      parts, rt_data.debug_calls_iact_gradient_interaction,
      "number of calls to this particle during the gradient interaction loop "
      "from the actual interaction function");
  list[5] = io_make_output_field("RTDebugCallsIactTransportInteractions", INT,
                                 1, UNIT_CONV_NO_UNITS, 0, parts,
                                 rt_data.debug_calls_iact_transport_interaction,
                                 "number of calls to this particle during the "
                                 "transport interaction loop from the actual "
                                 "interaction function");
  list[6] =
      io_make_output_field("RTDebugGradientsDone", INT, 1, UNIT_CONV_NO_UNITS,
                           0, parts, rt_data.debug_gradients_done,
                           "How many times finalise_gradients was called");
  list[7] =
      io_make_output_field("RTDebugTransportDone", INT, 1, UNIT_CONV_NO_UNITS,
                           0, parts, rt_data.debug_transport_done,
                           "How many times finalise_transport was called");
  list[8] = io_make_output_field(
      "RTDebugThermochemistryDone", INT, 1, UNIT_CONV_NO_UNITS, 0, parts,
      rt_data.debug_thermochem_done, "How many times rt_tchem was called");
  list[9] = io_make_output_field(
      "RTDebugRadAbsorbedTot", ULONGLONG, 1, UNIT_CONV_NO_UNITS, 0, parts,
      rt_data.debug_radiation_absorbed_tot,
      "Radiation absorbed by this part during its lifetime");
  list[10] = io_make_output_field(
      "RTDebugSubcycles", INT, 1, UNIT_CONV_NO_UNITS, 0, parts,
      rt_data.debug_nsubcycles, "How many times this part was subcycled");
#endif

  return num_elements;
}

/**
 * @brief Creates additional output fields for the radiative
 * transfer data of star particles.
 *
 * @param sparts The star particle array.
 * @param list The list of i/o properties to write.
 *
 * @return Returns the number of fields to write.
 */
INLINE static int rt_write_stars(const struct spart *sparts,
                                 struct io_props *list) {
  int num_elements = 0;

#ifdef SWIFT_RT_DEBUG_CHECKS
  num_elements += 4;
  list[0] = io_make_output_field("RTDebugHydroIact", INT, 1, UNIT_CONV_NO_UNITS,
                                 0, sparts, rt_data.debug_iact_hydro_inject,
                                 "number of interactions between this star "
                                 "particle and any particle during injection");
  list[1] = io_make_output_field(
      "RTDebugEmissionRateSet", INT, 1, UNIT_CONV_NO_UNITS, 0, sparts,
      rt_data.debug_emission_rate_set, "Stellar photon emission rates set?");
  list[2] = io_make_output_field(
      "RTDebugRadEmittedTot", ULONGLONG, 1, UNIT_CONV_NO_UNITS, 0, sparts,
      rt_data.debug_radiation_emitted_tot,
      "Total radiation emitted during the lifetime of this star");
  list[3] = io_make_output_field("RTDebugInjectedPhotonEnergy", FLOAT,
                                 RT_NGROUPS, UNIT_CONV_ENERGY, 0, sparts,
                                 rt_data.debug_injected_energy_tot,
                                 "Total radiation actually injected into gas");
#endif

  return num_elements;
}

/**
 * @brief Write the RT model properties to the snapshot.
 *
 * @param h_grp The HDF5 group in which to write
 * @param h_grp_columns The HDF5 group containing named columns
 * @param e The engine
 * @param internal_units The internal unit system
 * @param snapshot_units Units used for the snapshot
 * @param rtp The #rt_props
 */
INLINE static void rt_write_flavour(hid_t h_grp, hid_t h_grp_columns,
                                    const struct engine *e,
                                    const struct unit_system *internal_units,
                                    const struct unit_system *snapshot_units,
                                    const struct rt_props *rtp) {

#if defined(HAVE_HDF5)

  /* Write scheme name */
  /* ----------------- */
  io_write_attribute_s(h_grp, "RT Scheme", RT_IMPLEMENTATION);
  io_write_attribute_s(h_grp, "RT Riemann Solver", RT_RIEMANN_SOLVER_NAME);

  /* Write photon group counts */
  /* ------------------------- */
  io_write_attribute_i(h_grp, "PhotonGroupNumber", RT_NGROUPS);

  /* Write photon group bin edges */
  /* ---------------------------- */

  /* Note: photon frequency bin edges are kept in cgs. Convert them here to
   * internal units so we're still compatible with swiftsimio. */
  const float Hz_internal =
      units_cgs_conversion_factor(internal_units, UNIT_CONV_INV_TIME);
  const float Hz_internal_inv = 1.f / Hz_internal;
  float photon_groups_internal[RT_NGROUPS];
  for (int g = 0; g < RT_NGROUPS; g++)
    photon_groups_internal[g] = rtp->photon_groups[g] * Hz_internal_inv;

  hid_t type_float = H5Tcopy(io_hdf5_type(FLOAT));

  hsize_t dims[1] = {RT_NGROUPS};
  hid_t space = H5Screate_simple(1, dims, NULL);
  hid_t dset = H5Dcreate(h_grp, "PhotonGroupEdges", type_float, space,
                         H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
  H5Dwrite(dset, type_float, H5S_ALL, H5S_ALL, H5P_DEFAULT,
           photon_groups_internal);

  /* Write unit conversion factors for this data set */
  char buffer[FIELD_BUFFER_SIZE] = {0};
  units_cgs_conversion_string(buffer, snapshot_units, UNIT_CONV_INV_TIME,
                              /*scale_factor_exponent=*/0);
  float baseUnitsExp[5];
  units_get_base_unit_exponents_array(baseUnitsExp, UNIT_CONV_INV_TIME);
  io_write_attribute_f(dset, "U_M exponent", baseUnitsExp[UNIT_MASS]);
  io_write_attribute_f(dset, "U_L exponent", baseUnitsExp[UNIT_LENGTH]);
  io_write_attribute_f(dset, "U_t exponent", baseUnitsExp[UNIT_TIME]);
  io_write_attribute_f(dset, "U_I exponent", baseUnitsExp[UNIT_CURRENT]);
  io_write_attribute_f(dset, "U_T exponent", baseUnitsExp[UNIT_TEMPERATURE]);
  io_write_attribute_f(dset, "h-scale exponent", 0.f);
  io_write_attribute_f(dset, "a-scale exponent", 0.f);
  io_write_attribute_s(dset, "Expression for physical CGS units", buffer);
  io_write_attribute_b(dset, "Value stored as physical", 1);
  io_write_attribute_b(dset, "Property can be converted to comoving", 0);

  /* Write the actual number this conversion factor corresponds to */
  const double factor =
      units_cgs_conversion_factor(snapshot_units, UNIT_CONV_INV_TIME);
  io_write_attribute_d(
      dset, "Conversion factor to CGS (not including cosmological corrections)",
      factor);
  io_write_attribute_d(
      dset,
      "Conversion factor to physical CGS (including cosmological corrections)",
      factor * pow(e->cosmology->a, 0.f));

  H5Dclose(dset);
  /* H5Tclose(type_float); [> close this later <] */

  /* If without RT, we have nothing more to do. */
  const int with_rt = e->policy & engine_policy_rt;
  if (!with_rt) return;

  /* Write photon group names */
  /* -------------------------*/

  /* Generate Energy Group names */
  char names_energy[RT_NGROUPS * RT_LABELS_SIZE];
  for (int g = 0; g < RT_NGROUPS; g++) {
    char newEname[RT_LABELS_SIZE];
    sprintf(newEname, "Group%d", g + 1);
    strcpy(names_energy + g * RT_LABELS_SIZE, newEname);
  }

  /* Now write them down */
  hid_t type_string_label = H5Tcopy(H5T_C_S1);
  H5Tset_size(type_string_label, RT_LABELS_SIZE);

  hsize_t dimsE[1] = {RT_NGROUPS};
  hid_t spaceE = H5Screate_simple(1, dimsE, NULL);
  hid_t dsetE = H5Dcreate(h_grp_columns, "PhotonEnergies", type_string_label,
                          spaceE, H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
  H5Dwrite(dsetE, type_string_label, H5S_ALL, H5S_ALL, H5P_DEFAULT,
           names_energy);
  H5Dclose(dsetE);

  /* Generate Fluxes Group Names */
  char names_fluxes[3 * RT_NGROUPS * RT_LABELS_SIZE];
  int i = 0;
  for (int g = 0; g < RT_NGROUPS; g++) {
    char newFnameX[RT_LABELS_SIZE];
    sprintf(newFnameX, "Group%dX", g + 1);
    strcpy(names_fluxes + i * RT_LABELS_SIZE, newFnameX);
    i++;
    char newFnameY[RT_LABELS_SIZE];
    sprintf(newFnameY, "Group%dY", g + 1);
    strcpy(names_fluxes + i * RT_LABELS_SIZE, newFnameY);
    i++;
    char newFnameZ[RT_LABELS_SIZE];
    sprintf(newFnameZ, "Group%dZ", g + 1);
    strcpy(names_fluxes + i * RT_LABELS_SIZE, newFnameZ);
    i++;
  }

  /* Now write them down */
  hsize_t dimsF[1] = {3 * RT_NGROUPS};
  hid_t spaceF = H5Screate_simple(1, dimsF, NULL);
  hid_t dsetF = H5Dcreate(h_grp_columns, "PhotonFluxes", type_string_label,
                          spaceF, H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
  H5Dwrite(dsetF, type_string_label, H5S_ALL, H5S_ALL, H5P_DEFAULT,
           names_fluxes);
  H5Dclose(dsetF);

  /* H5Tclose(type_string_label); [> close this later <] */

  /* Write reduced speed of light */
  /* ---------------------------- */
  /* hid_t type2 = H5Tcopy(io_hdf5_type(FLOAT)); */

  hsize_t dims_cred[1] = {1};
  hid_t space_cred = H5Screate_simple(1, dims_cred, NULL);
  hid_t dset_cred =
      H5Dcreate(h_grp, "ReducedLightspeed", type_float, space_cred, H5P_DEFAULT,
                H5P_DEFAULT, H5P_DEFAULT);
  H5Dwrite(dset_cred, type_float, H5S_ALL, H5S_ALL, H5P_DEFAULT,
           &rt_params.reduced_speed_of_light);

  /* Write unit conversion factors for this data set */
  char buffer_cred[FIELD_BUFFER_SIZE] = {0};
  units_cgs_conversion_string(buffer_cred, snapshot_units, UNIT_CONV_VELOCITY,
                              /*scale_factor_exponent=*/0);
  float baseUnitsExp_cred[5];
  units_get_base_unit_exponents_array(baseUnitsExp_cred, UNIT_CONV_VELOCITY);
  io_write_attribute_f(dset_cred, "U_M exponent", baseUnitsExp_cred[UNIT_MASS]);
  io_write_attribute_f(dset_cred, "U_L exponent",
                       baseUnitsExp_cred[UNIT_LENGTH]);
  io_write_attribute_f(dset_cred, "U_t exponent", baseUnitsExp_cred[UNIT_TIME]);
  io_write_attribute_f(dset_cred, "U_I exponent",
                       baseUnitsExp_cred[UNIT_CURRENT]);
  io_write_attribute_f(dset_cred, "U_T exponent",
                       baseUnitsExp_cred[UNIT_TEMPERATURE]);
  io_write_attribute_f(dset_cred, "h-scale exponent", 0.f);
  io_write_attribute_f(dset_cred, "a-scale exponent", 0.f);
  io_write_attribute_s(dset_cred, "Expression for physical CGS units",
                       buffer_cred);
  io_write_attribute_b(dset_cred, "Value stored as physical", 1);
  io_write_attribute_b(dset_cred, "Property can be converted to comoving", 0);

  /* Write the actual number this conversion factor corresponds to */
  /* TODO Mladen: check cosmology. reduced_speed_of_light is physical only for
   * now. */
  const double factor_cred =
      units_cgs_conversion_factor(snapshot_units, UNIT_CONV_VELOCITY);
  io_write_attribute_d(
      dset_cred,
      "Conversion factor to CGS (not including cosmological corrections)",
      factor_cred);
  io_write_attribute_d(
      dset_cred,
      "Conversion factor to physical CGS (including cosmological corrections)",
      factor_cred * pow(e->cosmology->a, 0.f));

  H5Dclose(dset_cred);

  /* Write constituent species mass fractions */
  /* ---------------------------------------- */

  char names_mf[5 * RT_LABELS_SIZE];
  strcpy(names_mf + 0 * RT_LABELS_SIZE, "HI\0");
  strcpy(names_mf + 1 * RT_LABELS_SIZE, "HII\0");
  strcpy(names_mf + 2 * RT_LABELS_SIZE, "HeI\0");
  strcpy(names_mf + 3 * RT_LABELS_SIZE, "HeII\0");
  strcpy(names_mf + 4 * RT_LABELS_SIZE, "HeIII\0");

  hsize_t dims_mf[1] = {5};
  hid_t space_mf = H5Screate_simple(1, dims_mf, NULL);
  hid_t dset_mf =
      H5Dcreate(h_grp_columns, "IonMassFractions", type_string_label, space_mf,
                H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
  H5Dwrite(dset_mf, type_string_label, H5S_ALL, H5S_ALL, H5P_DEFAULT, names_mf);
  H5Dclose(dset_mf);

  /* Clean up after yourself */
  /* ----------------------- */

  /* Close up the types */
  H5Tclose(type_float);
  H5Tclose(type_string_label);

#endif /* HAVE_HDF5 */
}

#endif /* SWIFT_RT_IO_GEAR_H */
