/* Copyright (c) 2019, Anthony Latorre <tlatorre at uchicago>
 *
 * This program is free software: you can redistribute it and/or modify it
 * under the terms of the GNU General Public License as published by the Free
 * Software Foundation, either version 3 of the License, or (at your option)
 * any later version.

 * This program is distributed in the hope that it will be useful, but WITHOUT
 * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
 * FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for
 * more details.

 * You should have received a copy of the GNU General Public License along with
 * this program. If not, see <https://www.gnu.org/licenses/>.
 */

#include "solid_angle.h"
#include <math.h>
#include <stdio.h>
#include "optics.h"
#include "muon.h"
#include "mt19937ar.h"
#include <gsl/gsl_sf_gamma.h>
#include "misc.h"
#include "sno_charge.h"
#include <gsl/gsl_integration.h>
#include "path.h"
#include "random.h"
#include "pdg.h"
#include <gsl/gsl_spline.h>
#include "vector.h"
#include <gsl/gsl_statistics.h>
#include "electron.h"
#include "proton.h"
#include "likelihood.h"
#include "id_particles.h"
#include <gsl/gsl_spline2d.h>
#include <gsl/gsl_errno.h> /* for gsl_strerror() */
#include <gsl/gsl_randist.h>
#include <gsl/gsl_cdf.h>
#include "sno.h"
#include "quad.h"
#include "find_peaks.h"
#include <gsl/gsl_sort.h>
#include "util.h"

typedef int testFunction(char *err);

/* Table of some of the tabulated values of the refractive index of water as a
 * function of wavelength and temperature. In all cases, I just used the values
 * for standard atmospheric pressure and assume this corresponds approximately
 * to a density of 1000 kg/m^3.
 *
 * See Table 7 in https://aip.scitation.org/doi/pdf/10.1063/1.555859. */

struct refractive_index_results {
    double p;
    double T;
    double wavelength;
    double n;
} refractive_index_results[] = {
    {1.0, 0 ,  226.50, 1.39468},
    {1.0, 10,  226.50, 1.39439},
    {1.0, 20,  226.50, 1.39353},
    {1.0, 30,  226.50, 1.39224},
    {1.0, 0 ,  404.41, 1.34431},
    {1.0, 10,  404.41, 1.34404},
    {1.0, 20,  404.41, 1.34329},
    {1.0, 30,  404.41, 1.34218},
    {1.0, 0 ,  589.00, 1.33447},
    {1.0, 10,  589.00, 1.33422},
    {1.0, 20,  589.00, 1.33350},
    {1.0, 30,  589.00, 1.33243},
    {1.0, 0 ,  632.80, 1.33321},
    {1.0, 10,  632.80, 1.33296},
    {1.0, 20,  632.80, 1.33224},
    {1.0, 30,  632.80, 1.33118},
    {1.0, 0 , 1013.98, 1.32626},
    {1.0, 10, 1013.98, 1.32604},
    {1.0, 20, 1013.98, 1.32537},
    {1.0, 30, 1013.98, 1.32437},
    {1.0, 0 , 2325.42, 1.27663},
    {1.0, 10, 2325.42, 1.27663},
    {1.0, 20, 2325.42, 1.27627},
    {1.0, 30, 2325.42, 1.27563},
};

/* Table of the values of solid angle for various values of r0/r and L/r.
 *
 * See Table 1 in http://www.umich.edu/~ners312/CourseLibrary/SolidAngleOfADiskOffAxis.pdf. */

struct solid_angle_results {
    double L;
    double r0;
    double omega;
} solid_angle_results[] = {
    {0.5,0.0,3.4732594},
    {0.5,0.2,3.4184435},
    {0.5,0.4,3.2435434},
    {0.5,0.6,2.9185178},
    {0.5,0.8,2.4122535},
    {0.5,1.0,1.7687239},
    {0.5,1.2,1.1661307},
    {0.5,1.4,0.7428889},
    {0.5,1.6,0.4841273},
    {0.5,1.8,0.3287007},
    {0.5,2.0,0.2324189},
    {1.0,0.0,1.8403024},
    {1.0,0.2,1.8070933},
    {1.0,0.4,1.7089486},
    {1.0,0.6,1.5517370},
    {1.0,0.8,1.3488367},
    {1.0,1.0,1.1226876},
    {1.0,1.2,0.9003572},
    {1.0,1.4,0.7039130},
    {1.0,1.6,0.5436956},
    {1.0,1.8,0.4195415},
    {1.0,2.0,0.3257993},
    {1.5,0.0,1.0552591},
    {1.5,0.2,1.0405177},
    {1.5,0.4,0.9975504},
    {1.5,0.6,0.9301028},
    {1.5,0.8,0.8441578},
    {1.5,1.0,0.7472299},
    {1.5,1.2,0.6472056},
    {1.5,1.4,0.5509617},
    {1.5,1.6,0.4632819},
    {1.5,1.8,0.3866757},
    {1.5,2.0,0.3217142},
    {2.0,0.0,0.6633335},
    {2.0,0.2,0.6566352},
    {2.0,0.4,0.6370508},
    {2.0,0.6,0.6060694},
    {2.0,0.8,0.5659755},
    {2.0,1.0,0.5195359},
    {2.0,1.2,0.4696858},
    {2.0,1.4,0.4191714},
    {2.0,1.6,0.3702014},
    {2.0,1.8,0.3243908},
    {2.0,2.0,0.282707}
};

int test_proton_get_energy(char *err)
{
    /* A very simple test to make sure the energy as a function of distance
     * along the track makes sense. Should include more detailed tests later. */
    double T0, T;
    particle *p;

    optics_init();

    /* Assume initial kinetic energy is 1 GeV. */
    T0 = 1000.0;
    p = particle_init(IDP_PROTON, T0, 10000);
    T = particle_get_energy(1e-9,p);

    /* At the beginning of the track we should have roughly the same energy. */
    if (!isclose(T, T0, 1e-5, 0)) {
        sprintf(err, "KE = %.5f, but expected %.5f", T, T0);
        goto err;
    }

    T = particle_get_energy(p->range,p);

    /* At the end of the track the energy should be approximately 0. */
    if (!isclose(T, 0, 1e-5, 1e-5)) {
        sprintf(err, "KE = %.5f, but expected %.5f", T, 0.0);
        goto err;
    }

    particle_free(p);

    return 0;

err:
    particle_free(p);

    return 1;
}

int test_electron_get_energy(char *err)
{
    /* A very simple test to make sure the energy as a function of distance
     * along the track makes sense. Should include more detailed tests later. */
    double T0, T;
    particle *p;

    optics_init();

    /* Assume initial kinetic energy is 1 GeV. */
    T0 = 1000.0;
    p = particle_init(IDP_E_MINUS, T0, 10000);
    T = particle_get_energy(1e-9,p);

    /* At the beginning of the track we should have roughly the same energy. */
    if (!isclose(T, T0, 1e-5, 0)) {
        sprintf(err, "KE = %.5f, but expected %.5f", T, T0);
        goto err;
    }

    T = particle_get_energy(p->range,p);

    /* At the end of the track the energy should be approximately 0. */
    if (!isclose(T, 0, 1e-5, 1e-5)) {
        sprintf(err, "KE = %.5f, but expected %.5f", T, 0.0);
        goto err;
    }

    particle_free(p);

    return 0;

err:
    particle_free(p);

    return 1;
}

int test_muon_get_energy(char *err)
{
    /* A very simple test to make sure the energy as a function of distance
     * along the track makes sense. Should include more detailed tests later. */
    double T0, T;
    particle *p;

    optics_init();

    /* Assume initial kinetic energy is 1 GeV. */
    T0 = 1000.0;
    p = particle_init(IDP_MU_MINUS, T0, 10000);
    T = particle_get_energy(1e-9,p);

    /* At the beginning of the track we should have roughly the same energy. */
    if (!isclose(T, T0, 1e-5, 0)) {
        sprintf(err, "KE = %.5f, but expected %.5f", T, T0);
        goto err;
    }

    T = particle_get_energy(p->range,p);

    /* At the end of the track the energy should be approximately 0. */
    if (!isclose(T, 0, 1e-5, 1e-5)) {
        sprintf(err, "KE = %.5f, but expected %.5f", T, 0.0);
        goto err;
    }

    particle_free(p);

    return 0;

err:
    particle_free(p);

    return 1;
}

