#include <iostream>
#include <iomanip>
#include <math.h>

using namespace std;

double unirandom(float min, float max);

// state of a particle
class state {

public:
  void set_x_y_theta(float X, float Y, float Theta)
  { x = X; y = Y; theta = Theta; }

  float get_x()     { return x; }
  float get_y()     { return y; }
  float get_theta() { return theta; }

  void extrapolate(float new_x);
  void scatter();

private:
  float x;
  float y;
  float theta;
};

// extrapolate a state to a new x position (=detector plane)
void state::extrapolate(float new_x)
{
  y += tan(M_PI / 180. * theta) * (new_x-x);
  x = new_x;
  // theta remains unchanged
}


// scatter particle with probablity 50% by +/- 2 deg
void state::scatter()
{
    if (random() > pow(2., 31.)/2. ) {
      theta += unirandom(-2., 2.);
    }
}


// function to calculate a uniformly distributed random number 
// between min and max
double unirandom(float min, float max)
{
  return min + (max-min)*random() / (pow(2., 31.)-1.);
}

int main()
{

  // initialize random number generator with time, to produce different
  // results every time
  srandom(time(NULL));

  // allocate arrays to store positions and angles of particle at 
  // origin and all planes
  state states[11];

  // initilize first element at origin
  states[0].set_x_y_theta(0., 0., unirandom(-30., 30.));
  
  // loop over all detector planes
  for (int i=1; i<10; i++) {
    
    // copy state to next array position
    states[i] = states[i-1];

    // extrapolate track
    states[i].extrapolate(i);

    // scatter
    states[i].scatter();

  }

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

    cout << i 
         << ":  x = " << setw(9) << states[i].get_x()
         << "   y = " << setw(9) << states[i].get_y()
         << "   theta = " << setw(9) << states[i].get_theta()
         << endl;
  }
}