ガスパーティクルシミュレーション衝突計算C++

Question

これはかなり複雑な問題だと思いますが、私はそれを理解できるように十分小さなスペースに収めることができます。私は現在、コードを書いています

シミュレーションボックス内の理想気体粒子。 2つのパーティクルが衝突するかどうかを計算しています。パーティクルが最も近いポイントに到達するまでの時間を計算しています。 （彼らが衝突時に頭を持っている例を使用して）。私は、彼らがその後、どのような時に計算する前に、2個の粒子のためにすべてで衝突するかどうか確認する必要があり、コードのこのセクションでは

とそれらがどのように衝突するなど このように私の2 paricles用：

Main.cpp

私のプログラムの中で

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); 

私がする粒子を定義します。

ヘッダーファイル

は要するにx、y、z値を与える別のクラスです。また、これらを操作するための複数の関数が含まれています。

---

そして私は.cpp内、で助けが必要な部分（ETC coutの開始と文字を無視し、彼らは私のコードは、テストのために出て、簡単なチェックです。）

方程式を考えます：

私はすでに私のためのドット積とモジュラスを行うには、コードを書かれています

ここ

s時間tacで移動した距離です。

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 << "start" << endl; 

    if(modv < 0.0000001){ 
     cout << "a" << endl; 
     result = FLT_MAX; 
    }else{ 
     cout << "b" << 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 < 0) { 
      cout << "c" << endl; 
      result = FLT_MAX; 
     }else{ 
      cout << "d" << endl; 
      Vector s(v.GetX(), v.GetY(), v.GetZ()); 
      s.Scale(tca); 
      cout << getVelocity().GetX() << endl; 
      cout << getVelocity().GetY() << endl; 
      cout << getVelocity().GetZ() << 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 < 4 * radsqr) { 
       cout << "e" << endl; 
       result = FLT_MAX; 
      } else { 
      } 
      cout << "tca: " << tca << endl; 
     } 
    } 
    cout << "fin" << endl; 
    return result; 
} 

私はいくつかの側面を計算するための方程式を持って、tcaは最接近の時間を指します。

私はb > 4 r^2かどうかを確認するために必要なコードで記述されたように、私はいくつかの試みを行い、bアウトのX、YとZコンポーネントを書かれています。しかし、私はごみの回答を得ています。 私は既に間違いや方向づけをしているかどうかを確認するのに助けが必要です。

---

これまでのすべてのコードは、予期したとおりに動作していますが、それぞれをチェックするために複数のテストを記述しました。

あなたが感じるあらゆる情報については、コメントでお知らせください私は残してきたなど

大歓迎任意のヘルプ。

Answer 1

非常に貧弱な言葉遣いの質問に対して私は謝ります。幸運にも、私は自分の答えが当初予想されていたほうがはるかに簡単であることを発見しました。

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 = ((--r0).dot(v))/v.modulusSqr(); 
    double bsqr; 

    double result=0; 

    double rColTestx=r0.GetX()+v.GetX()*tca; 
    double rColTesty=r0.GetY()+v.GetY()*tca; 
    double rColTestz=r0.GetZ()+v.GetZ()*tca; 

    Vector rtColTest(rColTestx, rColTesty, rColTestz); 


    modv=getVelocity().modulus(); 

    cout << "start " << endl; 

    if(modv<0.0000001){ 
     cout << "a" << endl; 
     result=FLT_MAX; 
    }else{ 
     cout << "b" << endl; 

     if (tca < 0) { 
      cout << "c" << endl; 
      result=FLT_MAX; 
     }else{ 
      cout << "d" << endl; 
      Vector s(v.GetX(), v.GetY(), v.GetZ()); 
      s.Scale(tca); 

      cout << getVelocity().GetX() << endl; 
      cout << getVelocity().GetY() << endl; 
      cout << getVelocity().GetZ() << endl; 


      double radsqr= radius*radius; 

      bsqr=rtColTest.modulusSqr(); 

      if (bsqr < 4*radsqr) { 
       cout << "e" << endl; 
       cout << "collision occurs" << endl; 
       result = FLT_MAX; 
      } else { 
       cout << "collision does not occurs" << endl; 
      } 

     } 

    } 
    cout << "fin" << endl; 
    return result; 

} 

大変申し訳ございません。また、FLT_MAXはcfloat libのものです。私は私の質問でこれを統計しなかった。私はこれを調べるために紙で計算したいくつかの例で動作することがわかりました。

明確にすると、return resultとresult=0は任意です。私は後で時間を返すように編集するが、この部分は必要ない、または望んでいなかった。

Answer 2

コードにいくつか間違いがありました。 resultを0またはFLT_MAXと異なる値に設定することはありません。 bsqrも決して計算しません。そして、私はbsqr < 4r^2と他の方法ではない場合、衝突が起こると思います。 （まあ私はr^2の代わりに4r^2しかし、大丈夫と理解していない）。あなたがベクトル実装を隠すので、私は共通のベクトルライブラリを使用しました。とにかく手作りのものは使用しないことをお勧めします。アルマジロやアイゲンを見てみましょう。

ここでは、Eigenの試しに行っています。

#include <iostream> 
#include <limits> 
#include <type_traits> 
#include "Eigen/Dense" 

struct Particle { 
    double radius; 
    Eigen::Vector3d p; 
    Eigen::Vector3d v; 
}; 

template <class FloatingPoint> 
    std::enable_if_t<std::is_floating_point<FloatingPoint>::value, bool> 
     almost_equal(FloatingPoint x, FloatingPoint y, unsigned ulp=1) 
     { 
      FloatingPoint max = std::max(std::abs(x), std::abs(y)); 
      return std::abs(x-y) <= std::numeric_limits<FloatingPoint>::epsilon()*max*ulp; 
     } 

double timeToCollision(const Particle& left, const Particle& right){ 
    Eigen::Vector3d r0 = right.p - left.p; 
    Eigen::Vector3d v = right.v - left.v; 

    double result = std::numeric_limits<double>::infinity(); 

    double vv = v.dot(v); 
    if (!almost_equal(vv, 0.)) { 
     double tca = (-r0).dot(v)/vv; 
     if (tca >= 0) { 
      Eigen::Vector3d s = tca*v; 
      double bb = r0.dot(r0) - s.dot(s); 
      double radius = std::max(left.radius, right.radius); 
      if (bb < 4*radius*radius) 
       result = tca; 
     } 
    } 
    return result; 
} 

int main() 
{ 
    Eigen::Vector3d vp1 {0,0,0}; 
    Eigen::Vector3d vv1 {1,0,0}; 
    Eigen::Vector3d vp2 {12,0,0}; 
    Eigen::Vector3d vv2 {-1,0,0}; 
    Particle p1 {1, vp1, vv1}; 
    Particle p2 {1, vp2, vv2}; 
    std::cout << timeToCollision(p1, p2) << '\n'; 
}