int test_proton_get_range(char *err)
{
    /* A very simple test to make sure we read in the PDG table correctly. */
    double value;

    value = proton_get_range(1.0,1.0);

    if (!isclose(value, 2.458E-03, 1e-5, 0)) {
        sprintf(err, "range = %.5f, but expected %.5f", value, 2.458E-03);
        return 1;
    }

    return 0;
}

int test_electron_get_range(char *err)
{
    /* A very simple test to make sure we read in the PDG table correctly. */
    double value;

    value = electron_get_range(1.0,1.0);

    if (!isclose(value, 4.367e-01, 1e-5, 0)) {
        sprintf(err, "range = %.5f, but expected %.5f", value, 4.367e-01);
        return 1;
    }

    return 0;
}

int test_muon_get_range(char *err)
{
    /* A very simple test to make sure we read in the PDG table correctly. */
    double value;

    value = muon_get_range(1.0,1.0);

    if (!isclose(value, 2.071e-3, 1e-5, 0)) {
        sprintf(err, "range = %.5f, but expected %.5f", value, 1.863e-3);
        return 1;
    }

    return 0;
}

int test_proton_get_dEdx(char *err)
{
    /* A very simple test to make sure we read in the PDG table correctly. */
    double value;

    value = proton_get_dEdx(1.0,1.0);

    if (!isclose(value, 2.608E+02, 1e-5, 0)) {
        sprintf(err, "dE/dx = %.5f, but expected %.5f", value, 2.608E+02);
        return 1;
    }

    return 0;
}

int test_electron_get_dEdx(char *err)
{
    /* A very simple test to make sure we read in the PDG table correctly. */
    double value;

    value = electron_get_dEdx(1.0,1.0);

    if (!isclose(value, 1.862, 1e-5, 0)) {
        sprintf(err, "dE/dx = %.5f, but expected %.5f", value, 1.862);
        return 1;
    }

    return 0;
}

int test_muon_get_dEdx(char *err)
{
    /* A very simple test to make sure we read in the PDG table correctly. */
    double value;

    value = muon_get_dEdx(1.0,1.0);

    if (!isclose(value, 5.471, 1e-5, 0)) {
        sprintf(err, "dE/dx = %.5f, but expected %.5f", value, 6.097);
        return 1;
    }

    return 0;
}

int test_refractive_index(char *err)
{
    /* Tests the get_index() function. */
    int i;
    double n;
    struct refractive_index_results result;

    for (i = 0; i < LEN(refractive_index_results); i++) {
        result = refractive_index_results[i];
        n = get_index(result.p, result.wavelength, result.T);

        if (!isclose(n, result.n, 1e-2, 0)) {
            sprintf(err, "n = %.5f, but expected %.5f", n, result.n);
            return 1;
        }
    }

    return 0;
}

int test_solid_angle_approx(char *err)
{
    /* Tests the get_solid_angle_approx() function. */
    int i;
    double pmt[3] = {0,0,0};
    double pos[3] = {0,0,1};
    double n[3] = {0,0,1};
    double r = 1.0;
    double solid_angle;

    solid_angle = get_solid_angle_approx(pos,pmt,n,r);

    if (!isclose(solid_angle, 2*M_PI*(1-1/sqrt(2)), 1e-2, 0)) {
        sprintf(err, "solid angle = %.5f, but expected %.5f", solid_angle, 2*M_PI*(1-1/sqrt(2)));
        return 1;
    }

    for (i = 0; i < LEN(solid_angle_results); i++) {
        pos[0] = solid_angle_results[i].r0*r;
        pos[2] = solid_angle_results[i].L*r;

        solid_angle = get_solid_angle_approx(pos,pmt,n,r);

        if (!isclose(solid_angle, solid_angle_results[i].omega, 1e-2, 0)) {
            sprintf(err, "solid angle = %.5f, but expected %.5f", solid_angle, solid_angle_results[i].omega);
            return 1;
        }
    }

    return 0;
}

int test_solid_angle(char *err)
{
    /* Tests the get_solid_angle() function. */
    int i;
    double pmt[3] = {0,0,0};
    double pos[3] = {0,0,1};
    double n[3] = {0,0,1};
    double r = 1.0;
    double solid_angle;

    solid_angle = get_solid_angle(pos,pmt,n,r);

    if (!isclose(solid_angle, 2*M_PI*(1-1/sqrt(2)), 1e-9, 0)) {
        sprintf(err, "solid angle = %.5f, but expected %.5f", solid_angle, 2*M_PI*(1-1/sqrt(2)));
        return 1;
    }

    for (i = 0; i < LEN(solid_angle_results); i++) {
        pos[0] = solid_angle_results[i].r0*r;
        pos[2] = solid_angle_results[i].L*r;

        solid_angle = get_solid_angle(pos,pmt,n,r);

        if (!isclose(solid_angle, solid_angle_results[i].omega, 1e-4, 0)) {
            sprintf(err, "solid angle = %.5f, but expected %.5f", solid_angle, solid_angle_results[i].omega);
            return 1;
        }
    }

    return 0;
}

static double sno_charge(double q, void *params)
{
    int n = ((int *) params)[0];
    return get_pq(q,n);
}

int test_sno_charge_integral(char *err)
{
    /* Test that the single PE charge distribution integrates to one. */
    double result, error;
    gsl_function F;
    size_t nevals;
    gsl_integration_cquad_workspace *w;
    int params[1];

    w = gsl_integration_cquad_workspace_alloc(100);

    F.function = &sno_charge;
    params[0] = 1;
    F.params = params;

    init_charge();

    gsl_integration_cquad(&F, -10.0, 1000.0, 0, 1e-9, w, &result, &error, &nevals);

    if (!isclose(result, 1.0, 1e-9, 1e-3)) {
        sprintf(err, "integral of single PE charge distribution is %.2f", result);
        goto err;
    }

    gsl_integration_cquad_workspace_free(w);

    return 0;

err:
    gsl_integration_cquad_workspace_free(w);

    return 1;
}

int test_logsumexp(char *err)
{
    /* Tests the logsumexp() function. */
    int i, j;
    double logp[100], mu, result, expected;

    init_genrand(0);

    for (i = 0; i < 100; i++) {
        mu = genrand_real2();

        for (j = 0; j < LEN(logp); j++) {
            logp[j] = -mu + j*log(mu) - gsl_sf_lnfact(j);
        }
        result = logsumexp(logp, LEN(logp));

        expected = 0.0;

        if (!isclose(result, expected, 1e-9, 1e-9)) {
            sprintf(err, "logsumexp(logp) = %.5g, but expected %.5g", result, expected);
            return 1;
        }
    }

    return 0;
}

double getKineticEnergy(double x, void *params)
{
    return 1.0;
}

int test_path(char *err)
{
    /* Test the path code. This just does a basic check that when we create a
     * path without any KL coefficients, it is correctly rotated and
     * positioned. */
    int i, j, k;
    double pos0[3], dir0[3], T0, range, m;
    double beta, t;
    double pos_expected[3], t_expected;
    double pos[3], dir[3];
    double E, mom, beta0;
    double s;
    double theta0;
    path *p;

    T0 = 1.0;
    range = 1.0;
    m = 1.0;

    init_genrand(0);

    for (i = 0; i < 100000; i++) {
        pos0[0] = genrand_real2();
        pos0[1] = genrand_real2();
        pos0[2] = genrand_real2();

        rand_sphere(dir0);

        p = path_init(pos0,dir0,T0,range,1000,0.1,getKineticEnergy,NULL,NULL,NULL,0,m);

        for (j = 0; j < p->len; j++) {
            s = range*j/(p->len-1);
            pos_expected[0] = pos0[0] + dir0[0]*s;
            pos_expected[1] = pos0[1] + dir0[1]*s;
            pos_expected[2] = pos0[2] + dir0[2]*s;

            /* Calculate total energy */
            E = T0 + m;
            mom = sqrt(E*E - m*m);
            beta0 = mom/E;

            t_expected = s/(beta0*SPEED_OF_LIGHT);

            path_eval(p,s,pos,dir,&beta,&t,&theta0);

            for (k = 0; k < 3; k++) {
                if (!isclose(pos[k], pos_expected[k], 1e-9, 1e-9)) {
                    sprintf(err, "path_eval(%.2g) returned pos[%i] = %.5g, but expected %.5g", s, k, pos[k], pos_expected[k]);
                    goto err;
                }
            }

            for (k = 0; k < 3; k++) {
                if (!isclose(dir[k], dir0[k], 1e-9, 1e-9)) {
                    sprintf(err, "path_eval(%.2g) returned dir[%i] = %.5g, but expected %.5g", s, k, dir[k], dir0[k]);
                    goto err;
                }
            }

            if (!isclose(beta, beta0, 1e-9, 1e-9)) {
                sprintf(err, "path_eval(%.2g) returned beta = %.5g, but expected %.5g", s, beta, beta0);
                goto err;
            }

            if (!isclose(t, t_expected, 1e-9, 1e-9)) {
                sprintf(err, "path_eval(%.2g) returned t = %.5g, but expected %.5g", s, t, t_expected);
                goto err;
            }
        }

        path_free(p);
    }

    return 0;

err:
    path_free(p);

    return 1;
}

