#include "common.h"
#include <gsl/gsl_randist.h>

/* Create a sample of `n` points */
sample_t* sample_t_alloc(size_t n, struct par p) {
  sample_t *x = malloc(sizeof(sample_t));
  x->p    = p;
  x->data = gsl_matrix_alloc(n, 2);
  return x;
}

/* Delete a sample */
void sample_t_free(sample_t *x) {
  gsl_matrix_free(x->data);
  free(x);
}


/* `generate_normal(r, n, p)` will create
 * a sample of `n` points, distributed
 * according to a bivariate gaussian distribution
 * of parameters `p`.
 */
sample_t* generate_normal(
  gsl_rng *r, size_t n, struct par *p) {
  sample_t *s = sample_t_alloc(n, *p);

  for (size_t i = 0; i < n; i++) {
    /* Generate a vector (x,y) with
     * a standard (μ = 0) bivariate
     * gaussian PDF.
     */
    double *x = gsl_matrix_ptr(s->data, i, 0);
    double *y = gsl_matrix_ptr(s->data, i, 1);
    gsl_ran_bivariate_gaussian(
      r,
      p->sigma_x, p->sigma_y, p->rho,               
      x, y);

    /* Shift the vector to (x₀,y₀) */
    *x += p->x0;
    *y += p->y0;
  }

  return s;
}