Line data    Source code

       1             : #include "Reaction.hpp"
       2             : #include "ExceptionWithLocation.hpp"
       3             : #include "Hashing.hpp"
       4             : 
       5             : #include <numbers>
       6             : 
       7             : namespace cell
       8             : {
       9             : 
      10           4 : bool operator==(const Reaction& a, const Reaction& b)
      11             : {
      12           8 :     return a.getEduct1() == b.getEduct1() && a.hasEduct2() == b.hasEduct2() &&
      13           8 :            (!a.hasEduct2() || a.getEduct2() == b.getEduct2()) && a.getProduct1() == b.getProduct1() &&
      14          12 :            a.hasProduct2() == b.hasProduct2() && (!a.hasProduct2() || a.getProduct2() == b.getProduct2());
      15             : }
      16             : 
      17           5 : std::string toString(const Reaction& reaction, const DiscTypeRegistry& discTypeRegistry)
      18             : {
      19           5 :     std::string result = discTypeRegistry.getByID(reaction.getEduct1()).getName();
      20             : 
      21           5 :     if (reaction.hasEduct2())
      22           2 :         result += " + " + discTypeRegistry.getByID(reaction.getEduct2()).getName();
      23             : 
      24           5 :     result += " -> " + discTypeRegistry.getByID(reaction.getProduct1()).getName();
      25             : 
      26           5 :     if (reaction.hasProduct2())
      27           3 :         result += " + " + discTypeRegistry.getByID(reaction.getProduct2()).getName();
      28             : 
      29           5 :     return result;
      30           0 : }
      31             : 
      32           0 : bool contains(const Reaction& reaction, DiscTypeID discType)
      33             : {
      34           0 :     return reaction.getEduct1() == discType || reaction.getEduct2() == discType || reaction.getProduct1() == discType ||
      35           0 :            reaction.getProduct2() == discType;
      36             : }
      37             : 
      38          13 : Reaction::Type inferReactionType(bool educt2, bool product2)
      39             : {
      40          13 :     if (!educt2 && !product2)
      41           2 :         return Reaction::Type::Transformation;
      42          11 :     else if (educt2 && !product2)
      43           5 :         return Reaction::Type::Combination;
      44           6 :     else if (!educt2 && product2)
      45           4 :         return Reaction::Type::Decomposition;
      46             :     else
      47           2 :         return Reaction::Type::Exchange;
      48             : }
      49             : 
      50          13 : Reaction::Reaction(DiscTypeID educt1, const std::optional<DiscTypeID>& educt2, DiscTypeID product1,
      51          13 :                    const std::optional<DiscTypeID>& product2, double probability)
      52          13 :     : educt1_(educt1)
      53          13 :     , educt2_(educt2)
      54          13 :     , product1_(product1)
      55          13 :     , product2_(product2)
      56          13 :     , type_(inferReactionType(educt2.has_value(), product2.has_value()))
      57             : {
      58          13 :     setProbability(probability);
      59          13 : }
      60             : 
      61          13 : void Reaction::setProbability(double probability)
      62             : {
      63          13 :     if (probability < 0 || probability > 1)
      64           0 :         throw ExceptionWithLocation("Probability must be between 0 and 1");
      65             : 
      66          13 :     probability_ = probability;
      67          13 : }
      68             : 
      69          16 : void Reaction::validate(const DiscTypeRegistry& discTypeRegistry) const
      70             : {
      71          48 :     const auto getMass = [&](DiscTypeID discTypeID) { return discTypeRegistry.getByID(discTypeID).getMass(); };
      72             : 
      73          16 :     const auto eductMassSum = getMass(educt1_) + (educt2_ ? getMass(*educt2_) : 0);
      74          16 :     const auto productMassSum = getMass(product1_) + (product2_ ? getMass(*product2_) : 0);
      75             : 
      76          16 :     if (std::abs(eductMassSum - productMassSum) > 1e-6)
      77           0 :         throw ExceptionWithLocation(toString(*this, discTypeRegistry) +
      78           0 :                                     ": Product- and educt masses need to be identical");
      79             : 
      80          16 :     if (type_ == Type::Transformation && educt1_ == product1_)
      81           0 :         throw ExceptionWithLocation(toString(*this, discTypeRegistry) + ": Educt 1 and product 1 are identical");
      82          16 : }
      83             : 
      84           1 : void Reaction::validateAreaConservation(const DiscTypeRegistry& discTypeRegistry) const
      85             : {
      86           3 :     const auto getArea = [&](DiscTypeID discTypeID)
      87             :     {
      88           3 :         const auto R = discTypeRegistry.getByID(discTypeID).getRadius();
      89           3 :         return std::numbers::pi * R * R;
      90           1 :     };
      91           1 :     const auto eductAreaSum = getArea(educt1_) + (educt2_ ? getArea(*educt2_) : 0);
      92           1 :     const auto productAreaSum = getArea(product1_) + (product2_ ? getArea(*product2_) : 0);
      93             : 
      94           1 :     if (std::abs(eductAreaSum - productAreaSum) > 1e-6)
      95           1 :         throw ExceptionWithLocation(toString(*this, discTypeRegistry) + ": Area not conserved");
      96           0 : }
      97             : 
      98           0 : std::string Reaction::getTypeString() const
      99             : {
     100           0 :     switch (type_)
     101             :     {
     102           0 :     case Reaction::Type::Decomposition: return "Decomposition";
     103           0 :     case Reaction::Type::Transformation: return "Transformation";
     104           0 :     case Reaction::Type::Combination: return "Combination";
     105           0 :     case Reaction::Type::Exchange: return "Exchange";
     106           0 :     case Reaction::Type::None: return "None";
     107           0 :     default: throw ExceptionWithLocation("Invalid reaction type");
     108             :     }
     109             : }
     110             : 
     111             : } // namespace cell