int test_interp1d(char *err)
{
    /* Tests the interp1d() function. */
    size_t i;
    double xp[100], yp[100], range, x, y, expected;
    gsl_interp_accel *acc;
    gsl_spline *spline;

    range = 1.0;

    init_genrand(0);

    for (i = 0; i < LEN(xp); i++) {
        xp[i] = range*i/(100-1);
        yp[i] = genrand_real2();
    }

    acc = gsl_interp_accel_alloc();
    spline = gsl_spline_alloc(gsl_interp_linear,100);
    gsl_spline_init(spline,xp,yp,100);

    for (i = 0; i < 1000; i++) {
        x = genrand_real2()*range;
        expected = gsl_spline_eval(spline,x,acc);
        y = interp1d(x,xp,yp,100);
        if (!isclose(y, expected, 1e-9, 1e-9)) {
            sprintf(err, "interp1d(%.2g) returned %.5g, but expected %.5g", x, y, expected);
            goto err;
        }
    }

    gsl_interp_accel_free(acc);
    gsl_spline_free(spline);

    return 0;

err:
    gsl_interp_accel_free(acc);
    gsl_spline_free(spline);

    return 1;
}

int test_kahan_sum(char *err)
{
    /* Tests the kahan_sum() function. */
    size_t i;
    double x[100], sum, expected;

    init_genrand(0);

    expected = 0.0;
    for (i = 0; i < LEN(x); i++) {
        x[i] = genrand_real2();
        expected += x[i];
    }

    sum = kahan_sum(x,LEN(x));

    if (!isclose(sum, expected, 1e-9, 1e-9)) {
        sprintf(err, "kahan_sum returned %.5g, but expected %.5g", sum, expected);
        goto err;
    }

    return 0;

err:
    return 1;
}

int test_lnfact(char *err)
{
    /* Tests the lnfact() function. */
    size_t i;
    double x, expected;

    for (i = 0; i < 1000; i++) {
        x = lnfact(i);
        expected = gsl_sf_lnfact(i);

        if (!isclose(x, expected, 1e-9, 1e-9)) {
            sprintf(err, "lnfact(%zu) returned %.5g, but expected %.5g", i, x, expected);
            goto err;
        }
    }

    return 0;

err:
    return 1;
}

int test_ln(char *err)
{
    /* Tests the lnfact() function. */
    size_t i;
    double x, expected;

    for (i = 1; i < 1000; i++) {
        x = ln(i);
        expected = log(i);

        if (!isclose(x, expected, 1e-9, 1e-9)) {
            sprintf(err, "ln(%zu) returned %.5g, but expected %.5g", i, x, expected);
            goto err;
        }
    }

    return 0;

err:
    return 1;
}

int test_get_path_length(char *err)
{
    /* Tests the get_path_length() function. */
    size_t i;
    double pos1[3], pos2[3], tmp[3], r;
    double l1, l2, l1_expected, l2_expected;

    init_genrand(0);

    /* First we test two points within the sphere. */

    for (i = 1; i < 1000; i++) {
        rand_sphere(pos1);
        r = genrand_real2();
        MUL(pos1,r);
        rand_sphere(pos2);
        r = genrand_real2();
        MUL(pos2,r);

        SUB(tmp,pos2,pos1);

        l1_expected = NORM(tmp);
        l2_expected = 0.0;

        get_path_length(pos1,pos2,1.0,&l1,&l2);

        if (!isclose(l1, l1_expected, 1e-9, 1e-9)) {
            sprintf(err, "get_path_length() returned %.5g for l1, but expected %.5g", l1, l1_expected);
            goto err;
        }

        if (!isclose(l2, l2_expected, 1e-9, 1e-9)) {
            sprintf(err, "get_path_length() returned %.5g for l2, but expected %.5g", l2, l2_expected);
            goto err;
        }
    }

    /* Now we test one point at the origin and the other outside of the sphere. */

    pos1[0] = 0.0;
    pos1[1] = 0.0;
    pos1[2] = 0.0;

    for (i = 1; i < 1000; i++) {
        rand_sphere(pos2);
        r = genrand_real2() + 1.0;
        MUL(pos2,r);

        SUB(tmp,pos2,pos1);

        l1_expected = 1.0;
        l2_expected = NORM(tmp) - 1.0;

        get_path_length(pos1,pos2,1.0,&l1,&l2);

        if (!isclose(l1, l1_expected, 1e-9, 1e-9)) {
            sprintf(err, "get_path_length() returned %.5g for l1, but expected %.5g", l1, l1_expected);
            goto err;
        }

        if (!isclose(l2, l2_expected, 1e-9, 1e-9)) {
            sprintf(err, "get_path_length() returned %.5g for l2, but expected %.5g", l2, l2_expected);
            goto err;
        }
    }

    /* Now we test both points outside the sphere such that the ray doesn't
     * intersect the sphere. */

    for (i = 1; i < 1000; i++) {
        pos1[0] = genrand_real2()-0.5;
        pos1[1] = genrand_real2()-0.5;
        pos1[2] = 1.0 + genrand_real2();
        pos2[0] = genrand_real2()-0.5;
        pos2[1] = genrand_real2()-0.5;
        pos2[2] = 1.0 + genrand_real2();

        SUB(tmp,pos2,pos1);

        l1_expected = 0.0;
        l2_expected = NORM(tmp);

        get_path_length(pos1,pos2,1.0,&l1,&l2);

        if (!isclose(l1, l1_expected, 1e-9, 1e-9)) {
            sprintf(err, "get_path_length() returned %.5g for l1, but expected %.5g", l1, l1_expected);
            goto err;
        }

        if (!isclose(l2, l2_expected, 1e-9, 1e-9)) {
            sprintf(err, "get_path_length() returned %.5g for l2, but expected %.5g", l2, l2_expected);
            goto err;
        }
    }

    /* Now we test both points outside the sphere such that the ray intersects
     * the sphere. */

    for (i = 1; i < 1000; i++) {
        /* Pick a random point outside the sphere. */
        rand_sphere(pos1);
        r = genrand_real2() + 1.0;
        MUL(pos1,r);

        COPY(pos2,pos1);
        MUL(pos2,-1.0);
        normalize(pos2);

        r = genrand_real2() + 1.0;
        MUL(pos2,r);

        SUB(tmp,pos2,pos1);

        l1_expected = 2.0;
        l2_expected = NORM(tmp) - 2.0;

        get_path_length(pos1,pos2,1.0,&l1,&l2);

        if (!isclose(l1, l1_expected, 1e-9, 1e-9)) {
            sprintf(err, "get_path_length() returned %.5g for l1, but expected %.5g", l1, l1_expected);
            goto err;
        }

        if (!isclose(l2, l2_expected, 1e-9, 1e-9)) {
            sprintf(err, "get_path_length() returned %.5g for l2, but expected %.5g", l2, l2_expected);
            goto err;
        }
    }

    return 0;

err:
    return 1;
}

