path: root/src/optics.c

#include <math.h>
#include "optics.h"
#include "misc.h"
#include "quantum_efficiency.h"
#include <stdio.h> /* for fprintf() */
#include <gsl/gsl_integration.h>
#include <gsl/gsl_errno.h> /* for gsl_strerror() */
#include <gsl/gsl_spline.h>

/* Macro to compute the size of a static C array.
 *
 * See https://stackoverflow.com/questions/1598773. */
#define LEN(x) ((sizeof(x)/sizeof(0[x]))/((size_t)(!(sizeof(x) % sizeof(0[x])))))

/* Absorption coefficients for H2O and D2O as a function of wavelength from
 * SNOMAN.
 *
 * Note: These numbers come from prod/media.dat.
 *
 * FIXME: Should I be using the values from media_2_09.dat or media_qoca_*.cmd? */
static double absorption_coefficient_h2o_wavelengths[7] = {250.0, 350.0, 400.0, 450.0, 500.0, 548.0, 578.0};
static double absorption_coefficient_h2o[7] = {1e-5, 1e-5, 1.3e-5, 6.0e-5, 2.0e-4, 5.8e-4, 9.2e-4};

static double absorption_coefficient_d2o_wavelengths[7] = {254.0, 313.0, 366.0, 406.0, 436.0, 548.0, 578.0};
static double absorption_coefficient_d2o[7] = {1e-5, 1e-5, 1e-5, 1e-5, 1e-5, 1e-5, 1e-5};

/* From prod/media.dat for Acrylic (good). */
static double absorption_coefficient_acrylic_wavelengths[10] = {300.0, 310.0, 320.0, 330.0, 340.0, 350.0, 360.0, 380.0, 400.0, 450.0};
static double absorption_coefficient_acrylic[10] = {0.3641, 0.0966, 0.0667, 0.049, 0.035, 0.0261, 0.0194, 0.0115, 0.0086, 0.0071};

static int initialized = 0;

/* D2O absorption length weighted by the Cerenkov spectrum and the quantum
 * efficiency. */
static double absorption_length_d2o_weighted = 0.0;
/* H2O absorption length weighted by the Cerenkov spectrum and the quantum
 * efficiency. */
static double absorption_length_h2o_weighted = 0.0;
/* Acrylic absorption length weighted by the Cerenkov spectrum and the quantum
 * efficiency. */
static double absorption_length_acrylic_weighted = 0.0;


/* From Table 4 in the paper. */
static double A0 = 0.243905091;
static double A1 = 9.53518094e-3;
static double A2 = -3.64358110e-3;
static double A3 = 2.65666426e-4;
static double A4 = 1.59189325e-3;
static double A5 = 2.45733798e-3;
static double A6 = 0.897478251;
static double A7 = -1.63066183e-2;
static double UV = 0.2292020;
static double IR = 5.432937;

static double RIND_H2O_C1 = 1.302;
static double RIND_H2O_C2 = 0.01562;
static double RIND_H2O_C3 = 0.32;

static double RIND_D2O_C1 = 1.302;
static double RIND_D2O_C2 = 0.01333;
static double RIND_D2O_C3 = 0.32;

/* Wavelength of 1 eV particle in nm.
 *
 * From http://pdg.lbl.gov/2018/reviews/rpp2018-rev-phys-constants.pdf. */
static double HC = 1.239841973976e3;

static double gsl_absorption_length_snoman_d2o(double wavelength, void *params)
{
    /* Returns the absorption length in D2O weighted by the quantum efficiency
     * and the Cerenkov spectrum. */
    double qe;

    qe = get_quantum_efficiency(wavelength);

    return qe*get_absorption_length_snoman_d2o(wavelength)/pow(wavelength,2);
}

static double gsl_absorption_length_snoman_h2o(double wavelength, void *params)
{
    /* Returns the absorption length in H2O weighted by the quantum efficiency
     * and the Cerenkov spectrum. */
    double qe;

    qe = get_quantum_efficiency(wavelength);

    return qe*get_absorption_length_snoman_h2o(wavelength)/pow(wavelength,2);
}

