diff options
Diffstat (limited to 'src/optics.c')
| -rw-r--r-- | src/optics.c | 112 |
1 files changed, 110 insertions, 2 deletions
diff --git a/src/optics.c b/src/optics.c index 2d28ba4..4a0722a 100644 --- a/src/optics.c +++ b/src/optics.c @@ -7,6 +7,10 @@ #include <gsl/gsl_errno.h> /* for gsl_strerror() */ #include <gsl/gsl_spline.h> #include "misc.h" +#include "db.h" +#include "dict.h" + +char optics_err[256]; /* Absorption coefficients for H2O and D2O as a function of wavelength from * SNOMAN. @@ -15,6 +19,10 @@ static double absorption_coefficient_h2o_wavelengths[6] = {337.0, 365.0, 386.0, 420.0, 500.0, 620.0}; static double absorption_coefficient_h2o[6] = {8.410e-05, 2.880e-06, 1.431e-04, 1.034e-04, 5.239e-04, 2.482e-03}; +static double rayleigh_coefficient_wavelengths[51]; +static double rayleigh_coefficient_h2o[51]; +static double rayleigh_coefficient_d2o[51]; + static double absorption_coefficient_d2o_wavelengths[6] = {337.0, 365.0, 386.0, 420.0, 500.0, 620.0}; static double absorption_coefficient_d2o[6] = {6.057e-05, 2.680e-05, 3.333e-05, 3.006e-05, 2.773e-05, 2.736e-05}; @@ -32,7 +40,12 @@ 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; - +/* D2O Rayleigh scattering length weighted by the Cerenkov spectrum and the + * quantum efficiency. */ +static double rayleigh_scattering_length_d2o_weighted = 0.0; +/* H2O Rayleigh scattering length weighted by the Cerenkov spectrum and the + * quantum efficiency. */ +static double rayleigh_scattering_length_h2o_weighted = 0.0; /* From Table 4 in the paper. */ static double A0 = 0.243905091; @@ -59,6 +72,28 @@ 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_rayleigh_scattering_length_snoman_d2o(double wavelength, void *params) +{ + /* Returns the Rayleigh scattering length in D2O weighted by the quantum efficiency + * and the Cerenkov spectrum. */ + double qe; + + qe = get_quantum_efficiency(wavelength); + + return qe*get_rayleigh_scattering_length_snoman_d2o(wavelength)/pow(wavelength,2); +} + +static double gsl_rayleigh_scattering_length_snoman_h2o(double wavelength, void *params) +{ + /* Returns the Rayleigh scattering length in H2O weighted by the quantum efficiency + * and the Cerenkov spectrum. */ + double qe; + + qe = get_quantum_efficiency(wavelength); + + return qe*get_rayleigh_scattering_length_snoman_h2o(wavelength)/pow(wavelength,2); +} + static double gsl_absorption_length_snoman_d2o(double wavelength, void *params) { /* Returns the absorption length in D2O weighted by the quantum efficiency @@ -102,14 +137,16 @@ static double gsl_cerenkov(double wavelength, void *params) return qe/pow(wavelength,2); } -int optics_init(void) +int optics_init(dict *db) { + int i; double norm; double result, error; size_t nevals; int status; gsl_integration_cquad_workspace *w; gsl_function F; + dbval *rspr; w = gsl_integration_cquad_workspace_alloc(100); @@ -156,6 +193,41 @@ int optics_init(void) exit(1); } + rspr = get_bank(db, "RSPR", 1); + + if (!rspr) { + sprintf(optics_err, "failed to load RSPR bank\n"); + return -1; + } + + for (i = 0; i < 51; i++) { + rayleigh_coefficient_wavelengths[i] = rspr[30+i*3].u32; + rayleigh_coefficient_d2o[i] = rspr[30+i*3+1].f; + rayleigh_coefficient_h2o[i] = rspr[30+i*3+2].f; + } + + F.function = &gsl_rayleigh_scattering_length_snoman_d2o; + + status = gsl_integration_cquad(&F, 200, 800, 0, 1e-2, w, &result, &error, &nevals); + + rayleigh_scattering_length_d2o_weighted = result/norm; + + if (status) { + fprintf(stderr, "error integrating cerenkov distribution: %s\n", gsl_strerror(status)); + exit(1); + } + + F.function = &gsl_rayleigh_scattering_length_snoman_h2o; + + status = gsl_integration_cquad(&F, 200, 800, 0, 1e-2, w, &result, &error, &nevals); + + rayleigh_scattering_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; @@ -175,6 +247,42 @@ double get_weighted_absorption_length_snoman_h2o(void) return absorption_length_h2o_weighted; } +double get_rayleigh_scattering_length_snoman_d2o(double wavelength) +{ + /* Returns the Rayleigh scattering length in D2O in cm. */ + return 1.0/interp1d(wavelength,rayleigh_coefficient_wavelengths, rayleigh_coefficient_d2o,LEN(rayleigh_coefficient_wavelengths)); +} + +double get_rayleigh_scattering_length_snoman_h2o(double wavelength) +{ + /* Returns the Rayleigh scattering length in H2O in cm. */ + return 1.0/interp1d(wavelength,rayleigh_coefficient_wavelengths, rayleigh_coefficient_h2o,LEN(rayleigh_coefficient_wavelengths)); +} + +double get_weighted_rayleigh_scattering_length_snoman_d2o(void) +{ + /* Returns the weighted Rayleigh scattering length in D2O in cm. */ + + if (!initialized) { + fprintf(stderr, "weighted rayleigh scattering length hasn't been initialized!\n"); + exit(1); + } + + return rayleigh_scattering_length_d2o_weighted; +} + +double get_weighted_rayleigh_scattering_length_snoman_h2o(void) +{ + /* Returns the weighted Rayleigh scattering length in H2O in cm. */ + + if (!initialized) { + fprintf(stderr, "weighted rayleigh scattering length hasn't been initialized!\n"); + exit(1); + } + + return rayleigh_scattering_length_h2o_weighted; +} + double get_absorption_length_snoman_h2o(double wavelength) { /* Returns the absorption length in H2O in cm. */ |