int test_mean(char *err)
{
    /* Tests the mean() function. */
    size_t i, j;
    double x[100];
    double mu, expected;

    init_genrand(0);

    for (i = 0; i < 100; i++) {
        for (j = 0; j < LEN(x); j++) {
            x[j] = genrand_real2();
        }

        mu = mean(x,LEN(x));
        expected = gsl_stats_mean(x,1,LEN(x));

        if (!isclose(mu, expected, 0, 1e-9)) {
            sprintf(err, "mean() returned %.5g, but expected %.5g", mu, expected);
            goto err;
        }
    }

    return 0;

err:
    return 1;
}

int test_std(char *err)
{
    /* Tests the std() function. */
    size_t i, j;
    double x[100];
    double sigma, expected, mu;

    init_genrand(0);

    for (i = 0; i < 100; i++) {
        for (j = 0; j < LEN(x); j++) {
            x[j] = genrand_real2();
        }

        sigma = std(x,LEN(x));
        mu = gsl_stats_mean(x,1,LEN(x));
        expected = gsl_stats_sd_with_fixed_mean(x,1,LEN(x),mu);

        if (!isclose(sigma, expected, 0, 1e-9)) {
            sprintf(err, "std() returned %.5g, but expected %.5g", sigma, expected);
            goto err;
        }
    }

    return 0;

err:
    return 1;
}

int test_log_norm(char *err)
{
    /* Tests the log_norm() function. */
    size_t i;
    double x, mu, sigma, logp, expected;

    init_genrand(0);

    for (i = 0; i < 100; i++) {
        x = randn();
        mu = randn();
        sigma = genrand_real2() + 1.0;

        logp = log_norm(x,mu,sigma);
        expected = log(norm(x,mu,sigma));

        if (!isclose(logp, expected, 0, 1e-9)) {
            sprintf(err, "log_norm(%.2g,%.2g,%.2g) returned %.5g, but expected %.5g", x, mu, sigma, logp, expected);
            goto err;
        }
    }

    return 0;

err:
    return 1;
}

int test_trapz(char *err)
{
    /* Tests the trapz() function. */
    size_t i;
    double y[100], integral, expected;

    init_genrand(0);

    for (i = 0; i < LEN(y); i++) {
        y[i] = genrand_real2();
    }

    expected = 0.0;
    for (i = 1; i < LEN(y); i++) {
        expected += (y[i-1] + y[i])/2;
    }

    integral = trapz(y, 1.0, LEN(y));

    if (!isclose(integral, expected, 0, 1e-9)) {
        sprintf(err, "trapz() returned %.5g, but expected %.5g", integral, expected);
        goto err;
    }

    return 0;

err:
    return 1;
}

int test_interp2d(char *err)
{
    /* Tests the interp2d() function. */
    size_t i, j;
    double xp[100], yp[200], zp[20000], zp2[20000], range, x, y, z, expected;
    gsl_interp_accel *xacc;
    gsl_interp_accel *yacc;
    gsl_spline2d *spline;

    range = 1.0;

    init_genrand(0);

    for (i = 0; i < LEN(xp); i++) {
        xp[i] = range*i/(LEN(xp)-1);
    }

    for (i = 0; i < LEN(yp); i++) {
        yp[i] = range*i/(LEN(yp)-1);
    }

    spline = gsl_spline2d_alloc(gsl_interp2d_bilinear,LEN(xp),LEN(yp));

    for (i = 0; i < LEN(xp); i++) {
        for (j = 0; j < LEN(yp); j++) {
            zp[i*LEN(yp)+j] = genrand_real2();
            gsl_spline2d_set(spline, zp2, i, j, zp[i*LEN(yp)+j]);
        }
    }

    xacc = gsl_interp_accel_alloc();
    yacc = gsl_interp_accel_alloc();

    gsl_spline2d_init(spline,xp,yp,zp2,LEN(xp),LEN(yp));

    for (i = 0; i < 1000; i++) {
        x = genrand_real2()*range;
        y = genrand_real2()*range;
        expected = gsl_spline2d_eval(spline,x,y,xacc,yacc);
        z = interp2d(x,y,xp,yp,zp,LEN(xp),LEN(yp));
        if (!isclose(z, expected, 1e-9, 1e-9)) {
            sprintf(err, "interp2d() returned %.5g, but expected %.5g", z, expected);
            goto err;
        }
    }

    gsl_interp_accel_free(xacc);
    gsl_interp_accel_free(yacc);
    gsl_spline2d_free(spline);

    return 0;

err:
    gsl_interp_accel_free(xacc);
    gsl_interp_accel_free(yacc);
    gsl_spline2d_free(spline);

    return 1;
}

static double gsl_electron_get_angular_pdf(double cos_theta, void *params)
{
    double alpha = ((double *) params)[0];
    double beta = ((double *) params)[1];
    double mu = ((double *) params)[2];
    return electron_get_angular_pdf(cos_theta,alpha,beta,mu);
}

int test_electron_get_angular_pdf(char *err)
{
    /* Tests that the electron_get_angular_pdf() function integrates to 1. */
    size_t i;
    double params[3];
    gsl_integration_cquad_workspace *w;
    gsl_function F;
    double result, error;
    int status;
    size_t nevals;
    double T0;

    w = gsl_integration_cquad_workspace_alloc(100);

    F.function = &gsl_electron_get_angular_pdf;
    F.params = params;

    init_genrand(0);

    for (i = 0; i < 100; i++) {
        T0 = genrand_real2()*1000;

        params[0] = electron_get_angular_distribution_alpha(T0);
        params[1] = electron_get_angular_distribution_beta(T0);
        params[2] = genrand_real2();

        status = gsl_integration_cquad(&F, -1, 1, 0, 1e-9, w, &result, &error, &nevals);

        if (status) {
            sprintf(err, "error integrating electron angular distribution: %s\n", gsl_strerror(status));
            goto err;
        }

        if (!isclose(result, 1.0, 0, 1e-3)) {
            sprintf(err, "integral of electron_get_angular_pdf() returned %.5g, but expected %.5g", result, 1.0);
            goto err;
        }
    }

    gsl_integration_cquad_workspace_free(w);

    return 0;

err:
    gsl_integration_cquad_workspace_free(w);

    return 1;
}

int test_gamma_pdf(char *err)
{
    /* Tests the gamma_pdf() function. */
    size_t i;
    double x, k, theta, expected;
    double result;

    init_genrand(0);

    for (i = 0; i < 100; i++) {
        x = genrand_real2();
        k = genrand_real2();
        theta = genrand_real2();

        result = gamma_pdf(x,k,theta);

        expected = gsl_ran_gamma_pdf(x,k,theta);

        if (!isclose(result, expected, 0, 1e-9)) {
            sprintf(err, "gamma_pdf() returned %.5g, but expected %.5g", result, expected);
            goto err;
        }
    }

    return 0;

err:
    return 1;
}

int test_solid_angle2(char *err)
{
    /* Tests the get_solid_angle()2 function. */
    size_t i;
    double pmt[3] = {0,0,0};
    double pos[3] = {0,0,1};
    double n[3] = {0,0,1};
    double r = 1.0;
    double solid_angle;
    double dir[3];
    double L, r0, R;
    double expected;

    init_genrand(0);

    for (i = 0; i < 100; i++) {
        rand_sphere(pmt);
        rand_sphere(n);
        r = genrand_real2();

        dir[0] = pos[0] - pmt[0];
        dir[1] = pos[1] - pmt[1];
        dir[2] = pos[2] - pmt[2];

        L = fabs(dir[0]*n[0] + dir[1]*n[1] + dir[2]*n[2]);
        R = sqrt(dir[0]*dir[0] + dir[1]*dir[1] + dir[2]*dir[2]);
        r0 = sqrt(R*R - L*L);

        solid_angle = get_solid_angle2(L/r,r0/r);

        expected = get_solid_angle(pos,pmt,n,r);

        if (!isclose(solid_angle, expected, 1e-4, 0)) {
            sprintf(err, "solid angle = %.5f, but expected %.5f", solid_angle, expected);
            return 1;
        }
    }

    return 0;
}

