title: Gas Particle simulation Collision calculation C++
body: <p>This I feel is a rather complicated problem, I hope I can fit it in to small enough of a space to make it understandable. I'm presently writing code to</p>
  
  <blockquote>
    <p>simulate Ideal gas particles inside a box. I'm calculating if two particles will collide having calculated the time taken for them to reach their closest point. (using an example where they have head on collision).</p>
  </blockquote>
  
  <p>In this section of code I need to find if they will collide at all for two particles, before then calculating at what time and how they collide etc.
  Thus for my two paricles:</p>
  
  <p><code>Main.cpp</code> </p>
  
  <pre><code>Vector vp1(0,0,0);
  Vector vv1(1,0,0);
  Vector vp2(12,0,0);
  Vector vv2(-1,0,0);
  Particle Particle1(1, vp1, vv1);
  Particle Particle2(1, vp2, vv2);
  Particle1.timeToCollision(Particle2);
  </code></pre>
  
  <p>Within my program I define a particle to be:</p>
  
  <p>Header File</p>
  
  <pre><code>class Particle {
  private:
  
      Vector p;   //position
      Vector v;   //velocity
      double radius;   //radius
  
  public:
  
      Particle();
  
      Particle(double r, const Vector Vecp, const Vector Vecv);
  
      void setPosition(Vector);
      void setVelocity(Vector);
  
      Vector getPosition() const;
      Vector getVelocity() const;
  
      double getRadius() const;
  
      void move(double t);
  
      double timeToCollision(const Particle particle);
  
      void collideParticles(Particle);
  
      ~Particle();
  };
  </code></pre>
  
  <p><code>Vector</code> is another class that in short gives <code>x</code>, <code>y</code>, <code>z</code> values. It also contains multiple functions for manipulating these.</p>
  
  <hr>
  
  <p><strong>And the part that I need help with</strong>, within the <code>.cpp</code> (Ignore the cout start and letters etc, they are simple checks where my code exits for tests.)</p>
  
  <p>Given the equations:</p>
  
  <p><a href="https://i.stack.imgur.com/qs5Eb.png" rel="nofollow noreferrer"><img src="https://i.stack.imgur.com/qs5Eb.png" alt="enter image description here"></a></p>
  
  <p>I have already written code to do the dot product and modulus for me and: </p>
  
  <p><a href="https://i.stack.imgur.com/YEEg7.png" rel="nofollow noreferrer"><img src="https://i.stack.imgur.com/YEEg7.png" alt="b^2=r(0)^2-s^2"></a> </p>
  
  <p>where </p>
  
  <p><a href="https://i.stack.imgur.com/Il1dg.png" rel="nofollow noreferrer"><img src="https://i.stack.imgur.com/Il1dg.png" alt="s^2=((-r(0) dot product v)^2 / modulus squared v)"></a></p>
  
  <p><code>s</code> is distance travelled in time <code>tac</code>.</p>
  
  <pre><code>double Particle::timeToCollision(const Particle particle){
      Vector r2 = particle.getPosition();
      Vector r1 = p;
      Vector v2 = particle.getVelocity();
      Vector v1 = v;
      Vector r0 = r2 - r1;
      Vector v = v2 - v1;
  
      double modv;
      double tca;
      double result = 0;
      double bsqr;
  
      modv = getVelocity().modulus();
  
      cout &lt;&lt; "start" &lt;&lt; endl;
  
      if(modv &lt; 0.0000001){
          cout &lt;&lt; "a" &lt;&lt; endl;
          result = FLT_MAX;
      }else{
          cout &lt;&lt; "b" &lt;&lt; endl;
          tca = ((--r0).dot(v)) / v.modulusSqr();
  
  // -- is an overridden operator that gives the negation ( eg (2, 3, 4) to (-2, -3, -4) )
          if (tca &lt; 0) {
              cout &lt;&lt; "c" &lt;&lt; endl;
              result = FLT_MAX;
          }else{
              cout &lt;&lt; "d" &lt;&lt; endl;
              Vector s(v.GetX(), v.GetY(), v.GetZ());
              s.Scale(tca);
              cout &lt;&lt; getVelocity().GetX() &lt;&lt; endl;
              cout &lt;&lt; getVelocity().GetY() &lt;&lt; endl;
              cout &lt;&lt; getVelocity().GetZ() &lt;&lt; endl;
  
              double radsqr = radius * radius;
              double bx = (r0.GetX() * r0.GetX() - (((r0).dot(v)) *((r0).dot(v)) / v.modulusSqr()));
              double by = (r0.GetY() * r0.GetY() - (((r0).dot(v)) *((r0).dot(v)) / v.modulusSqr()));
              double bz=(r0.GetZ() * r0.GetZ() - (((r0).dot(v)) * ((r0).dot(v)) / v.modulusSqr()));
  
              if (bsqr &lt; 4 * radsqr) {
                  cout &lt;&lt; "e" &lt;&lt; endl;
                  result = FLT_MAX;
              } else {
              }
              cout &lt;&lt; "tca: " &lt;&lt; tca &lt;&lt; endl;
          }
      }
      cout &lt;&lt; "fin" &lt;&lt; endl;
      return result;
  }
  </code></pre>
  
  <p>I have equations for calculating several aspects, <code>tca</code> refers to Time of closest approach. </p>
  
  <p>As written in the code I need to check if <code>b &gt; 4 r^2</code>, I Have made some attempts and written the <code>X</code>, <code>Y</code> and <code>Z</code> components of <code>b</code> out. But I'm getting rubbish answers.
  I just need help to establish if I've already made mistakes or the sort of direction I should be heading. </p>
  
  <hr>
  
  <p>All my code prior to this works as expected and I've written multiple tests for each to check.</p>
  
  <p>Please inform me in a comment for any information you feel I've left out etc.</p>
  
  <p>Any help greatly appreciated. </p>