#include <cstdlib>
#include <ctime>
#include <iostream>
#include <fstream>
#include <math.h>
#define NUM_POINTS 10000
 
using namespace std;
 
const long double PI = 4.0 * atan (1.0);
ofstream out;
long double data[NUM_POINTS][6];
 
//return random long double between 0 and number specified
//changed denominator to RAND_MAX + 1, supposedly more accurate
long double rand (long double max){
  return (long double) rand () / ((long double) (RAND_MAX+1.0)) * max;
}
 
//generates data for distribution
//also calculates and prints <v_x>, <v_y>, <v_z>
void generate (){
  long double x_total = 0;
  long double y_total = 0;
  long double z_total = 0;
 
 
  //generating theta in column 0 of matrix
  for(int i = 0; i < NUM_POINTS; ++i)
    data[i][0] = rand(2*PI);
 
  //generating phi's in column 1 of matrix
  for (int i = 0; i < NUM_POINTS; ++i){
    for (data[i][1] = rand(1); rand(1) > 2*acos(pow(data[i][1],.5));)
      data[i][1] = rand(1);
  }
 
  //generating speeds in column 2 of matrix
  for (int i = 0; i < NUM_POINTS; ++i){
    for (data[i][2] = rand(4); rand(1) > 4*pow(PI,-.5)*pow(data[i][2],2)*exp(-pow(data[i][2],2));)
      data[i][2] = rand(4);
  }
 
  //deriving x,y,z components in columns 3,4,5 respectively
  //also calculating averages for components
  for (int i = 0; i < NUM_POINTS; ++i){
    data[i][3] = data[i][2] * sin(data[i][1]) * cos(data[i][0]);
    data[i][4] = data[i][2] * sin(data[i][1]) * sin(data[i][0]);
    data[i][5] = data[i][2] * cos(data[i][1]);
    x_total += data[i][3];
    y_total += data[i][4];
    z_total += data[i][5];
  }
 
 
  cout << "Mean for range (0-4): " << endl 
    << "\tx component: " << (x_total/NUM_POINTS) << endl
    << "\ty component: " << (y_total/NUM_POINTS) << endl
    << "\tz component: " << (z_total/NUM_POINTS) << endl;
}
 
int main () {
  //seeding with sys time
  srand ((unsigned) time (0));
 
  generate();
}