#include <math.h>
#include <stdlib.h>
#include <string.h>
#include <gsl/gsl_rng.h>

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

  // Set default options.
  //
  size_t N = 50000;          // number of events.
  size_t n = 50;             // number of bins.
  double p_max = 10;         // maximum value of momentum module.

  // Process CLI arguments.
  //
  for (size_t i = 1; i < argc; i++)
    {
         if (!strcmp(argv[i], "-N")) N     = atol(argv[++i]);
    else if (!strcmp(argv[i], "-n")) n     = atol(argv[++i]);
    else if (!strcmp(argv[i], "-p")) p_max = atof(argv[++i]);
    else
      {
      fprintf(stderr, "Usage: %s -[hiIntp]\n", argv[0]);
      fprintf(stderr, "\t-h\tShow this message.\n");
      fprintf(stderr, "\t-N integer\tThe number of events to generate.\n");
      fprintf(stderr, "\t-n integer\tThe number of bins of the histogram.\n");
      fprintf(stderr, "\t-p float\tThe maximum value of momentum.\n");
      return EXIT_FAILURE;
      }
    }

  // Initialize an RNG.
  //
  gsl_rng_env_setup();
  gsl_rng *r = gsl_rng_alloc(gsl_rng_default);

  // Generate the angle θ uniformly distributed on a sphere using the
  // inverse transform:
  //
  //   θ = acos(1 - 2X)
  //
  // where X is a random uniform variable in [0,1), and the module p of
  // the vector:
  //
  //   p² = p_v² + p_h²
  // 
  // uniformly distributed between 0 and p_max. The two components are
  // then computed as:
  //
  //   p_v = p⋅cos(θ)
  //   p_h = p⋅sin(θ)
  //
  // The histogram will be updated this way.
  // The j-th bin where p_h goes in is given by:
  //
  //   step = p_max / n
  //   j = floor(p_h / step)
  //
  // Thus an histogram was created and a structure containing the number of
  // entries in each bin and the sum of |p_v| in each of them is created and
  // filled while generating the events.
  //
  struct bin
    {
    size_t amo;  // Amount of events in the bin.
    double sum;  // Sum of |p_v|s of all the events in the bin.
    };
  struct bin *histo = calloc(n, sizeof(struct bin));

  // Some useful variables.
  //
  double step = p_max / n;
  struct bin *b;
  double theta;
  double p;
  double p_v;
  double p_h;
  size_t j;

  for (size_t i = 0; i < N; i++)
    {
    // Generate the event.
    // 
    theta = acos(1 - 2*gsl_rng_uniform(r));
    p = p_max * gsl_rng_uniform(r);

    // Compute the components.
    //
    p_v = p * cos(theta);
    p_h = p * sin(theta);

    // Update the histogram.
    //
    j = floor(p_h / step);
    b = &histo[j];
    b->amo++;
    b->sum += fabs(p_v);
    }

  // Compute the mean value of each bin and print it to stodut
  // together with other useful things to make the histogram.
  //
  printf("bins: \t%ld\n", n);
  printf("step: \t%.5f\n", step);
  for (size_t i = 0; i < n; i++)
    {
    histo[i].sum = histo[i].sum / histo[i].amo;
    printf("\n%.5f", histo[i].sum);
    };

  // free memory
  gsl_rng_free(r);
  free(histo);

  return EXIT_SUCCESS;
  }