#pragma once

#include <math.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_matrix.h>

/* Structure that holds the θ,φ angles 
 * of a single event of the sample.
 */
struct event {
  double th;
  double ph;
};

/* Structure that contains the sample events
 * and its size.
 */
struct sample {
  const size_t size;
  struct event *events;
};

/* This structure contains the result of a
 * minimisation method: estimated parameters,
 * standard error and their covariance.
 */
typedef struct {
  gsl_vector *par;
  gsl_vector *err;
  gsl_matrix *cov;
} min_result;


/* `vector_fprint(fp, x)` prints to file `fp`
 * the components of a gsl_vector.
 */
void vector_fprint(FILE *fp, const gsl_vector *x);


/* `matrix_fprint(fp, x)` prints to file `fp`
 * the components of a gsl_vector.
 */
void matrix_fprint(FILE *fp, const gsl_matrix *x);


/* `vector_map(f, v)` applies the function
 * `f` to each element of the vector `v`, inplace.
 * The function return a GSL exit code to either
 * report success or a failure.
 */
int vector_map(double (*f)(double x), gsl_vector *v);


/* Angular distribution function F(θ,φ; α,β,γ)
 * Takes three parameters α,β,γ and the angles θ,φ.
 */
double distr(const gsl_vector *par,
             const struct event *event);


/* `grad_distr(par, event, grad)`
 * computes the gradient of F(α,β,γ; α,β) as:
 *
 * ⎡ 1/2 - 3/2⋅cos²(θ) ⎤
 * ⎢                   ⎥ 
 * ⎢  sin²(θ)⋅cos(2φ)  ⎥ 3/(4π) 
 * ⎢                   ⎥
 * ⎣ √2⋅sin(2θ)⋅cos(φ) ⎦
 *
 * where `par` is [α,β,γ], `event` is the angle
 * structure. The results is stored in the
 * gsl_vector `grad`
 */
void grad_distr(const gsl_vector *par,
                const struct event *event,
                gsl_vector *grad);


/* Test the compatibility between the original
 * parameters `def` and the estimates `est.par`
 * with standard errors `est.err`.
 */
void compatibility(gsl_vector *def, min_result est);


/* Free a min_result structure.
 */
void min_result_free(min_result x);