diff options
Diffstat (limited to 'src/optics.c')
| -rw-r--r-- | src/optics.c | 119 |
1 files changed, 119 insertions, 0 deletions
diff --git a/src/optics.c b/src/optics.c index af89802..b01ce7d 100644 --- a/src/optics.c +++ b/src/optics.c @@ -1,6 +1,10 @@ #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() */ /* Absorption coefficients for H2O and D2O as a function of wavelength from * SNOMAN. @@ -14,6 +18,15 @@ static double absorption_coefficient_h2o[7] = {1e-5, 1e-5, 1.3e-5, 6.0e-5, 2.0e- 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}; +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; + /* From Table 4 in the paper. */ static double A0 = 0.243905091; static double A1 = 9.53518094e-3; @@ -39,12 +52,118 @@ static double RIND_D2O_C3 = 0.32; * 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_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); + } + + 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. */ return 1.0/interp1d(wavelength, absorption_coefficient_h2o_wavelengths, absorption_coefficient_h2o, sizeof(absorption_coefficient_h2o_wavelengths)/sizeof(absorption_coefficient_h2o_wavelengths[0])); } +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. */ |