aboutsummaryrefslogtreecommitdiff
path: root/src/optics.c
blob: 95698653ff29fc4920794fb5cbd08ed5b55a6a67 (plain)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
#include <math.h>
#include "optics.h"

/* 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;

double get_absorption_length_snoman_h2o(double wavelength)
{
    /* Returns the absorption length in H2O in cm. */
    return 1e4;
}

double get_absorption_length_snoman_d2o(double wavelength)
{
    /* Returns the absorption length in D2O in cm. */
    return 1e4;
}

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);
}