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

using namespace std;

// 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
  float x[11], y[11], theta[11];

  // initilize first element at origin
  x[0] = 0.;
  y[0] = 0.;
  theta[0] = unirandom(-30., 30.);
  
  // loop over all detector planes
  for (int i=1; i<10; i++) {
    
    // planes are located at unit distances
    x[i] = i;

    // extrapolate from last plane 
    y[i] = y[i-1] + tan(M_PI / 180. * theta[i-1])*(x[i]-x[i-1]);

    // scatter by +-2 deg with 50% probability
    if (random() > pow(2., 31.)/2. ) {
      theta[i] = theta[i-1] + unirandom(-2., 2.);
    } else {
      theta[i] = theta[i-1];
    }
      
  }

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

    cout << i 
         << ":  x = " << setw(8) << x[i]
         << "   y = " << setw(8) << y[i]
         << "   theta = " << setw(8) << theta[i]
         << endl;
  }
}