int test_solid_angle_lookup(char *err)
{
    /* Tests the get_solid_angle_lookup() function. */
    size_t i;
    double pmt[3] = {0,0,0};
    double pos[3] = {0,0,1};
    double n[3] = {0,0,1};
    double r = 1.0;
    double solid_angle;
    double dir[3];
    double L, r0, R;
    double expected;

    init_genrand(0);

    for (i = 0; i < 100; i++) {
        rand_sphere(pmt);
        MUL(pmt,800);
        rand_sphere(n);
        r = genrand_real2();

        dir[0] = pos[0] - pmt[0];
        dir[1] = pos[1] - pmt[1];
        dir[2] = pos[2] - pmt[2];

        L = fabs(dir[0]*n[0] + dir[1]*n[1] + dir[2]*n[2]);
        R = sqrt(dir[0]*dir[0] + dir[1]*dir[1] + dir[2]*dir[2]);
        r0 = sqrt(R*R - L*L);

        solid_angle = get_solid_angle_lookup(L/r,r0/r);

        expected = get_solid_angle(pos,pmt,n,r);

        if (!isclose(solid_angle, expected, 1e-2, 0)) {
            sprintf(err, "solid angle = %.5f, but expected %.5f", solid_angle, expected);
            return 1;
        }
    }

    return 0;
}

int test_solid_angle_fast(char *err)
{
    /* Tests the get_solid_angle_lookup() function. */
    size_t i;
    double pmt[3] = {0,0,0};
    double pos[3] = {0,0,1};
    double n[3] = {0,0,1};
    double r = 1.0;
    double solid_angle;
    double expected;

    init_genrand(0);

    for (i = 0; i < 100; i++) {
        rand_sphere(pmt);
        MUL(pmt,800);
        rand_sphere(n);
        r = genrand_real2();

        solid_angle = get_solid_angle_fast(pos,pmt,n,r);

        expected = get_solid_angle(pos,pmt,n,r);

        if (!isclose(solid_angle, expected, 1e-2, 0)) {
            sprintf(err, "solid angle = %.5f, but expected %.5f", solid_angle, expected);
            return 1;
        }
    }

    return 0;
}

int test_norm(char *err)
{
    /* Tests the norm() function. */
    size_t i;
    double x, mean, sigma, value, expected;

    init_genrand(0);

    for (i = 0; i < 100; i++) {
        x = genrand_real2();
        mean = genrand_real2();
        sigma = genrand_real2();

        value = norm(x,mean,sigma);

        expected = gsl_ran_ugaussian_pdf((x-mean)/sigma)/sigma;

        if (!isclose(value, expected, 1e-9, 0)) {
            sprintf(err, "norm = %.5f, but expected %.5f", value, expected);
            return 1;
        }
    }

    return 0;
}

int test_norm_cdf(char *err)
{
    /* Tests the norm_cdf() function. */
    size_t i;
    double x, mean, sigma, value, expected;

    init_genrand(0);

    for (i = 0; i < 100; i++) {
        x = genrand_real2();
        mean = genrand_real2();
        sigma = genrand_real2();

        value = norm_cdf(x,mean,sigma);

        expected = gsl_cdf_gaussian_P(x-mean,sigma);

        if (!isclose(value, expected, 1e-9, 0)) {
            sprintf(err, "norm_cdf = %.5f, but expected %.5f", value, expected);
            return 1;
        }
    }

    return 0;
}

static double gsl_time_pdf(double x, void *params)
{
    double mu_noise, mu_indirect, tmean;
    double mu[2], ts[2], ts_sigma[2];

    double *pars = (double *) params;

    mu_noise = pars[0];
    mu_indirect = pars[1];
    mu[0] = pars[2];
    mu[1] = pars[3];
    ts[0] = pars[4];
    ts[1] = pars[5];
    tmean = pars[6];
    ts_sigma[0] = pars[7];
    ts_sigma[1] = pars[8];

    return time_pdf(x,mu_noise,mu_indirect,mu,2,ts,tmean,ts_sigma);
}

int test_time_pdf_norm(char *err)
{
    /* Tests the time_pdf() function. */
    size_t i;
    gsl_integration_cquad_workspace *w;
    gsl_function F;
    double result, error;
    int status;
    size_t nevals;
    double params[9];
    double expected;

    w = gsl_integration_cquad_workspace_alloc(100);

    params[0] = 0.1;
    params[1] = 0.5;
    params[2] = 1.0;
    params[3] = 1.0;
    params[4] = 100.0;
    params[5] = 120.0;
    params[6] = 100.0;
    params[7] = PMT_TTS;
    params[8] = 10.0;

    F.function = &gsl_time_pdf;
    F.params = params;

    init_genrand(0);

    for (i = 0; i < 100; i++) {
        params[0] = genrand_real2();
        params[1] = genrand_real2();
        params[2] = genrand_real2();
        params[3] = genrand_real2();

        status = gsl_integration_cquad(&F, 0, GTVALID, 0, 1e-9, w, &result, &error, &nevals);

        if (status) {
            sprintf(err, "error integrating time PDF: %s\n", gsl_strerror(status));
            goto err;
        }

        expected = 1.0;

        if (!isclose(result, expected, 1e-2, 0)) {
            sprintf(err, "integral of time_pdf = %.5f, but expected %.5f", result, expected);
            goto err;
        }
    }

    gsl_integration_cquad_workspace_free(w);

    return 0;

err:
    gsl_integration_cquad_workspace_free(w);

    return 1;
}

int test_time_cdf(char *err)
{
    /* Tests the time_cdf() function. */
    size_t i;
    gsl_integration_cquad_workspace *w;
    gsl_function F;
    double result, error;
    int status;
    size_t nevals;
    double params[9];
    double expected;

    w = gsl_integration_cquad_workspace_alloc(100);

    params[0] = 0.1;
    params[1] = 0.5;
    params[2] = 1.0;
    params[3] = 1.0;
    params[4] = 100.0;
    params[5] = 120.0;
    params[6] = 100.0;
    params[7] = PMT_TTS;
    params[8] = 10.0;

    F.function = &gsl_time_pdf;
    F.params = params;

    init_genrand(0);

    for (i = 0; i < 100; i++) {
        params[0] = genrand_real2();
        params[1] = genrand_real2();
        params[2] = genrand_real2();
        params[3] = genrand_real2();

        double t = genrand_real2()*GTVALID;

        status = gsl_integration_cquad(&F, 0, t, 0, 1e-9, w, &result, &error, &nevals);

        if (status) {
            sprintf(err, "error integrating time PDF: %s\n", gsl_strerror(status));
            goto err;
        }

        expected = time_cdf(t,params[0],params[1],&params[2],2,&params[4],params[6],&params[7]);

        if (!isclose(result, expected, 1e-2, 0)) {
            sprintf(err, "integral of time_pdf = %.5f, but expected %.5f", result, expected);
            goto err;
        }
    }

    gsl_integration_cquad_workspace_free(w);

    return 0;

err:
    gsl_integration_cquad_workspace_free(w);

    return 1;
}

