/**
 * @file Backup file containing the now-defunct Newton integration scheme
 * for the cooling.
 */

/**
 * @brief Computes the log_10 of the temperature corresponding to a given
 * internal energy, hydrogen number density, Helium fraction and redshift.
 *
 * Note that the redshift is implicitly passed in via the currently loaded
 * tables in the #cooling_function_data.
 *
 * For the low-z case, we interpolate the flattened 4D table 'u_to_temp' that
 * is arranged in the following way:
 * - 1st dim: redshift, length = eagle_cooling_N_loaded_redshifts
 * - 2nd dim: Hydrogen density, length = eagle_cooling_N_density
 * - 3rd dim: Helium fraction, length = eagle_cooling_N_He_frac
 * - 4th dim: Internal energy, length = eagle_cooling_N_temperature
 *
 * For the high-z case, we interpolate the flattened 3D table 'u_to_temp' that
 * is arranged in the following way:
 * - 1st dim: Hydrogen density, length = eagle_cooling_N_density
 * - 2nd dim: Helium fraction, length = eagle_cooling_N_He_frac
 * - 3rd dim: Internal energy, length = eagle_cooling_N_temperature
 *
 * @param log_10_u_cgs Log base 10 of internal energy in cgs.
 * @param redshift Current redshift.
 * @param n_H_index Index along the Hydrogen density dimension.
 * @param He_index Index along the Helium fraction dimension.
 * @param d_n_H Offset between Hydrogen density and table[n_H_index].
 * @param d_He Offset between helium fraction and table[He_index].
 * @param cooling #cooling_function_data structure.
 *
 * @param compute_dT_du Do we want to compute dT/du ?
 * @param dT_du (return) The value of dT/du
 *
 * @return log_10 of the temperature.
 */
__attribute__((always_inline)) INLINE double eagle_convert_u_to_temp(
    const double log_10_u_cgs, const float redshift, const int compute_dT_du,
    float *dT_du, int n_H_index, int He_index, float d_n_H, float d_He,
    const struct cooling_function_data *restrict cooling) {

  /* Get index of u along the internal energy axis */
  int u_index;
  float d_u;
  get_index_1d(cooling->Therm, eagle_cooling_N_temperature, log_10_u_cgs,
               &u_index, &d_u);

  /* Interpolate temperature table to return temperature for current
   * internal energy (use 3D interpolation for high redshift table,
   * otherwise 4D) */
  float log_10_T;
  if (redshift > cooling->Redshifts[eagle_cooling_N_redshifts - 1]) {

    log_10_T = interpolation_3d(cooling->table.temperature,   /* */
                                n_H_index, He_index, u_index, /* */
                                d_n_H, d_He, d_u,             /* */
                                eagle_cooling_N_density,      /* */
                                eagle_cooling_N_He_frac,      /* */
                                eagle_cooling_N_temperature); /* */
  } else {

    log_10_T =
        interpolation_4d(cooling->table.temperature,                  /* */
                         /*z_index=*/0, n_H_index, He_index, u_index, /* */
                         cooling->dz, d_n_H, d_He, d_u,               /* */
                         eagle_cooling_N_loaded_redshifts,            /* */
                         eagle_cooling_N_density,                     /* */
                         eagle_cooling_N_He_frac,                     /* */
                         eagle_cooling_N_temperature);                /* */
  }

  if (compute_dT_du) {

    float log_10_T_high, log_10_T_low;

    /* Interpolate temperature table to return temperature for internal energy
     * at grid point above current internal energy for computing dT_du used for
     * calculation of dlambda_du in cooling.c (use 3D interpolation for high
     * redshift table, otherwise 4D) */
    if (redshift > cooling->Redshifts[eagle_cooling_N_redshifts - 1]) {

      log_10_T_high = interpolation_3d(cooling->table.temperature,   /* */
                                       n_H_index, He_index, u_index, /* */
                                       d_n_H, d_He, /*delta_u=*/1.f, /* */
                                       eagle_cooling_N_density,      /* */
                                       eagle_cooling_N_He_frac,      /* */
                                       eagle_cooling_N_temperature); /* */

    } else {

      log_10_T_high =
          interpolation_4d(cooling->table.temperature,                  /* */
                           /*z_index=*/0, n_H_index, He_index, u_index, /* */
                           cooling->dz, d_n_H, d_He, /*delta_u=*/1.f,   /* */
                           eagle_cooling_N_loaded_redshifts,            /* */
                           eagle_cooling_N_density,                     /* */
                           eagle_cooling_N_He_frac,                     /* */
                           eagle_cooling_N_temperature);                /* */
    }

    /* Interpolate temperature table to return temperature for internal energy
     * at grid point below current internal energy for computing dT_du used for
     * calculation of dlambda_du in cooling.c (use 3D interpolation for high
     * redshift table, otherwise 4D) */
    if (redshift > cooling->Redshifts[eagle_cooling_N_redshifts - 1]) {

      log_10_T_low = interpolation_3d(cooling->table.temperature,   /* */
                                      n_H_index, He_index, u_index, /* */
                                      d_n_H, d_He, /*delta_u=*/0.f, /* */
                                      eagle_cooling_N_density,      /* */
                                      eagle_cooling_N_He_frac,      /* */
                                      eagle_cooling_N_temperature); /* */

    } else {

      log_10_T_low =
          interpolation_4d(cooling->table.temperature,                  /* */
                           /*z_index=*/0, n_H_index, He_index, u_index, /* */
                           cooling->dz, d_n_H, d_He, /*delta_u=*/0.f,   /* */
                           eagle_cooling_N_loaded_redshifts,            /* */
                           eagle_cooling_N_density,                     /* */
                           eagle_cooling_N_He_frac,                     /* */
                           eagle_cooling_N_temperature);                /* */
    }

    /* Calculate dT/du */
    const float delta_u = exp(cooling->Therm[u_index + 1] * M_LN10) -
                          exp(cooling->Therm[u_index] * M_LN10);
    *dT_du =
        (exp(M_LN10 * log_10_T_high) - exp(M_LN10 * log_10_T_low)) / delta_u;
  }

  /* Special case for temperatures below the start of the table */
  if (u_index == 0 && d_u == 0.f) {

    /* The temperature is multiplied by u / 10^T[0]
     * where T[0] is the first entry in the table */
    log_10_T += log_10_u_cgs - cooling->Temp[0];
  }

  return log_10_T;
}

