Back port some fixes from James and Sylvia on the calculation of the...

Back port some fixes from James and Sylvia on the calculation of the equilibirum pressure on the EoS when using cooling from the Ploeckinger tables

src/cooling/COLIBRE/cooling_subgrid.h

@@ -430,7 +430,7 @@ double compute_subgrid_property(
 
     for (int i = iden + 1; i < colibre_cooling_N_density; i++) {
 
-      float log10_Peq_current;
+      float log10_Peq_current, log10_n_at_Peq;
 
       if ((log10_Teq <= log10_T) &&
           (logZZsol <=
@@ -458,6 +458,16 @@ double compute_subgrid_property(
 
       if (log10_Peq_current >= log10_P_cgs) {
 
+        if (log10_Peq_current == log10_Peq_prev) {
+
+          /* the interpolated pressure exactly equals the input pressure
+           * on the first density iteration (where log10_Peq_prev =
+           * log10_P_cgs), no interpolation necessary */
+
+          log10_n_at_Peq = cooling->nH[i];
+
+        } else {
+
           /* first density bin where log P_eq >= log P_SPH
            * the solution is therefore between log10_Peq_prev and
            * log10_Peq_current */
@@ -466,9 +476,9 @@ double compute_subgrid_property(
                                    (log10_Peq_current - log10_Peq_prev);
 
           /* Interpolate to get the density at equilibrium */
-        const float log10_n_at_Peq =
-            cooling->nH[i - 1] +
+          log10_n_at_Peq = cooling->nH[i - 1] +
                            delta_P_eq * (cooling->nH[i] - cooling->nH[i - 1]);
+        }
 
         /* Return case 5:
          * SPH density and temperatures are projected onto the thermal