// Collision.cc
// This program simulates the interaction of two
// particles colliding.
// Usage:
//  mkdir output
//  ./Collision > output/data1

#include <iostream>
#include <cmath>
using namespace std;

struct particle {
    double x    ; // x coordinate
    double im   ; // inverse mass
    double v    ; // velocity
} ;

void collide (particle p1, particle p2)
// two particles collide elastically
{
  // find the relative velocity
  double velocity_along_line =  p1.v - p2.v;
  // find mass fraction for p1
  double p1f = p1.im / ( p1.im + p2.im ) ;
  double p2f = p2.im / ( p1.im + p2.im ) ;
  // reverse the c.o.m. velocity of each along the line of collision
  p1.v -= 2.0 * p1f * velocity_along_line ;
  p2.v += 2.0 * p2f * velocity_along_line ;
}

void ShowState(double t, particle &a, particle &b, particle &w1, particle &w2) {
  // Displays data in the order time, Wall 1, Wall 2, Position A, Position B
  cout << t << "\t" << w1.x << "\t" << w2.x << "\t" << a.x << "\t" << b.x << endl;
}

void LeapForward(particle &a, particle &b, particle &w1, particle &w2, double &t, double dt){

  double X = t/dt; // Number of iterations

  for(int i=0; i<X; i++){

    a.x += dt * a.v ;
    b.x += dt * b.v ;

  }

  ShowState(t, a, b, w1, w2) ;
}

// Find next collision time
void nct(particle &a, particle &b, particle &w1, particle &w2, double &t, double &Y, double &tab, double &taw, double &tbw, int &collisiontype) {





  if(a.v > b.v){
    tab = (a.x - b.x)/(b.v - a.v) ; // A to B collision
  }
  if(a.v < 0){
    taw = - a.x/a.v ; // A to Wall 1 collision
  }
  if(b.v > 0){
    tbw = (w2.x - b.x)/b.v ; // B to Wall 2 collision
  }
  if(tab < Y){
    Y = tab;
    collisiontype = 1;
  }
  if(taw < Y){
    Y = taw;
    collisiontype = 2;
  }
  if(tbw < Y){
    Y = tbw;
    collisiontype = 3;
  }
  tab = taw = tbw = 1e100;

  t = t + Y ;
}

int main(){

  double tab = 1e100;
  double taw = 1e100;
  double tbw = 1e100;

  particle a;
  particle b;
  particle w1;
  particle w2;
  int L = 10 ; // Length of Box
  double dt = 0.1;
  double t = 0.0;
  double Y = 1000000000000;
  int collisiontype = 0;

  // Set wall conditions
  w1.x = 0;
  w2.x = L;
  w1.v = 0;
  w2.v = 0;
  w1.im = 0;
  w2.im = 0;

  // Set default conditions
  a.x = 3.0;
  a.v = 0.0;
  b.x = 7.0;
  b.v = 2.0;
  a.im = 1.0;
  b.im = 1.0;

  for(int i=0; i<4; i++){
    nct(a, b, w1, w2, t, Y, tab, taw, tbw, collisiontype);
    LeapForward(a, b, w1, w2, t, dt);

    if(collisiontype==1){
      collide(a,b);
    }

    if(collisiontype==2){
      collide(w1, a);
    }

    if(collisiontype==3){
      collide(b, w2);
    }
  }

  return 0 ;
}