#pragma once

#include <gsl/gsl_histogram.h>
#include <gsl/gsl_vector_double.h>

/* Intensity function parameters */
struct param {
  double a;  // aperture radius
  double k;  // wavenumber
  double e;  // wave amplitude
  double l;  // screen distance
};


/* Intensity of the EM field at a large distance
 * from a circular aperture diffraction of a plane
 * wave.
 *
 * The intensity is given as a function of θ,
 * the diffraction angle. All the other physical
 * parameters are controlled by the structure `param`.
 *
 * See Fraunhöfer diffraction on "Hect - Optics" for
 * a derivation of the formula.
 */
double intensity(double theta, struct param p);


/* `gaussian_kernel(bins, sigma)` generates a convolution
 * kernel for the intensity histogram with `bins` bins.
 *
 * The kernel is a gaussian PDF with μ at the central bin
 * and σ equal to `sigma`. The size is 6% of `bins` and
 * the bin width is the same of the intensity histogram:
 * π/2 / `bins`.
 */
gsl_histogram* gaussian_kernel(size_t bins, double sigma);


/* `gsl_vector_convolve(a, b, res)` computes the linear
 * convolution of two vectors. The resulting vector
 * has size `a->n + b->n - 1` and is stored in `res`.
 */
int gsl_vector_convolve(gsl_vector *a, gsl_vector *b, gsl_vector *res);


/* Convolves two histograms while keeping the
 * original bin edges. This is a simple wrapper
 * function around `gsl_vector_convolve`.
 */
gsl_histogram* histogram_convolve(gsl_histogram *a, gsl_histogram *b);