int test_quad(char *err)
{
    /* Tests the quad fitter without noise. We draw 1000 random hits, compute
     * the time each PMT is hit (without noise) and then make sure that the
     * positions returned by the quad fitter are all within 1 cm and the
     * initial time is within 1 ns. */
    size_t i, j;
    double x0[3];
    double t0;
    const gsl_rng_type *T;
    gsl_rng *r;
    event ev;
    int index[MAX_PMTS];
    int hits[1000];
    size_t valid_pmts;
    double fit_pos[3];
    double pmt_dir[3];
    double fit_t0;
    double wavelength0;
    double n_d2o;

    init_genrand(0);

    T = gsl_rng_default;
    r = gsl_rng_alloc(T);

    load_pmt_info();

    valid_pmts = 0;
    for (i = 0; i < MAX_PMTS; i++) {
        if (pmts[i].pmt_type != PMT_NORMAL) continue;

        index[valid_pmts++] = i;
    }

    wavelength0 = 400.0;
    n_d2o = get_index_snoman_d2o(wavelength0);

    for (i = 0; i < 100; i++) {
        /* Generate a random position within a cube the size of the AV.
         *
         * Note: This does produce some points which are outside of the PSUP,
         * but I think the quad fitter should still be able to fit these. */
        x0[0] = (genrand_real2()*2 - 1)*AV_RADIUS;
        x0[1] = (genrand_real2()*2 - 1)*AV_RADIUS;
        x0[2] = (genrand_real2()*2 - 1)*AV_RADIUS;
        t0 = genrand_real2()*GTVALID;

        /* Zero out all PMTs. */
        for (j = 0; j < LEN(ev.pmt_hits); j++) {
            ev.pmt_hits[j].hit = 0;
            ev.pmt_hits[j].flags = 0;
        }

        /* Choose LEN(hits) random PMTs which are hit. */
        gsl_ran_choose(r,hits,LEN(hits),index,valid_pmts,sizeof(int));

        /* Calculate the time the PMT got hit. */
        for (j = 0; j < LEN(hits); j++) {
            SUB(pmt_dir,pmts[hits[j]].pos,x0);
            ev.pmt_hits[hits[j]].hit = 1;
            ev.pmt_hits[hits[j]].q = genrand_real2();
            ev.pmt_hits[hits[j]].t = t0 + NORM(pmt_dir)*n_d2o/SPEED_OF_LIGHT;
        }

        if (quad(&ev, fit_pos, &fit_t0, 10000)) {
            sprintf(err, "%s", quad_err);
            goto err;
        }

        /* Since there's no noise, these results should be exact. We test to
         * see that all of the positions are within 1 cm and the time is within
         * 1 ns. */

        if (!isclose(fit_pos[0], x0[0], 0, 1.0)) {
            sprintf(err, "quad returned x = %.5f, but expected %.5f", fit_pos[0], x0[0]);
            goto err;
        }

        if (!isclose(fit_pos[1], x0[1], 0, 1.0)) {
            sprintf(err, "quad returned y = %.5f, but expected %.5f", fit_pos[1], x0[1]);
            goto err;
        }

        if (!isclose(fit_pos[2], x0[2], 0, 1.0)) {
            sprintf(err, "quad returned z = %.5f, but expected %.5f", fit_pos[2], x0[2]);
            goto err;
        }

        if (!isclose(fit_t0, t0, 0, 1.0)) {
            sprintf(err, "quad returned t0 = %.5f, but expected %.5f", fit_t0, t0);
            goto err;
        }
    }

    gsl_rng_free(r);

    return 0;

err:
    gsl_rng_free(r);

    return 1;
}

int test_quad_noise(char *err)
{
    /* Tests the quad fitter with noise. We draw 1000 random hits, compute
     * the time each PMT is hit, add a gaussian sample with standard deviation
     * PMT_TTS and then make sure that the positions returned by the quad
     * fitter are all within 1 m and the initial time is within 10 ns. */
    size_t i, j;
    double x0[3];
    double t0;
    const gsl_rng_type *T;
    gsl_rng *r;
    event ev;
    int index[MAX_PMTS];
    int hits[1000];
    size_t valid_pmts;
    double fit_pos[3];
    double pmt_dir[3];
    double fit_t0;
    double wavelength0;
    double n_d2o;

    init_genrand(0);

    T = gsl_rng_default;
    r = gsl_rng_alloc(T);

    load_pmt_info();

    valid_pmts = 0;
    for (i = 0; i < MAX_PMTS; i++) {
        if (pmts[i].pmt_type != PMT_NORMAL) continue;

        index[valid_pmts++] = i;
    }

    wavelength0 = 400.0;
    n_d2o = get_index_snoman_d2o(wavelength0);

    for (i = 0; i < 100; i++) {
        /* Generate a random position within a cube the size of the AV.
         *
         * Note: This does produce some points which are outside of the PSUP,
         * but I think the quad fitter should still be able to fit these. */
        x0[0] = (genrand_real2()*2 - 1)*AV_RADIUS;
        x0[1] = (genrand_real2()*2 - 1)*AV_RADIUS;
        x0[2] = (genrand_real2()*2 - 1)*AV_RADIUS;
        t0 = genrand_real2()*GTVALID;

        /* Zero out all PMTs. */
        for (j = 0; j < LEN(ev.pmt_hits); j++) {
            ev.pmt_hits[j].hit = 0;
            ev.pmt_hits[j].flags = 0;
        }

        /* Choose LEN(hits) random PMTs which are hit. */
        gsl_ran_choose(r,hits,LEN(hits),index,valid_pmts,sizeof(int));

        /* Calculate the time the PMT got hit. */
        for (j = 0; j < LEN(hits); j++) {
            SUB(pmt_dir,pmts[hits[j]].pos,x0);
            ev.pmt_hits[hits[j]].hit = 1;
            ev.pmt_hits[hits[j]].q = genrand_real2();
            ev.pmt_hits[hits[j]].t = t0 + NORM(pmt_dir)*n_d2o/SPEED_OF_LIGHT + randn()*PMT_TTS;
        }

        if (quad(&ev, fit_pos, &fit_t0, 10000)) {
            sprintf(err, "%s", quad_err);
            goto err;
        }

        if (!isclose(fit_pos[0], x0[0], 0, 100.0)) {
            sprintf(err, "quad returned x = %.5f, but expected %.5f", fit_pos[0], x0[0]);
            goto err;
        }

        if (!isclose(fit_pos[1], x0[1], 0, 100.0)) {
            sprintf(err, "quad returned y = %.5f, but expected %.5f", fit_pos[1], x0[1]);
            goto err;
        }

        if (!isclose(fit_pos[2], x0[2], 0, 100.0)) {
            sprintf(err, "quad returned z = %.5f, but expected %.5f", fit_pos[2], x0[2]);
            goto err;
        }

        if (!isclose(fit_t0, t0, 0, 10.0)) {
            sprintf(err, "quad returned t0 = %.5f, but expected %.5f", fit_t0, t0);
            goto err;
        }
    }

    gsl_rng_free(r);

    return 0;

err:
    gsl_rng_free(r);

    return 1;
}

int test_find_peaks_array(char *err)
{
    /* Tests the find_peaks_array() function. */
    double x[4] = {0,1,0,0};
    size_t imax[10], jmax[10], npeaks;

    find_peaks_array(x,2,2,imax,jmax,&npeaks,LEN(imax),0.1);

    if (npeaks != 1) {
        sprintf(err, "number of peaks = %zu, but expected %i", npeaks, 1);
        return 1;
    }

    if (imax[0] != 0) {
        sprintf(err, "imax[0] = %zu, but expected %i", imax[0], 0);
        return 1;
    }

    if (jmax[0] != 1) {
        sprintf(err, "jmax[0] = %zu, but expected %i", jmax[0], 1);
        return 1;
    }

    double y[10][10] = {
        {0,0,2,0,0,0,0,0,0,0},
        {0,0,0,0,0,0,0,0,0,0},
        {4,0,0,0,0,0,0,0,0,5},
        {0,0,0,0,0,0,0,0,0,0},
        {0,0,0,0,0,0,0,0,0,0},
        {0,0,0,0,0,0,0,0,0,0},
        {0,0,0,0,0,0,0,0,0,0},
        {0,0,0,0,0,0,0,0,0,0},
        {0,0,0,0,0,0,0,0,0,0},
        {0,0,3,0,0,0,0,0,0,0},
    };

    find_peaks_array(&y[0][0],10,10,imax,jmax,&npeaks,LEN(imax),0.1);

    if (npeaks != 2) {
        sprintf(err, "number of peaks = %zu, but expected %i", npeaks, 2);
        return 1;
    }

    if (imax[0] != 2) {
        sprintf(err, "imax[0] = %zu, but expected %i", imax[0], 2);
        return 1;
    }

    if (jmax[0] != 9) {
        sprintf(err, "jmax[0] = %zu, but expected %i", jmax[0], 9);
        return 1;
    }

    if (imax[1] != 9) {
        sprintf(err, "imax[1] = %zu, but expected %i", imax[1], 9);
        return 1;
    }

    if (jmax[1] != 2) {
        sprintf(err, "jmax[1] = %zu, but expected %i", jmax[1], 2);
        return 1;
    }

    return 0;
}

int test_ipow(char *err)
{
    /* Tests the ipow() function. */
    size_t result;

    result = ipow(2,2);

    if (result != 4) {
        sprintf(err, "ipow(2,2) returned %zu, but expected %i", result, 4);
        return 1;
    }

    result = ipow(2,3);

    if (result != 8) {
        sprintf(err, "ipow(2,3) returned %zu, but expected %i", result, 8);
        return 1;
    }

    return 0;
}