/**
 * @brief Computes the cooling rate corresponding to a given internal energy,
 * hydrogen number density, Helium fraction, redshift and metallicity from
 * all the possible channels.
 *
 * 1) Metal-free cooling:
 * We interpolate the flattened 4D table 'H_and_He_net_heating' that is
 * arranged in the following way:
 * - 1st dim: redshift, length = eagle_cooling_N_loaded_redshifts
 * - 2nd dim: Hydrogen density, length = eagle_cooling_N_density
 * - 3rd dim: Helium fraction, length = eagle_cooling_N_He_frac
 * - 4th dim: Internal energy, length = eagle_cooling_N_temperature
 *
 * 2) Electron abundance
 * We compute the electron abundance by interpolating the flattened 4d table
 * 'H_and_He_electron_abundance' that is arranged in the following way:
 * - 1st dim: redshift, length = eagle_cooling_N_loaded_redshifts
 * - 2nd dim: Hydrogen density, length = eagle_cooling_N_density
 * - 3rd dim: Helium fraction, length = eagle_cooling_N_He_frac
 * - 4th dim: Internal energy, length = eagle_cooling_N_temperature
 *
 * 3) Compton cooling is applied via the analytic formula.
 *
 * 4) Solar electron abudance
 * We compute the solar electron abundance by interpolating the flattened 3d
 * table 'solar_electron_abundance' that is arranged in the following way:
 * - 1st dim: redshift, length = eagle_cooling_N_loaded_redshifts
 * - 2nd dim: Hydrogen density, length = eagle_cooling_N_density
 * - 3rd dim: Internal energy, length = eagle_cooling_N_temperature
 *
 * 5) Metal-line cooling
 * For each tracked element we interpolate the flattened 4D table
 * 'table_metals_net_heating' that is arrange in the following way:
 * - 1st dim: element, length = eagle_cooling_N_metal
 * - 2nd dim: redshift, length = eagle_cooling_N_loaded_redshifts
 * - 3rd dim: Hydrogen density, length = eagle_cooling_N_density
 * - 4th dim: Internal energy, length = eagle_cooling_N_temperature
 *
 * Note that this is a fake 4D interpolation as we do not interpolate
 * along the 1st dimension. We just do this once per element.
 *
 * Since only the temperature changes when cooling a given particle,
 * the redshift, hydrogen number density and helium fraction indices
 * and offsets passed in.
 *
 * If the arguement dlambda_du is non-NULL, the routine also
 * calculates derivative of cooling rate with respect to internal
 * energy.
 *
 * If the argument element_lambda is non-NULL, the routine also
 * returns the cooling rate per element in the array.
 *
 * @param log10_u_cgs Log base 10 of internal energy per unit mass in CGS units.
 * @param redshift The current redshift
 * @param n_H_cgs The Hydrogen number density in CGS units.
 * @param solar_ratio Array of ratios of particle metal abundances
 * to solar metal abundances
 *
 * @param n_H_index Particle hydrogen number density index
 * @param d_n_H Particle hydrogen number density offset
 * @param He_index Particle helium fraction index
 * @param d_He Particle helium fraction offset
 * @param cooling Cooling data structure
 *
 * @param dlambda_du (return) Derivative of the cooling rate with respect to u.
 * @param element_lambda (return) Cooling rate from each element
 *
 * @return The cooling rate
 */