static double gsl_absorption_length_snoman_acrylic(double wavelength, void *params)
{
    /* Returns the absorption length in acrylic weighted by the quantum
     * efficiency and the Cerenkov spectrum. */
    double qe;

    qe = get_quantum_efficiency(wavelength);

    return qe*get_absorption_length_snoman_acrylic(wavelength)/pow(wavelength,2);
}

static double gsl_cerenkov(double wavelength, void *params)
{
    /* Returns the quantum efficiency multiplied by the Cerenkov spectrum. */
    double qe;

    qe = get_quantum_efficiency(wavelength);

    return qe/pow(wavelength,2);
}

int optics_init(void)
{
    double norm;
    double result, error;
    size_t nevals;
    int status;
    gsl_integration_cquad_workspace *w;
    gsl_function F;

    w = gsl_integration_cquad_workspace_alloc(100);

    F.function = &gsl_cerenkov;
    F.params = NULL;

    status = gsl_integration_cquad(&F, 200, 800, 0, 1e-2, w, &norm, &error, &nevals);

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

    F.function = &gsl_absorption_length_snoman_d2o;

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

    absorption_length_d2o_weighted = result/norm;

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

    F.function = &gsl_absorption_length_snoman_h2o;

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

    absorption_length_h2o_weighted = result/norm;

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

    F.function = &gsl_absorption_length_snoman_acrylic;

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

    absorption_length_acrylic_weighted = result/norm;

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

    gsl_integration_cquad_workspace_free(w);

    initialized = 1;

    return 0;
}

double get_weighted_absorption_length_snoman_h2o(void)
{
    /* Returns the weighted absorption length in H2O in cm. */

    if (!initialized) {
        fprintf(stderr, "weighted absorption length hasn't been initialized!\n");
        exit(1);
    }

    return absorption_length_h2o_weighted;
}

double get_absorption_length_snoman_h2o(double wavelength)
{
    /* Returns the absorption length in H2O in cm. */
    static gsl_spline *spline;
    static gsl_interp_accel *acc;

    if (!spline) {
        spline = gsl_spline_alloc(gsl_interp_linear, LEN(absorption_coefficient_h2o_wavelengths));
        gsl_spline_init(spline, absorption_coefficient_h2o_wavelengths, absorption_coefficient_h2o, LEN(absorption_coefficient_h2o_wavelengths));
        acc = gsl_interp_accel_alloc();
    }

    if (wavelength < absorption_coefficient_h2o_wavelengths[0])
        wavelength = absorption_coefficient_h2o_wavelengths[0];

    if (wavelength > absorption_coefficient_h2o_wavelengths[LEN(absorption_coefficient_h2o_wavelengths)-1])
        wavelength = absorption_coefficient_h2o_wavelengths[LEN(absorption_coefficient_h2o_wavelengths)-1];

    return 1.0/gsl_spline_eval(spline, wavelength, acc);
}

double get_weighted_absorption_length_snoman_d2o(void)
{
    /* Returns the weighted absorption length in D2O in cm. */

    if (!initialized) {
        fprintf(stderr, "weighted absorption length hasn't been initialized!\n");
        exit(1);
    }

    return absorption_length_d2o_weighted;
}

double get_absorption_length_snoman_d2o(double wavelength)
{
    /* Returns the absorption length in D2O in cm. */
    static gsl_spline *spline;
    static gsl_interp_accel *acc;

    if (!spline) {
        spline = gsl_spline_alloc(gsl_interp_linear, LEN(absorption_coefficient_d2o_wavelengths));
        gsl_spline_init(spline, absorption_coefficient_d2o_wavelengths, absorption_coefficient_d2o, LEN(absorption_coefficient_d2o_wavelengths));
        acc = gsl_interp_accel_alloc();
    }

    if (wavelength < absorption_coefficient_d2o_wavelengths[0])
        wavelength = absorption_coefficient_d2o_wavelengths[0];

    if (wavelength > absorption_coefficient_d2o_wavelengths[LEN(absorption_coefficient_d2o_wavelengths)-1])
        wavelength = absorption_coefficient_d2o_wavelengths[LEN(absorption_coefficient_d2o_wavelengths)-1];

    return 1.0/gsl_spline_eval(spline, wavelength, acc);
}