int test_product(char *err)
{
    /* Tests the product() function. */
    size_t i, j;
    size_t result[1000];

    size_t expected1[4][2] = {{0,0},{0,1},{1,0},{1,1}};

    product(2,2,result);

    for (i = 0; i < 4; i++) {
        for (j = 0; j < 2; j++) {
            if (result[i*2+j] != expected1[i][j]) {
                sprintf(err, "result[%zu,%zu] = %zu, but expected %zu", i, j, result[i*2+j], expected1[i][j]);
                return 1;
            }
        }
    }

    size_t expected2[8][3] = {{0,0,0},{0,0,1},{0,1,0},{0,1,1},{1,0,0},{1,0,1},{1,1,0},{1,1,1}};

    product(2,3,result);

    for (i = 0; i < 8; i++) {
        for (j = 0; j < 3; j++) {
            if (result[i*3+j] != expected2[i][j]) {
                sprintf(err, "result[%zu,%zu] = %zu, but expected %zu", i, j, result[i*3+j], expected2[i][j]);
                return 1;
            }
        }
    }

    return 0;
}

int test_unique_vertices(char *err)
{
    /* Tests the unique_vertices() function. */
    size_t i, j;
    size_t result[1000];
    size_t nvertices;

    int id1[2] = {IDP_E_MINUS, IDP_E_MINUS};
    size_t expected1[3][2] = {{0,0},{0,1},{1,1}};

    unique_vertices(id1,LEN(id1),2,result,&nvertices);

    if (nvertices != 3) {
        sprintf(err, "unique vertices returned nvertices = %zu, but expected %i", nvertices, 3);
        return 1;
    }

    for (i = 0; i < nvertices; i++) {
        for (j = 0; j < 2; j++) {
            if (result[i*2+j] != expected1[i][j]) {
                sprintf(err, "result[%zu,%zu] = %zu, but expected %zu", i, j, result[i*2+j], expected1[i][j]);
                return 1;
            }
        }
    }

    int id2[3] = {IDP_E_MINUS, IDP_MU_MINUS, IDP_TAU_MINUS};
    size_t expected2[8][3] = {{0,0,0},{0,0,1},{0,1,0},{0,1,1},{1,0,0},{1,0,1},{1,1,0},{1,1,1}};

    unique_vertices(id2,LEN(id2),2,result,&nvertices);

    if (nvertices != 8) {
        sprintf(err, "unique vertices returned nvertices = %zu, but expected %i", nvertices, 8);
        return 1;
    }

    for (i = 0; i < 8; i++) {
        for (j = 0; j < 3; j++) {
            if (result[i*3+j] != expected2[i][j]) {
                sprintf(err, "result[%zu,%zu] = %zu, but expected %zu", i, j, result[i*3+j], expected2[i][j]);
                return 1;
            }
        }
    }

    int id3[3] = {IDP_E_MINUS, IDP_E_MINUS, IDP_MU_MINUS};
    size_t expected3[6][3] = {{0,0,0},{0,0,1},{0,1,0},{0,1,1},{1,1,0},{1,1,1}};

    unique_vertices(id3,LEN(id3),2,result,&nvertices);

    if (nvertices != 6) {
        sprintf(err, "unique vertices returned nvertices = %zu, but expected %i", nvertices, 6);
        return 1;
    }

    for (i = 0; i < nvertices; i++) {
        for (j = 0; j < 3; j++) {
            if (result[i*3+j] != expected3[i][j]) {
                sprintf(err, "result[%zu,%zu] = %zu, but expected %zu", i, j, result[i*3+j], expected3[i][j]);
                return 1;
            }
        }
    }

    return 0;
}

int test_find_peaks_highest(char *err)
{
    /* Tests that the find_peaks_array() function returns the *highest* n peaks
     * assuming there are more than n total peaks. */
    size_t imax[2], jmax[2], npeaks;

    double x[10][10] = {
        {0,0,2,0,0,0,0,0,0,0},
        {0,0,0,0,0,0,0,0,0,0},
        {4,0,0,0,0,0,0,0,0,5},
        {0,0,0,0,0,0,0,0,0,0},
        {0,0,0,0,0,0,0,0,0,0},
        {0,0,0,0,0,0,0,0,0,0},
        {0,0,0,0,0,0,0,0,0,0},
        {0,0,0,0,0,0,0,0,0,0},
        {0,0,0,0,0,0,0,0,0,0},
        {0,0,10,0,0,0,0,0,0,10.1},
    };

    find_peaks_array(&x[0][0],10,10,imax,jmax,&npeaks,LEN(imax),0.1);

    if (npeaks != 2) {
        sprintf(err, "number of peaks = %zu, but expected %i", npeaks, 2);
        return 1;
    }

    if (imax[0] != 9) {
        sprintf(err, "imax[0] = %zu, but expected %i", imax[0], 9);
        return 1;
    }

    if (jmax[0] != 9) {
        sprintf(err, "jmax[0] = %zu, but expected %i", jmax[0], 9);
        return 1;
    }

    if (imax[1] != 9) {
        sprintf(err, "imax[0] = %zu, but expected %i", imax[1], 9);
        return 1;
    }

    if (jmax[1] != 2) {
        sprintf(err, "jmax[0] = %zu, but expected %i", jmax[1], 2);
        return 1;
    }

    return 0;
}

int test_find_peaks_sorted(char *err)
{
    /* Tests that the find_peaks_array() function returns peaks in order from
     * highest to lowest. */
    size_t i, j, k;
    size_t imax[10], jmax[10], npeaks;
    double x[10][10];

    init_genrand(0);

    for (i = 1; i < 100; i++) {
        /* Randomly initialize array. */
        for (j = 0; j < 10; j++) {
            for (k = 0; k < 10; k++) {
                x[j][k] = genrand_real2();
            }
        }

        find_peaks_array(&x[0][0],10,10,imax,jmax,&npeaks,LEN(imax),0.1);

        /* Test that each peak is greater than the previous peak. */
        for (j = 1; j < npeaks; j++) {
            if (x[imax[j]][jmax[j]] > x[imax[j-1]][jmax[j-1]]) {
                sprintf(err, "peak %zu is higher than peak %zu", j, j-1);
                return 1;
            }
        }
    }

    return 0;
}

int test_combinations_with_replacement(char *err)
{
    /* Tests the combinations_with_replacement() function. */
    size_t i, j;
    size_t result[100];
    size_t len;

    size_t expected1[3][2] = {{0,0},{0,1},{1,1}};

    combinations_with_replacement(2,2,result,&len);

    if (len != 3) {
        sprintf(err, "combinations_with_replacement() returned %zu combinations but expected %i", len, 3);
        return 1;
    }

    for (i = 0; i < len; i++) {
        for (j = 0; j < 2; j++) {
            if (result[i*2 + j] != expected1[i][j]) {
                sprintf(err, "result[%zu][%zu] = %zu but expected %zu", i, j, result[i*2+j], expected1[i][j]);
                return 1;
            }
        }
    }

    size_t expected2[6][2] = {{0,0},{0,1},{0,2},{1,1},{1,2},{2,2}};

    combinations_with_replacement(3,2,result,&len);

    if (len != 6) {
        sprintf(err, "combinations_with_replacement() returned %zu combinations but expected %i", len, 6);
        return 1;
    }

    for (i = 0; i < len; i++) {
        for (j = 0; j < 2; j++) {
            if (result[i*2 + j] != expected2[i][j]) {
                sprintf(err, "result[%zu][%zu] = %zu but expected %zu", i, j, result[i*2+j], expected1[i][j]);
                return 1;
            }
        }
    }

    return 0;
}

int test_get_dir(char *err)
{
    /* Tests the get_dir() function. */
    size_t i, j;
    double dir[3], expected_dir[3];
    double theta, phi;

    init_genrand(0);

    for (i = 0; i < 100; i++) {
        theta = genrand_real2()*M_PI;
        phi = genrand_real2()*2*M_PI;

        get_dir(dir,theta,phi);

        expected_dir[0] = sin(theta)*cos(phi);
        expected_dir[1] = sin(theta)*sin(phi);
        expected_dir[2] = cos(theta);

        for (j = 0; j < 3; j++) {
            if (!isclose(dir[j],expected_dir[j],0,1e-9)) {
                sprintf(err, "dir[%zu] returned %.2f but expected %.2f", j, dir[j], expected_dir[j]);
                return 1;
            }
        }
    }

    return 0;
}

