Particle collision

Question

I'm trying to write a code which given the masses, velocity, displacement and radii of N particles in a 1D box will return their positions at small increments of 0.03s. The programme complies. However does not produce the correct values.

#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;
}

There is a very obvious problem when you only have one ball. Radii, mass and velocity all 1. The balls velocity changes at position 10.03 for some reason instead of bouncing of the wall at 20.

I believe the problem is in the "time" function part of the code. This problem is made more obvious by outputting "k" into example.txt in the bonkers function instead of "l".

Thank you very much :)

Answer 1

The problem is this line.

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

You're recalculating the time to collision in each iteration of the while loop and also incrementing ellapsed time, k, and then comparing them. Time until event and elapsed time doesn't make sense to compare. You can either calculate time to collision and increment time, k, until you get there or you can recalculate time to collision and compare with 0 in each iteration.