fix: use AKN's formula for rescaling phi

also formatting

Author: yoctoyotta1024

libs/superdrops/collisions/coalbure.hpp

@@ -54,13 +54,6 @@ struct DoCoalBuRe {
   DoBreakup<NFrags> bu; /**< Instance of DoBreakup with specified no. of fragments calculation. */
   Flag coalbure_flag;   /**< Instance of CoalBuReFlag indicating the action to perform. */
 
-  /*
-  rescale random number phi be in desired range of [0, 1] to account for fact that if a
-  collision occurs (i.e. if gamma != 0) then phi lies in range [0, prob - floor(prob)] rather than
-  [0, 1].
-  * _Note:_ This function is assumed to be consitent with collision_gamma(...) and must be.
-  */
-
   /**
    * @brief Rescales a random number phi to be in the desired range [0, 1].
    *
@@ -76,7 +69,12 @@ struct DoCoalBuRe {
    */
   KOKKOS_FUNCTION
   uint64_t rescale_phi(const double prob, const double phi) const {
-    return phi / (prob - Kokkos::floor(prob));
+    const auto denom = prob - Kokkos::floor(prob);
+    if (phi < denom) {
+      return phi / (prob - Kokkos::floor(prob));
+    } else {
+      return (1.0 - phi) / (prob - Kokkos::floor(prob));
+    }
   }
 
   /**
@@ -149,7 +147,7 @@ struct DoCoalBuRe {
    * @return True if the resulting superdroplet is null, otherwise false.
    */
   KOKKOS_INLINE_FUNCTION
-  bool operator()(Superdrop &drop1, Superdrop &drop2, const double prob, const double phi) const;
+  bool operator()(Superdrop& drop1, Superdrop& drop2, const double prob, const double phi) const;
 };
 
 /**
@@ -198,7 +196,7 @@ inline MicrophysicalProcess auto CoalBuRe(const unsigned int interval,
  * @return True if the resulting superdroplet is null, otherwise false.
  */
 template <NFragments NFrags, CoalBuReFlag Flag>
-KOKKOS_FUNCTION bool DoCoalBuRe<NFrags, Flag>::operator()(Superdrop &drop1, Superdrop &drop2,
+KOKKOS_FUNCTION bool DoCoalBuRe<NFrags, Flag>::operator()(Superdrop& drop1, Superdrop& drop2,
                                                           const double prob,
                                                           const double phi) const {
   /* 1. calculate gamma factor for collision  */