Line data    Source code

       1             : #include "ReactionEngine.hpp"
       2             : #include "Disc.hpp"
       3             : #include "MathUtils.hpp"
       4             : #include "Reaction.hpp"
       5             : #include "ReactionTable.hpp"
       6             : 
       7             : #include <numbers>
       8             : 
       9             : namespace cell
      10             : {
      11             : 
      12          12 : ReactionEngine::ReactionEngine(const DiscTypeRegistry& discTypeRegistry, const ReactionTable& reactionTable)
      13          12 :     : discTypeRegistry_(discTypeRegistry)
      14             : {
      15          12 :     combineReactionsIntoSingleMaps(reactionTable);
      16          12 : }
      17             : 
      18           1 : Disc ReactionEngine::transformationReaction(Disc* educt, DiscTypeID productID) const
      19             : {
      20           1 :     Disc product(*educt);
      21           1 :     product.setType(productID);
      22           1 :     educt->markDestroyed();
      23             : 
      24           1 :     return product;
      25             : }
      26             : 
      27           1 : std::pair<Disc, Disc> ReactionEngine::decompositionReaction(Disc* educt, DiscTypeID product1ID,
      28             :                                                             DiscTypeID product2ID) const
      29             : {
      30           1 :     double v = mathutils::abs(educt->getVelocity());
      31           1 :     if (v == 0)
      32             :     {
      33           0 :         const auto angle = mathutils::getRandomNumber<double>(0, 2 * std::numbers::pi);
      34           0 :         educt->setVelocity(Vector2d{std::cos(angle), std::sin(angle)});
      35           0 :         v = mathutils::abs(educt->getVelocity());
      36             :     }
      37             : 
      38           1 :     const Vector2d eductNormalizedVelocity = educt->getVelocity() / v;
      39           1 :     const Vector2d n{-eductNormalizedVelocity.y, eductNormalizedVelocity.x};
      40             : 
      41           1 :     Disc product1(*educt);
      42           1 :     Disc product2(product2ID);
      43             : 
      44           1 :     product1.setType(product1ID);
      45           1 :     product1.setVelocity(v * n);
      46             : 
      47           1 :     product2.setPosition(educt->getPosition());
      48           1 :     product2.setVelocity(-v * n);
      49             : 
      50           1 :     const auto R1 = discTypeRegistry_.getByID(product1.getTypeID()).getRadius();
      51           1 :     const auto R2 = discTypeRegistry_.getByID(product2.getTypeID()).getRadius();
      52           1 :     const auto overlap = R1 + R2 + 1e-6; // Discs at same position always have maximum overlap R1 + R2
      53             : 
      54           1 :     product1.move(0.5 * overlap * n);
      55           1 :     product2.move(-0.5 * overlap * n);
      56             : 
      57           1 :     educt->markDestroyed();
      58             : 
      59           2 :     return std::make_pair(std::move(product1), std::move(product2));
      60             : }
      61             : 
      62           5 : Disc ReactionEngine::combinationReaction(Disc* educt1, Disc* educt2, DiscTypeID productID) const
      63             : {
      64           5 :     const double m = discTypeRegistry_.getByID(productID).getMass();
      65           5 :     const double m1 = discTypeRegistry_.getByID(educt1->getTypeID()).getMass();
      66           5 :     const double m2 = discTypeRegistry_.getByID(educt2->getTypeID()).getMass();
      67             : 
      68           5 :     const Vector2d v1 = educt1->getVelocity();
      69           5 :     const Vector2d v2 = educt2->getVelocity();
      70           5 :     Vector2d v = (m1 * v1 + m2 * v2) / m;
      71             : 
      72           5 :     const double kineticEnergyBefore = educt1->getKineticEnergy(m1) + educt2->getKineticEnergy(m2);
      73           5 :     const double kineticEnergyAfter = 0.5 * m * (v.x * v.x + v.y * v.y);
      74           5 :     const double e = 1e-12;
      75             : 
      76           5 :     if (kineticEnergyAfter > e && kineticEnergyBefore > e)
      77           4 :         v *= std::sqrt(kineticEnergyBefore / kineticEnergyAfter);
      78             :     else
      79             :     {
      80           1 :         const auto angle = mathutils::getRandomNumber<double>(0, 2 * std::numbers::pi);
      81           1 :         const double speed = std::sqrt(2 * kineticEnergyBefore / m);
      82           1 :         v = Vector2d{std::cos(angle), std::sin(angle)} * speed;
      83             :     }
      84             : 
      85           5 :     Disc newDisc(productID);
      86           5 :     newDisc.setVelocity(v);
      87           5 :     newDisc.setPosition((educt1->getPosition() + educt2->getPosition()) / 2.0);
      88             : 
      89           5 :     educt1->markDestroyed();
      90           5 :     educt2->markDestroyed();
      91             : 
      92          10 :     return newDisc;
      93             : }
      94             : 
      95           1 : std::pair<Disc, Disc> ReactionEngine::exchangeReaction(Disc* educt1, Disc* educt2, DiscTypeID product1ID,
      96             :                                                        DiscTypeID product2ID) const
      97             : {
      98           1 :     const auto* d1Type = &discTypeRegistry_.getByID(educt1->getTypeID());
      99           1 :     const auto* d2Type = &discTypeRegistry_.getByID(educt2->getTypeID());
     100           1 :     const auto* product1Type = &discTypeRegistry_.getByID(product1ID);
     101           1 :     const auto* product2Type = &discTypeRegistry_.getByID(product2ID);
     102             : 
     103             :     // Sort both product types and educt discs by radius
     104             :     // Now the smallest/largest disc gets the smallest/largest product type
     105             : 
     106           1 :     if (product1Type->getRadius() > product2Type->getRadius())
     107             :     {
     108           0 :         std::swap(product1Type, product2Type);
     109           0 :         std::swap(product1ID, product2ID);
     110             :     }
     111           1 :     if (d1Type->getRadius() > d2Type->getRadius())
     112             :     {
     113           0 :         std::swap(educt1, educt2);
     114           0 :         std::swap(d1Type, d2Type);
     115             :     }
     116             : 
     117             :     // Prefer the assignment that keeps
     118             :     // as many discs as possible with their original type.
     119             : 
     120           1 :     int leaveAsIs = (educt1->getTypeID() == product1ID) + (educt2->getTypeID() == product2ID);
     121           1 :     int swapAgain = (educt1->getTypeID() == product2ID) + (educt2->getTypeID() == product1ID);
     122             : 
     123           1 :     if (swapAgain > leaveAsIs)
     124             :     {
     125           0 :         std::swap(product1Type, product2Type);
     126           0 :         std::swap(product1ID, product2ID);
     127             :     }
     128             : 
     129           1 :     Disc product1(*educt1);
     130           1 :     Disc product2(*educt2);
     131             : 
     132           1 :     product1.scaleVelocity(std::sqrt(d1Type->getMass() / product1Type->getMass()));
     133           1 :     product1.setType(product1ID);
     134             : 
     135           1 :     product2.scaleVelocity(std::sqrt(d2Type->getMass() / product2Type->getMass()));
     136           1 :     product2.setType(product2ID);
     137             : 
     138           1 :     educt1->markDestroyed();
     139           1 :     educt2->markDestroyed();
     140             : 
     141           2 :     return std::make_pair(std::move(product1), std::move(product2));
     142             : }
     143             : 
     144          26 : void ReactionEngine::applyUnimolecularReactions(Disc& disc, double dt, std::vector<Disc>& newDiscs) const
     145             : {
     146          26 :     const Reaction* reaction = selectUnimolecularReaction(disc.getTypeID(), dt);
     147          26 :     if (!reaction)
     148          24 :         return;
     149             : 
     150           2 :     if (reaction->getType() == Reaction::Type::Transformation)
     151             :     {
     152           1 :         auto product = transformationReaction(&disc, reaction->getProduct1());
     153           1 :         newDiscs.push_back(std::move(product));
     154             :     }
     155             :     else
     156             :     {
     157           1 :         auto products = decompositionReaction(&disc, reaction->getProduct1(), reaction->getProduct2());
     158           1 :         newDiscs.push_back(std::move(products.first));
     159           1 :         newDiscs.push_back(std::move(products.second));
     160             :     }
     161             : }
     162             : 
     163          14 : void ReactionEngine::applyBimolecularReactions(const std::vector<CollisionDetector::Collision>& collisions,
     164             :                                                std::vector<Disc>& newDiscs) const
     165             : {
     166          21 :     for (const auto& collision : collisions)
     167             :     {
     168           7 :         if (collision.invalidatedByDestroyedDiscs())
     169           0 :             continue;
     170             : 
     171             :         const Reaction* reaction =
     172           7 :             selectBimolecularReaction(std::minmax(collision.disc->getTypeID(), collision.otherDisc->getTypeID()));
     173           7 :         if (!reaction)
     174           1 :             continue;
     175             : 
     176           6 :         if (reaction->getType() == Reaction::Type::Combination)
     177             :         {
     178           5 :             auto product = combinationReaction(collision.disc, collision.otherDisc, reaction->getProduct1());
     179           5 :             newDiscs.push_back(std::move(product));
     180             :         }
     181             :         else
     182             :         {
     183             :             auto products =
     184           1 :                 exchangeReaction(collision.disc, collision.otherDisc, reaction->getProduct1(), reaction->getProduct2());
     185           1 :             newDiscs.push_back(std::move(products.first));
     186           1 :             newDiscs.push_back(std::move(products.second));
     187             :         }
     188             :     }
     189          14 : }
     190             : 
     191          26 : const Reaction* ReactionEngine::selectUnimolecularReaction(const DiscTypeID& key, double dt) const
     192             : {
     193          26 :     return selectReaction(
     194          26 :         unimolecularReactions_, key, [&](const Reaction& reaction)
     195          54 :         { return mathutils::getRandomNumber<double>(0, 1) <= 1 - std::pow(1 - reaction.getProbability(), dt); });
     196             : }
     197             : 
     198           7 : const Reaction* ReactionEngine::selectBimolecularReaction(const std::pair<DiscTypeID, DiscTypeID>& key) const
     199             : {
     200          14 :     return selectReaction(bimolecularReactions_, key, [](const Reaction& reaction)
     201          20 :                           { return mathutils::getRandomNumber<double>(0, 1) <= reaction.getProbability(); });
     202             : }
     203             : 
     204          12 : void ReactionEngine::combineReactionsIntoSingleMaps(const ReactionTable& reactionTable)
     205             : {
     206          60 :     for (const auto& table : {reactionTable.getTransformations(), reactionTable.getDecompositions()})
     207             :     {
     208          27 :         for (const auto& [educt, reactions] : table)
     209           3 :             unimolecularReactions_[educt].insert(unimolecularReactions_[educt].end(), reactions.begin(),
     210             :                                                  reactions.end());
     211          36 :     }
     212             : 
     213          60 :     for (const auto& table : {reactionTable.getCombinations(), reactionTable.getExchanges()})
     214             :     {
     215          29 :         for (const auto& [educts, reactions] : table)
     216           5 :             bimolecularReactions_[educts].insert(bimolecularReactions_[educts].end(), reactions.begin(),
     217             :                                                  reactions.end());
     218          36 :     }
     219             : 
     220             :     // Random shuffle all reactions to avoid a bias
     221          12 :     std::minstd_rand rng{std::random_device{}()};
     222             : 
     223          15 :     for (auto& [educt, reactions] : unimolecularReactions_)
     224           3 :         std::shuffle(reactions.begin(), reactions.end(), rng);
     225             : 
     226          17 :     for (auto& [educts, reactions] : bimolecularReactions_)
     227           5 :         std::shuffle(reactions.begin(), reactions.end(), rng);
     228          12 : }
     229             : 
     230             : } // namespace cell