analistica/ex-6/common.c

#include <math.h>
#include <gsl/gsl_sf.h>
#include <gsl/gsl_randist.h>
#include <gsl/gsl_blas.h>

#include "common.h"


/* 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) {
  double x = p.k * p.a * sin(theta);
  double R = p.l / cos(theta);

  /* The intensity function has an eliminable
   * discontinuoty at the origin, which must
   * be handled separately.
   */
  if (fabs(x) < 1e-15) {
    return 0.5 * pow(p.e/R * M_PI * pow(p.a, 2), 2);
  }

  double y = 2*M_PI*pow(p.a, 2) * gsl_sf_bessel_J1(x) / x;
  return 0.5 * pow(p.e * y/R, 2);
}


/* `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) {
  /* Set the size to 6% of the histogram.
   * Always choose n even to correctly center the
   * gaussian in the middle.
   */
  size_t n = ((double) bins * 6) / 100;
  n = n % 2 ? n : n+1;

  double dx = (M_PI/2) / bins;

  gsl_histogram *res = gsl_histogram_alloc(n);
  gsl_histogram_set_ranges_uniform(res, 0, dx * n);

  /* The histogram will be such that
   * the maximum falls in the central bin.
   */
  long int offset = (res->n)/2;

  for (long int i = 0; i < (long int)res->n; i++)
    res->bin[i] = gsl_ran_gaussian_pdf(
        ((double)(i - offset)) * dx, sigma * dx);

  return res;
}


/* `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) {
  size_t m = a->size;
  size_t n = b->size;

  /* Vector views for the overlap dot product */
  gsl_vector_view va, vb;

  /* Reverse `b` before convolving */
  gsl_vector_reverse(b);

  double dot;
  size_t start_a, start_b, many;
  for (size_t i = 0; i < m + n - 1; i++) {
    /* Calculate the starting indices
     * of the two vectors a,b overlap.
     */
    if (i < n) start_a = 0;
    else       start_a = i - (n - 1);

    if (i < n) start_b = (n - 1) - i;
    else       start_b = 0;

    /* Calculate how many elements overlap */
         if (i < n)  many = i + 1;
    else if (i >= m) many = n - (i - m) - 1;
    else             many = n;

    /* Take the dot product of the overlap
     * and store it in the resulting histogram
     */
    va = gsl_vector_subvector(a, start_a, many);
    vb = gsl_vector_subvector(b, start_b, many);
    gsl_blas_ddot(&va.vector, &vb.vector, &dot);
    gsl_vector_set(res, i, dot);
  }

  /* Put b back together */
  gsl_vector_reverse(b);

  return GSL_SUCCESS;
}


/* 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) {
  gsl_histogram *res = gsl_histogram_calloc(a->n + b->n - 1);

  /* Set the same edges as `a`*/
  double max = gsl_histogram_max(a);
  double min = gsl_histogram_min(a);
  gsl_histogram_set_ranges_uniform(res, min, max);

  /* Create vector views for everybody */
  gsl_vector_view va  = gsl_vector_view_array(a->bin, a->n);
  gsl_vector_view vb  = gsl_vector_view_array(b->bin, b->n);
  gsl_vector_view vres = gsl_vector_view_array(res->bin, res->n);

  gsl_vector_convolve(&va.vector, &vb.vector, &vres.vector);

  return res;
}