diff options
Diffstat (limited to 'src')
| -rw-r--r-- | src/optics.c | 40 |
1 files changed, 38 insertions, 2 deletions
diff --git a/src/optics.c b/src/optics.c index b01ce7d..57c9d57 100644 --- a/src/optics.c +++ b/src/optics.c @@ -5,6 +5,12 @@ #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. @@ -149,7 +155,22 @@ double get_weighted_absorption_length_snoman_h2o(void) 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])); + 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) @@ -167,7 +188,22 @@ double get_weighted_absorption_length_snoman_d2o(void) double get_absorption_length_snoman_d2o(double wavelength) { /* Returns the absorption length in D2O in cm. */ - return 1.0/interp1d(wavelength, absorption_coefficient_d2o_wavelengths, absorption_coefficient_d2o, sizeof(absorption_coefficient_d2o_wavelengths)/sizeof(absorption_coefficient_d2o_wavelengths[0])); + 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_index(double p, double wavelength, double T) |