double get_weighted_absorption_length_snoman_acrylic(void)
{
    /* Returns the weighted absorption length in acrylic in cm. */

    if (!initialized) {
        fprintf(stderr, "weighted absorption length hasn't been initialized!\n");
        exit(1);
    }

    return absorption_length_acrylic_weighted;
}

double get_absorption_length_snoman_acrylic(double wavelength)
{
    /* Returns the absorption length in acrylic in cm. */
    static gsl_spline *spline;
    static gsl_interp_accel *acc;

    if (!spline) {
        spline = gsl_spline_alloc(gsl_interp_linear, LEN(absorption_coefficient_acrylic_wavelengths));
        gsl_spline_init(spline, absorption_coefficient_acrylic_wavelengths, absorption_coefficient_acrylic, LEN(absorption_coefficient_acrylic_wavelengths));
        acc = gsl_interp_accel_alloc();
    }

    if (wavelength < absorption_coefficient_acrylic_wavelengths[0])
        wavelength = absorption_coefficient_acrylic_wavelengths[0];

    if (wavelength > absorption_coefficient_acrylic_wavelengths[LEN(absorption_coefficient_acrylic_wavelengths)-1])
        wavelength = absorption_coefficient_acrylic_wavelengths[LEN(absorption_coefficient_acrylic_wavelengths)-1];

    return 1.0/gsl_spline_eval(spline, wavelength, acc);
}

double get_index(double p, double wavelength, double T)
{
    /* Returns the index of refraction of pure water for a given density,
     * wavelength, and temperature. The density should be in units of g/cm^3,
     * the wavelength in nm, and the temperature in Celsius.
     *
     * See "Refractive Index of Water and Steam as a function of Wavelength,
     * Temperature, and Density" by Schiebener et al. 1990. */
    /* normalize the temperature and pressure */
    wavelength = wavelength/589.0;
    T = (T+273.15)/273.15;

    /* first we compute the right hand side of Equation 7 */
    double c = A0 + A1*p + A2*T + A3*pow(wavelength,2)*T + A4/pow(wavelength,2) + A5/(pow(wavelength,2)-pow(UV,2)) + A6/(pow(wavelength,2)-pow(IR,2)) + A7*pow(p,2);
    c *= p;

    return sqrt((2*c+1)/(1-c));
}

double get_index_snoman_h2o(double wavelength)
{
    /* Returns the index of refraction of light water that SNOMAN uses for a
     * given wavelength.  The wavelength should be in units of nm.
     *
     * The coefficients come from media.dat in SNOMAN version 5.0294 and the
     * formula used to compute the index of refraction comes from the SNOMAN
     * companion page for titles MEDA. */
    double E;

    /* Calculate the energy of the photon in eV. */
    E = HC/wavelength;

    return RIND_H2O_C1 + RIND_H2O_C2*exp(RIND_H2O_C3*E);
}

double get_index_snoman_d2o(double wavelength)
{
    /* Returns the index of refraction of heavy water that SNOMAN uses for a
     * given wavelength.  The wavelength should be in units of nm.
     *
     * The coefficients come from media.dat in SNOMAN version 5.0294 and the
     * formula used to compute the index of refraction comes from the SNOMAN
     * companion page for titles MEDA. */
    double E;

    /* Calculate the energy of the photon in eV. */
    E = HC/wavelength;

    return RIND_D2O_C1 + RIND_D2O_C2*exp(RIND_D2O_C3*E);
}