INLINE static double eagle_metal_cooling_rate(
    double log10_u_cgs, double redshift, double n_H_cgs,
    const float solar_ratio[eagle_cooling_N_abundances], int n_H_index,
    float d_n_H, int He_index, float d_He,
    const struct cooling_function_data *restrict cooling, double *dlambda_du,
    double *element_lambda) {

  /* used for calculating dlambda_du */
  double temp_lambda_high = 0, temp_lambda_low = 0;
  double h_plus_he_electron_abundance_high = 0;
  double h_plus_he_electron_abundance_low = 0;
  double solar_electron_abundance_high = 0;
  double solar_electron_abundance_low = 0;
  double elem_cool_low = 0, elem_cool_high = 0;

  /* We only need dT_du if dLambda_du is non-NULL */
  const int compute_dT_du = (dlambda_du != NULL) ? 1 : 0;

  /* Temperature */
  float dT_du = -1.f;
  const double log_10_T =
      eagle_convert_u_to_temp(log10_u_cgs, redshift, compute_dT_du, &dT_du,
                              n_H_index, He_index, d_n_H, d_He, cooling);

  /* Get index along temperature dimension of the tables */
  int T_index;
  float d_T;
  get_index_1d(cooling->Temp, eagle_cooling_N_temperature, log_10_T, &T_index,
               &d_T);

  /* Difference between entries on the temperature table around u */
  const float delta_T = exp(M_LN10 * cooling->Temp[T_index + 1]) -
                        exp(M_LN10 * cooling->Temp[T_index]);

  /**********************/
  /* Metal-free cooling */
  /**********************/

  double Lambda_free;

  if (redshift > cooling->Redshifts[eagle_cooling_N_redshifts - 1]) {

    /* If we're using the high redshift tables then we don't interpolate
     * in redshift */
    Lambda_free = interpolation_3d(cooling->table.H_plus_He_heating, /* */
                                   n_H_index, He_index, T_index,     /* */
                                   d_n_H, d_He, d_T,                 /* */
                                   eagle_cooling_N_density,          /* */
                                   eagle_cooling_N_He_frac,          /* */
                                   eagle_cooling_N_temperature);     /* */

    /* compute values at temperature gridpoints above and below input
     * temperature for calculation of dlambda_du. Pass in NULL pointer for
     * dlambda_du in order to skip */
    if (dlambda_du != NULL) {
      temp_lambda_high = interpolation_3d(
          cooling->table.H_plus_He_heating, n_H_index, He_index, T_index, d_n_h,
          d_He, 1.f, cooling->N_nH, cooling->N_He, cooling->N_Temp);
      temp_lambda_low = interpolation_3d(
          cooling->table.H_plus_He_heating, n_H_index, He_index, T_index, d_n_h,
          d_He, 0.f, cooling->N_nH, cooling->N_He, cooling->N_Temp);
    }

  } else {

    /* Using normal tables, have to interpolate in redshift */
    Lambda_free =
        interpolation_4d(cooling->table.H_plus_He_heating,            /* */
                         /*z_index=*/0, n_H_index, He_index, T_index, /* */
                         cooling->dz, d_n_H, d_He, d_T,               /* */
                         eagle_cooling_N_loaded_redshifts,            /* */
                         eagle_cooling_N_density,                     /* */
                         eagle_cooling_N_He_frac,                     /* */
                         eagle_cooling_N_temperature);                /* */

    /* compute values at temperature gridpoints above and below input
     * temperature for calculation of dlambda_du */
    if (dlambda_du != NULL) {
      temp_lambda_high =
          interpolation_4d(cooling->table.H_plus_He_heating, 0, n_H_index,
                           He_index, T_index, cooling->dz, d_n_h, d_He, 1.f, 2,
                           cooling->N_nH, cooling->N_He, cooling->N_Temp);
      temp_lambda_low =
          interpolation_4d(cooling->table.H_plus_He_heating, 0, n_H_index,
                           He_index, T_index, cooling->dz, d_n_h, d_He, 0.f, 2,
                           cooling->N_nH, cooling->N_He, cooling->N_Temp);
    }
  }

  if (dlambda_du != NULL) {
    *dlambda_du += (temp_lambda_high - temp_lambda_low) / delta_T * dT_du;
  }

  /* If we're testing cooling rate contributions write to array */
  if (element_lambda != NULL) {
    element_lambda[0] = Lambda_free;
  }

  /**********************/
  /* Electron abundance */
  /**********************/

  double H_plus_He_electron_abundance;

  if (redshift > cooling->Redshifts[eagle_cooling_N_redshifts - 1]) {

    H_plus_He_electron_abundance =
        interpolation_3d(cooling->table.H_plus_He_electron_abundance, /* */
                         n_H_index, He_index, T_index,                /* */
                         d_n_H, d_He, d_T,                            /* */
                         eagle_cooling_N_density,                     /* */
                         eagle_cooling_N_He_frac,                     /* */
                         eagle_cooling_N_temperature);                /* */
    /* compute values at temperature gridpoints above and below input
     * temperature for calculation of dlambda_du. Pass in NULL pointer for
     * dlambda_du in order to skip */

    h_plus_he_electron_abundance_high =
        interpolation_3d(cooling->table.H_plus_He_electron_abundance, n_H_index,
                         He_index, T_index, d_n_h, d_He, 1.f, cooling->N_nH,
                         cooling->N_He, cooling->N_Temp);
    h_plus_he_electron_abundance_low =
        interpolation_3d(cooling->table.H_plus_He_electron_abundance, n_H_index,
                         He_index, T_index, d_n_h, d_He, 0.f, cooling->N_nH,
                         cooling->N_He, cooling->N_Temp);

  } else {

    H_plus_He_electron_abundance =
        interpolation_4d(cooling->table.H_plus_He_electron_abundance, /* */
                         /*z_index=*/0, n_H_index, He_index, T_index, /* */
                         cooling->dz, d_n_H, d_He, d_T,               /* */
                         eagle_cooling_N_loaded_redshifts,            /* */
                         eagle_cooling_N_density,                     /* */
                         eagle_cooling_N_He_frac,                     /* */
                         eagle_cooling_N_temperature);                /* */

    /* compute values at temperature gridpoints above and below input
     * temperature for calculation of dlambda_du */
    h_plus_he_electron_abundance_high =
        interpolation_4d(cooling->table.H_plus_He_electron_abundance, 0,
                         n_H_index, He_index, T_index, cooling->dz, d_n_h, d_He,
                         1.f, 2, cooling->N_nH, cooling->N_He, cooling->N_Temp);
    h_plus_he_electron_abundance_low =
        interpolation_4d(cooling->table.H_plus_He_electron_abundance, 0,
                         n_H_index, He_index, T_index, cooling->dz, d_n_h, d_He,
                         0.f, 2, cooling->N_nH, cooling->N_He, cooling->N_Temp);
  }

  /**********************/
  /* Compton cooling    */
  /**********************/

  double Lambda_Compton = 0.;

  /* Do we need to add the inverse Compton cooling? */
  /* It is *not* stored in the tables before re-ionisation */
  if ((redshift > cooling->Redshifts[eagle_cooling_N_redshifts - 1]) ||
      (redshift > cooling->H_reion_z)) {

    const double T = exp10(log_10_T);

    /* Note the minus sign */
    Lambda_Compton -= eagle_Compton_cooling_rate(cooling, redshift, n_H_cgs, T,
                                                 H_plus_He_electron_abundance);
  }

  /* If we're testing cooling rate contributions write to array */
  if (element_lambda != NULL) {
    element_lambda[1] = Lambda_Compton;
  }

  /*******************************/
  /* Solar electron abundance    */
  /*******************************/

  double solar_electron_abundance;

  if (redshift > cooling->Redshifts[eagle_cooling_N_redshifts - 1]) {

    /* If we're using the high redshift tables then we don't interpolate
     * in redshift */
    solar_electron_abundance =
        interpolation_2d(cooling->table.electron_abundance, /* */
                         n_H_index, T_index,                /* */
                         d_n_H, d_T,                        /* */
                         eagle_cooling_N_density,           /* */
                         eagle_cooling_N_temperature);      /* */

    /* compute values at temperature gridpoints above and below input
     * temperature for calculation of dlambda_du */
    if (dlambda_du != NULL) {
      solar_electron_abundance_high =
          interpolation_2d(cooling->table.electron_abundance, n_H_index,
                           T_index, d_n_h, 1.f, cooling->N_nH, cooling->N_Temp);
      solar_electron_abundance_low =
          interpolation_2d(cooling->table.electron_abundance, n_H_index,
                           T_index, d_n_h, 0.f, cooling->N_nH, cooling->N_Temp);
    }

  } else {

    /* Using normal tables, have to interpolate in redshift */
    solar_electron_abundance =
        interpolation_3d(cooling->table.electron_abundance, /* */
                         /*z_index=*/0, n_H_index, T_index, /* */
                         cooling->dz, d_n_H, d_T,           /* */
                         eagle_cooling_N_loaded_redshifts,  /* */
                         eagle_cooling_N_density,           /* */
                         eagle_cooling_N_temperature);      /* */

    /* compute values at temperature gridpoints above and below input
     * temperature for calculation of dlambda_du */
    if (dlambda_du != NULL) {
      solar_electron_abundance_high = interpolation_3d(
          cooling->table.electron_abundance, 0, n_H_index, T_index, cooling->dz,
          d_n_h, 1.f, 2, cooling->N_nH, cooling->N_Temp);
      solar_electron_abundance_low = interpolation_3d(
          cooling->table.electron_abundance, 0, n_H_index, T_index, cooling->dz,
          d_n_h, 0.f, 2, cooling->N_nH, cooling->N_Temp);
    }

    const double electron_abundance_ratio =
        H_plus_He_electron_abundance / solar_electron_abundance;

    /**********************/
    /* Metal-line cooling */
    /**********************/

    /* for each element the cooling rate is multiplied by the ratio of H, He
     * electron abundance to solar electron abundance then by the ratio of the
     * particle metal abundance to solar metal abundance. */

    double lambda_metal[eagle_cooling_N_metal + 2] = {0.};

    if (redshift > cooling->Redshifts[eagle_cooling_N_redshifts - 1]) {

      /* Loop over the metals (ignore H and He) */
      for (int elem = 2; elem < eagle_cooling_N_metal + 2; elem++) {

        if (solar_ratio[elem] > 0.) {

          /* Note that we do not interpolate along the x-axis
           * (element dimension) */
          lambda_metal[elem] =
              interpolation_3d_no_x(cooling->table.metal_heating,   /* */
                                    elem - 2, n_H_index, T_index,   /* */
                                    /*delta_elem=*/0.f, d_n_H, d_T, /* */
                                    eagle_cooling_N_metal,          /* */
                                    eagle_cooling_N_density,        /* */
                                    eagle_cooling_N_temperature);   /* */

          lambda_metal[elem] *= electron_abundance_ratio;
          lambda_metal[elem] *= solar_ratio[elem];
        }

        /* compute values at temperature gridpoints above and below input
         * temperature for calculation of dlambda_du */
        if (dlambda_du != NULL) {
          elem_cool_high = interpolation_3d_no_x(
              cooling->table.metal_heating, elem, n_H_index, T_index, 0.f,
              d_n_h, 1.f, cooling->N_Elements, cooling->N_nH, cooling->N_Temp);

          elem_cool_low = interpolation_3d_no_x(
              cooling->table.metal_heating, elem, n_H_index, T_index, 0.f,
              d_n_h, 0.f, cooling->N_nH, cooling->N_Temp, cooling->N_Elements);

          *dlambda_du += (elem_cool_high * h_plus_he_electron_abundance_high /
                              solar_electron_abundance_high -
                          elem_cool_low * h_plus_he_electron_abundance_low /
                              solar_electron_abundance_low) /
                         delta_T * dT_du * solar_ratio[elem + 2];
        }
      }

    } else {

      /* Loop over the metals (ignore H and He) */
      for (int elem = 2; elem < eagle_cooling_N_metal + 2; elem++) {

        if (solar_ratio[elem] > 0.) {

          /* Note that we do not interpolate along the x-axis
           * (element dimension) */
          lambda_metal[elem] = interpolation_4d_no_x(
              cooling->table.metal_heating,                /* */
              elem - 2, /*z_index=*/0, n_H_index, T_index, /* */
              /*delta_elem=*/0.f, cooling->dz, d_n_H, d_T, /* */
              eagle_cooling_N_metal,                       /* */
              eagle_cooling_N_loaded_redshifts,            /* */
              eagle_cooling_N_density,                     /* */
              eagle_cooling_N_temperature);                /* */

          lambda_metal[elem] *= electron_abundance_ratio;
          lambda_metal[elem] *= solar_ratio[elem];
        }

        /* compute values at temperature gridpoints above and below input
         * temperature for calculation of dlambda_du */
        if (dlambda_du != NULL) {
          elem_cool_high = interpolation_4d_no_x(
              cooling->table.metal_heating, elem, 0, n_H_index, T_index, 0.,
              cooling->dz, d_n_h, 1.f, cooling->N_Elements, 2, cooling->N_nH,
              cooling->N_Temp);

          elem_cool_low = interpolation_4d_no_x(
              cooling->table.metal_heating, elem, 0, n_H_index, T_index, 0.,
              cooling->dz, d_n_h, 0.f, cooling->N_Elements, 2, cooling->N_nH,
              cooling->N_Temp);

          *dlambda_du += (elem_cool_high * h_plus_he_electron_abundance_high /
                              solar_electron_abundance_high -
                          elem_cool_low * h_plus_he_electron_abundance_low /
                              solar_electron_abundance_low) /
                         delta_T * dT_du * solar_ratio[elem + 2];
        }
      }
    }

    if (element_lambda != NULL) {
      for (int elem = 2; elem < eagle_cooling_N_metal + 2; ++elem) {
        element_lambda[elem] = lambda_metal[elem];
      }
    }

    /* Sum up all the contributions */
    double Lambda_net = Lambda_free + Lambda_Compton;
    for (int elem = 2; elem < eagle_cooling_N_metal + 2; ++elem) {
      Lambda_net += lambda_metal[elem];
    }

    return Lambda_net;
  }

  /**
   * @brief Wrapper function used to calculate cooling rate and dLambda_du.
   * Table indices and offsets for redshift, hydrogen number density and
   * helium fraction are passed it so as to compute them only once per particle.
   *
   * @param log_u_cgs Natural log of internal energy per unit mass in CGS units.
   * @param redshift The current redshift.
   * @param n_H_cgs Hydrogen number density in CGS units.
   * @param abundance_ratio Ratio of element abundance to solar.
   *
   * @param n_H_index Particle hydrogen number density index
   * @param d_n_H Particle hydrogen number density offset
   * @param He_index Particle helium fraction index
   * @param d_He Particle helium fraction offset
   * @param cooling #cooling_function_data structure
   *
   * @param dLambdaNet_du (return) Derivative of the cooling rate with respect
   * to u.
   *
   * @return The cooling rate
   */
  INLINE static double eagle_cooling_rate(
      double log_u_cgs, double redshift, double n_H_cgs,
      const float abundance_ratio[eagle_cooling_N_abundances], int n_H_index,
      float d_n_H, int He_index, float d_He,
      const struct cooling_function_data *restrict cooling,
      double *dLambdaNet_du) {

    return eagle_metal_cooling_rate(log_u_cgs / M_LN10, redshift, n_H_cgs,
                                    abundance_ratio, n_H_index, d_n_H, He_index,
                                    d_He, cooling, dLambdaNet_du,
                                    /*element_lambda=*/NULL);
  }

  /**
   * @brief Newton Raphson integration scheme to calculate particle cooling over
   * timestep. This replaces bisection scheme used in EAGLE to minimize the
   * number of array accesses. Integration defaults to bisection scheme (see
   * function bisection_iter) if this function does not converge within a
   * specified number of steps
   *
   * @param logu_init Initial guess for log(internal energy)
   * @param u_ini Internal energy at beginning of hydro step
   * @param n_H_index Particle hydrogen number density index
   * @param d_n_H Particle hydrogen number density offset
   * @param He_index Particle helium fraction index
   * @param d_He Particle helium fraction offset
   * @param He_reion_heat Heating due to helium reionization
   * (only depends on redshift, so passed as parameter)
   * @param p #part structure
   * @param cosmo #cosmology structure
   * @param cooling #cooling_function_data structure
   * @param phys_const #phys_const data structure
   * @param abundance_ratio Array of ratios of metal abundance to solar
   * @param dt timestep
   * @param bisection_flag Flag to identify if scheme failed to converge
   */
  float newton_iter(float logu_init, double u_ini, int n_H_index, float d_n_H,
                    int He_index, float d_He, float He_reion_heat,
                    struct part *restrict p,
                    const struct cosmology *restrict cosmo,
                    const struct cooling_function_data *restrict cooling,
                    const struct phys_const *restrict phys_const,
                    const float abundance_ratio[eagle_cooling_N_abundances],
                    float dt, int *bisection_flag) {

    double logu, logu_old;
    double dLambdaNet_du = 0.0, LambdaNet;

    /* table bounds */
    const float log_table_bound_high =
        (cooling->Therm[eagle_cooling_N_temperature - 1] - 0.05) / M_LOG10E;
    const float log_table_bound_low = (cooling->Therm[0] + 0.05) / M_LOG10E;

    /* convert Hydrogen mass fraction in Hydrogen number density */
    float const *metal_fraction =
        chemistry_get_metal_mass_fraction_for_cooling(p);
    const float XH = metal_fraction[chemistry_element_H];
    const double n_H = hydro_get_physical_density(p, cosmo) * XH /
                       phys_const->const_proton_mass;
    const double n_H_cgs = n_H * cooling->number_density_to_cgs;

    /* compute ratefact = n_H * n_H / rho; Might lead to round-off error:
     * replaced by equivalent expression below */
    const double ratefact_cgs = n_H_cgs * XH * cooling->inv_proton_mass_cgs;

    logu_old = logu_init;
    logu = logu_old;
    int i = 0;

    float LambdaNet_old = 0;
    LambdaNet = 0;
    do /* iterate to convergence */
    {
      logu_old = logu;
      LambdaNet_old = LambdaNet;
      LambdaNet = (He_reion_heat / (dt * ratefact_cgs)) +
                  eagle_cooling_rate(logu_old, cosmo->z, n_H_cgs,
                                     abundance_ratio, n_H_index, d_n_H,
                                     He_index, d_He, cooling, &dLambdaNet_du);

      /* Newton iteration. For details on how the cooling equation is integrated
       * see documentation in theory/Cooling/ */
      logu = logu_old - (1.0 - u_ini * exp(-logu_old) -
                         LambdaNet * ratefact_cgs * dt * exp(-logu_old)) /
                            (1.0 - dLambdaNet_du * ratefact_cgs * dt);
      /* Check if first step passes over equilibrium solution, if it does adjust
       * next guess */
      if (i == 1 && LambdaNet_old * LambdaNet < 0)
        logu = newton_log_u_guess_cgs;

      /* check whether iterations go within about 10% of the table bounds,
       * if they do default to bisection method */
      if (logu > log_table_bound_high) {
        i = newton_max_iterations;
        break;
      } else if (logu < log_table_bound_low) {
        i = newton_max_iterations;
        break;
      }

      i++;
    } while (fabs(logu - logu_old) > newton_tolerance &&
             i < newton_max_iterations);
    if (i >= newton_max_iterations) {
      /* flag to trigger bisection scheme */
      *bisection_flag = 1;
    }

    return logu;
  }
