// Scattering in 2-D  //

                   // Michel Vallieres           //

				//Notation [0] --> x     [1] --> y	[2] --> vx	[3] --> vy

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

                // dimension of system
#define ndim 4
                // model parameters
int    NS;
double k, b, xs[3], ys[3], mass;


                // total energy //
double total_energy( double f[] )
{
  double energy, x, y, vx, vy, temp;
  int    source;

  vx = f[2];
  vy = f[3];
  energy = 0.5*mass*( vx*vx + vy*vy );
  for ( source=0; source<NS ; source++ )
    {
      x = f[0] - xs[source];
      y = f[1] - ys[source];
      temp = k * exp( -( x*x + y*y )/(b*b) );
      energy = energy + 0.5*b*b*temp;
    }

  return energy;
}


void define_model()
{
  k = 1.0;
  b = 0.25;
  mass = 1.0;
  NS = 3;
  xs[0] = 0;
  ys[0] = 0.5;
  xs[1] = 0;
  ys[1] = -0.5;
  xs[2] = sqrt(3.0)/2.0;
  ys[2] = 0.0;
}




void rhs_function(double t, double f_t[], double rhsf[] ) 
{
	int source;
	double x, y, temp;
	
	rhsf[0] = f_t[2] ;	//Equation for x
	rhsf[1] = f_t[3] ;	//Equation for y

	rhsf[2] = 0.0;		//Equation for vx
	rhsf[3] = 0.0;		//Equation for vy

  	for ( source=0; source<NS ; source++ )
    	{
      		x = f_t[0] - xs[source];
     		y = f_t[1] - ys[source];
      		temp = k * exp( -( x*x + y*y )/(b*b) );
		rhsf[2] = rhsf[2] + x*temp;
		rhsf[3] = rhsf[3] + y*temp;
    	}

	rhsf[2] = rhsf[2]/mass;
	rhsf[3] = rhsf[3]/mass;
	
}
                // Runge-Kutta 4th order time step     //
                // f_t[] is the function to advance    //
                // requires a global definition for    //
                // int ndim: dimension of system       //
                // void rhs_function (t, f, rhsf )     //
                                                                                              
void rk4_step ( double t, double dt, double f_t[] )
{
  double k1[ndim], k2[ndim], k3[ndim],
            k4, yt[ndim], rhsf[ndim];
  int i;
                                                                                              
  rhs_function( t, f_t, rhsf );
  for ( i=0 ; i<ndim ; i++ )
    {
      k1[i] = dt * rhsf[i];
      yt[i] = f_t[i] + 0.5*k1[i];
    }
  rhs_function( t+0.5*dt, yt, rhsf );
  for ( i=0 ; i<ndim ; i++ )
    {
      k2[i] = dt * rhsf[i];
      yt[i] = f_t[i] + 0.5*k2[i];
    }
  rhs_function( t+0.5*dt, yt, rhsf );
  for ( i=0 ; i<ndim ; i++ )
    {
      k3[i] = dt * rhsf[i];
      yt[i] = f_t[i] + k3[i];
    }
  rhs_function( t+dt, yt, rhsf );
  for ( i=0 ; i<ndim ; i++ )
    {
      k4 = dt * rhsf[i];
      f_t[i] = f_t[i] + ( 0.5*k1[i] + k2[i] + k3[i] + 0.5*k4 ) / 3.0;
    }
}

int main( int argc, char * argv[] )
{

  double dt, t, t_min, t_max, asymptotic;
  double energy;
  double acc;
  double energy_0;
  int    i;
  int init_from_file;
  double f_t[ndim] ;
                           // time grid
  t_min = 0.0;
  t_max = 250.0;
  dt = 0.01;
 			   //Big radius, asymptotic region

  asymptotic = 4.1 * 4.1;

                           // model parameters
  define_model();
                           // initial values
  t = t_min;
  f_t[0] = -4.0;
  f_t[1] = 0.05;
  f_t[2] = 0.15;
  f_t[3] = 0.0;
                           // reset parameters

  init_from_file = 0;
  for ( i=1; i<argc; i++ )
    {
      if (argv[i][0] == '-') {
	switch (argv[i][1]) {
	case 'd':   dt = atof(argv[++i]); break;
	case 't':   t_max = atof(argv[++i]); break;
	case 'b':   f_t[1] = atof(argv[++i]); break;
	case 'i':   init_from_file = 1; break;
	}
      }
    }

  for(;;){

  if(cin.eof())
	exit(0);

  if(init_from_file == 1)
	cin >>  f_t[0] >> f_t[1] >> f_t[2] >> f_t[3];

  

                           // energy check
  energy = total_energy(f_t);
  energy_0 = energy;
                           // print results

	
	cout << t << ' ' << f_t[0] << ' ' << f_t[1] << ' ' << f_t[2] << ' ' << f_t[3] << ' ' << energy 
             << ' ' << energy_0 << endl ;
                           // loop over time
  while ( t < t_max )
    {
							//Runge-Kutta 4 Step
    	rk4_step (t, dt, f_t) ;

                           // time grid
      t = t + dt;

                           // energy check
      energy = total_energy(f_t);

                           // print results
	cout << t << ' ' << f_t[0] << ' ' << f_t[1] << ' ' << f_t[2] << ' ' << f_t[3] << ' ' << energy 
             << ' ' << energy_0 << endl ;

			  //Reach Asymptotic region
	if (f_t[0]*f_t[0] + f_t[1]*f_t[1] > asymptotic)
		break;
    }
  //break after all cases done

  if(init_from_file == 0)
	break;
  }
}