int test_argmax(char *err)
{
    /* Tests the argmax() function. */
    size_t i, j;
    double a[100];
    size_t max, p;

    init_genrand(0);

    for (i = 0; i < 100; i++) {
        for (j = 0; j < 100; j++)
            a[j] = genrand_real2();

        max = argmax(a,LEN(a));

        gsl_sort_largest_index(&p,1,a,1,LEN(a));

        if (max != p) {
            sprintf(err, "argmax() returned %zu but expected %zu", max, p);
            return 1;
        }
    }

    return 0;
}

int test_argmin(char *err)
{
    /* Tests the argmin() function. */
    size_t i, j;
    double a[100];
    size_t min, p;

    init_genrand(0);

    for (i = 0; i < 100; i++) {
        for (j = 0; j < 100; j++)
            a[j] = genrand_real2();

        min = argmin(a,LEN(a));

        gsl_sort_smallest_index(&p,1,a,1,LEN(a));

        if (min != p) {
            sprintf(err, "argmin() returned %zu but expected %zu", min, p);
            return 1;
        }
    }

    return 0;
}

static double gsl_electron_get_angular_pdf_no_norm(double cos_theta, void *params)
{
    double alpha = ((double *) params)[0];
    double beta = ((double *) params)[1];
    double mu = ((double *) params)[2];
    return electron_get_angular_pdf_no_norm(cos_theta,alpha,beta,mu);
}

static double electron_get_angular_pdf_norm_test(double alpha, double beta, double mu)
{
    double params[3];
    gsl_integration_cquad_workspace *w;
    gsl_function F;
    double result, error;
    int status;
    size_t nevals;

    w = gsl_integration_cquad_workspace_alloc(100);

    F.function = &gsl_electron_get_angular_pdf_no_norm;
    F.params = params;

    params[0] = alpha;
    params[1] = beta;
    params[2] = mu;

    status = gsl_integration_cquad(&F, -1, 1, 0, 1e-9, w, &result, &error, &nevals);

    if (status) {
        fprintf(stderr, "error integrating electron angular distribution: %s\n", gsl_strerror(status));
        exit(1);
    }

    gsl_integration_cquad_workspace_free(w);

    return result;
}

int test_electron_get_angular_pdf_norm(char *err)
{
    /* Tests that the electron_get_angular_pdf_norm() function returns the correct value. */
    size_t i;
    double a, b, mu, result, expected;

    init_genrand(0);

    for (i = 0; i < 100; i++) {
        a = fmax(0.1,genrand_real2());
        b = genrand_real2();
        mu = genrand_real2()*2 - 1;

        result = electron_get_angular_pdf_norm(a,b,mu);
        expected = electron_get_angular_pdf_norm_test(a,b,mu);

        if (!isclose(result, expected, 0, 1e-9)) {
            sprintf(err, "electron_get_angular_pdf_norm() returned %.5g, but expected %.5g", result, expected);
            goto err;
        }
    }

    return 0;

err:
    return 1;
}

int test_fast_acos(char *err)
{
    /* Tests that the fast_acos() function returns values within 0.1% of acos(). */
    size_t i;
    double x, result, expected;

    init_genrand(0);

    for (i = 0; i < 100; i++) {
        x = genrand_real2()*2 - 1;

        result = fast_acos(x);
        expected = acos(x);

        if (!isclose(result, expected, 0, 1e-3)) {
            sprintf(err, "fast_acos() returned %.5g, but expected %.5g", result, expected);
            goto err;
        }
    }

    return 0;

err:
    return 1;
}

int test_fast_sqrt(char *err)
{
    /* Tests that the fast_sqrt() function returns values within 0.1% of sqrt(). */
    size_t i;
    double x, result, expected;

    init_genrand(0);

    for (i = 0; i < 100; i++) {
        x = genrand_real2()*2;

        result = fast_sqrt(x);
        expected = sqrt(x);

        if (!isclose(result, expected, 0, 1e-3)) {
            sprintf(err, "fast_sqrt() returned %.5g, but expected %.5g", result, expected);
            goto err;
        }
    }

    return 0;

err:
    return 1;
}

int test_get_most_likely_mean_pe(char *err)
{
    /* Tests that the get_most_likely_mean_pe() function returns values which
     * are monotonically increasing. */
    size_t i;
    double q, result, last_result;
    size_t n = 100;

    init_charge();

    for (i = 0; i < n; i++) {
        q = i*1000.0/(n-1);

        result = get_most_likely_mean_pe(q);

        if (i > 0 && result < last_result) {
            sprintf(err, "get_most_likely_mean_pe() returned %.5g for q = %.2f, but expected something bigger than %.5g", result, q, last_result);
            goto err;
        }

        last_result = result;
    }

    return 0;

err:
    return 1;
}

struct tests {
    testFunction *test;
    char *name;
} tests[] = {
    {test_solid_angle, "test_solid_angle"},
    {test_solid_angle_approx, "test_solid_angle_approx"},
    {test_refractive_index, "test_refractive_index"},
    {test_muon_get_dEdx, "test_muon_get_dEdx"},
    {test_muon_get_range, "test_muon_get_range"},
    {test_muon_get_energy, "test_muon_get_energy"},
    {test_logsumexp, "test_logsumexp"},
    {test_sno_charge_integral, "test_sno_charge_integral"},
    {test_path, "test_path"},
    {test_interp1d, "test_interp1d"},
    {test_kahan_sum, "test_kahan_sum"},
    {test_lnfact, "test_lnfact"},
    {test_ln, "test_ln"},
    {test_get_path_length, "test_get_path_length"},
    {test_mean, "test_mean"},
    {test_std, "test_std"},
    {test_electron_get_dEdx, "test_electron_get_dEdx"},
    {test_electron_get_range, "test_electron_get_range"},
    {test_electron_get_energy, "test_electron_get_energy"},
    {test_proton_get_dEdx, "test_proton_get_dEdx"},
    {test_proton_get_range, "test_proton_get_range"},
    {test_proton_get_energy, "test_proton_get_energy"},
    {test_log_norm, "test_log_norm"},
    {test_trapz, "test_trapz"},
    {test_interp2d, "test_interp2d"},
    {test_electron_get_angular_pdf, "test_electron_get_angular_pdf"},
    {test_gamma_pdf, "test_gamma_pdf"},
    {test_solid_angle2, "test_solid_angle2"},
    {test_solid_angle_lookup, "test_solid_angle_lookup"},
    {test_solid_angle_fast, "test_solid_angle_fast"},
    {test_norm, "test_norm"},
    {test_norm_cdf, "test_norm_cdf"},
    {test_time_pdf_norm, "test_time_pdf_norm"},
    {test_time_cdf, "test_time_cdf"},
    {test_quad, "test_quad"},
    {test_quad_noise, "test_quad_noise"},
    {test_find_peaks_array, "test_find_peaks_array"},
    {test_ipow, "test_ipow"},
    {test_product, "test_product"},
    {test_unique_vertices, "test_unique_vertices"},
    {test_find_peaks_highest, "test_find_peaks_highest"},
    {test_find_peaks_sorted, "test_find_peaks_sorted"},
    {test_combinations_with_replacement, "test_combinations_with_replacement"},
    {test_get_dir, "test_get_dir"},
    {test_argmax, "test_argmax"},
    {test_argmin, "test_argmin"},
    {test_electron_get_angular_pdf_norm, "test_electron_get_angular_pdf_norm"},
    {test_fast_acos, "test_fast_acos"},
    {test_fast_sqrt, "test_fast_sqrt"},
    {test_get_most_likely_mean_pe, "test_get_most_likely_mean_pe"},
};

int main(int argc, char **argv)
{
    int i;
    char err[256];
    int retval = 0;
    struct tests test;

    for (i = 0; i < LEN(tests); i++) {
        test = tests[i];

        if (!test.test(err)) {
            printf("[\033[92mok\033[0m] %s\n", test.name);
        } else {
            printf("[\033[91mfail\033[0m] %s: %s\n", test.name, err);
            retval = 1;
        }
    }

    return retval;
}