2014-12-02 211 views
0

我試圖寫一個代碼,它給出了1D框中N粒子的質量,速度,位移和半徑將以0.03s的小增量返回它們的位置。 程序符合。但是不會產生正確的值。粒子碰撞

#include <iostream> 
#include <cmath> 
#include <fstream> 
#include <vector> 


using namespace std; 



struct particle { 
    double x; //position of particle 
    double im; //inverse mass 
    double v; //velocity of particle 
    double r; //radius of particle 
}; 

void bonkers(vector<particle> &array, int NN, double t, double ddt, int T); 
double time(vector<particle> &a, int NN, int &whichn); 
void collision(particle &a, particle &b); 
void position(particle &a, double ddt); 

void position(particle &a, double ddt)//position of particles 
{ 
    a.x += ddt * a.v; //position = t * velocity 
} 

void collision(particle &a, particle &b)//velocities after collision 
{ 
    double realativeV = a.v - b.v; 
    double af = a.im/(a.im + b.im); 
    double bf = b.im/(a.im + b.im); 
    a.v -= 2.0 * af * realativeV; 
    b.v += 2.0 * bf * realativeV; 
} 

double time(vector<particle> &a, int NN, int &whichn)//Finds time of first collision. 
{ 
    double dt = 1e100; //really large number 
    for (int n = 0; n < NN - 1 ; n++) 
    { 
     double RelativeV = a[n].v - a[n + 1].v; 
     if (RelativeV > 0)//this means dt won't be 0, no negative times 
     { 
      double Collisiont = ((a[n + 1].x - a[n + 1].r) - (a[n].x + a[n].r))/ RelativeV;//gives time of collision 
      //time of collision worked out as time when the distance between the two balls centre is equal to their combined radi 
      if (Collisiont < dt)//finds smallest possible time of collision 
      { 
       dt = Collisiont; //returns time of first collision 
       whichn = n; //gives which particles collide. Therefore which particles velocities need to change 
      } 
     } 
    } 
    return dt; 
} 

void bonkers(vector<particle> &array, int NN, double t, double ddt, int T) 
{ 

    ofstream myfile; 
    myfile.open("example.txt");//sends all information to file called example.txt 
    int whichn = 0; 
    double k; 
    double l = 0; 
    for (int i = 0; i <= T; i++)//set arbitrary number of collision i.e. T 
    { 
     k = t;//after every collision k reset to 0 as time function works from 0 
     while(k <= time(array, NN, whichn)) 
     { 
      myfile << l; //prints the time 
      k += ddt; //incriments time 
      l += ddt; //increments time, but doesn't reset time so gives total time 

      for (int n = 1; n < NN -1 ; n++) 
      { 
       position(array[n], ddt);//calculates new position 
       myfile << "\t" << array[n].x;//prints the positions of all the particles at said time 
      } 
      myfile << endl; 
     } 
     collision(array[whichn], array[whichn + 1]);//asings new velocities to the two particles which collided 
    } 
    myfile.close(); 

} 

int main() 
{ 
    int N; //number of balls 
    double m; //mass of balls 
    int T = 50; //number of collisions 
    double t = 0.0; //starts from 0 time 
    double ddt = 0.03; //invrements by 0.03 seconds each time 

    cout << "Input number of balls" << endl; 
    cin >> N; 

    vector<particle> a(N + 2); //extra two for the two walls, wall1 at 0 and wall2 at 20 

    for (int k = 1; k <= N; k++) 
    { 
     cout << "Input mass of ball " << k << endl; 
     cin >> m; 
     a[k].im = 1.0/m; //finds inverse mass 
     cout << "Input position of ball " << k << " between 0 and 20" << endl; 
     cin >> a[k].x; 
     cout << "Input speed of ball " << k << endl; 
     cin >> a[k].v; 
     cout << "Input radius of ball " << k << endl; 
     cin >> a[k].r; 
    } 

    //asign wall variables 
    a[0].im = 0; //infinte mass, so walls don't move 
    a[0].r = 0; 
    a[0].x = 0; //wall1 at 0 
    a[0].v = 0; 
    a[N + 1].im = 0;  
    a[N + 1].r = 0; 
    a[N + 1].x = 20; //wall2 at 20 
    a[N + 1].v = 0; 

    bonkers(a, N + 2, t, ddt, T); 
    return 0; 
} 

當你只有一個球時,存在一個非常明顯的問題。半徑,質量和速度全部爲1. 由於某種原因,球速度在位置10.03變化,而不是在20處跳動。

我相信問題出在代碼的「時間」功能部分。 通過在bonkers函數中輸出「k」到example.txt而不是「l」,這個問題變得更加明顯。

非常感謝你:)

+0

您是否使用過調試器和單步執行代碼?如果是這樣,問題在哪裏?如果這是實時的,您是否考慮將值打印到控制檯或日誌文件? – 2014-12-02 22:12:27

+0

你也可以提供你的初始變量嗎?否則,我無法重新創建您的結果... – 2014-12-02 22:13:35

回答

0

問題是這條線。

while(k <= time(array, NN, whichn)) 

你重新計算時間在while循環的每次迭代的碰撞,也增加ellapsed時間,k,然後比較它們。直到事件和經過的時間沒有意義進行比較。您可以計算碰撞時間和增量時間,直到您到達那裏,或者您可以重新計算碰撞時間並在每次迭代中與0進行比較。

+0

非常感謝:D – Emily 2014-12-02 23:47:17

+0

沒問題。我注意到了一些其他的事情:你將std :: vector的大小作爲單獨的變量傳遞,但最好使用std :: vector.size()。使用const將有助於代碼的清晰度,*特別是*如果您正在修改通過引用傳遞的變量(如time()函數)。 – Praxeolitic 2014-12-02 23:50:26