139 lines
3.8 KiB
C
139 lines
3.8 KiB
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;
|
||
}
|