diff options
| -rw-r--r-- | Makefile | 2 | ||||
| -rw-r--r-- | refractive_index.c | 30 | ||||
| -rw-r--r-- | refractive_index.h | 6 | ||||
| -rw-r--r-- | test.c | 73 |
4 files changed, 108 insertions, 3 deletions
@@ -5,7 +5,7 @@ calculate_limits: calculate_limits.c solid_angle.o: solid_angle.c -test: test.c solid_angle.c +test: test.c solid_angle.c refractive_index.c clean: rm calculate_limits diff --git a/refractive_index.c b/refractive_index.c new file mode 100644 index 0000000..d7d07f9 --- /dev/null +++ b/refractive_index.c @@ -0,0 +1,30 @@ +#include <math.h> +#include "refractive_index.h" + +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; + +double get_index(double p, double wavelength, double T) +{ + /* Returns the index of pure water for a given density, wavelength, and + * temperature. The density should be in units of kg/m^3, the wavelength in + * nm, and the temperature in Celsius. */ + /* normalize the density, temperature, and pressure */ + p = p/1000.0; + 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)); +} diff --git a/refractive_index.h b/refractive_index.h new file mode 100644 index 0000000..5976eca --- /dev/null +++ b/refractive_index.h @@ -0,0 +1,6 @@ +#ifndef REFRACTIVE_INDEX_H +#define REFRACTIVE_INDEX_H + +double get_index(double p, double wavelength, double T); + +#endif @@ -1,6 +1,46 @@ #include "solid_angle.h" #include <math.h> #include <stdio.h> +#include "refractive_index.h" + +/* Table of some of the tabulated values of the refractive index of water as a + * function of wavelength and temperature. In all cases, I just used the values + * for standard atmospheric pressure and assume this corresponds approximately + * to a density of 1000 kg/m^3. + * + * See Table 7 in https://aip.scitation.org/doi/pdf/10.1063/1.555859. */ + +struct refractive_index_results { + double p; + double T; + double wavelength; + double n; +} refractive_index_results[] = { + {1000.0, 0 , 361.05, 1.39468}, + {1000.0, 10, 361.05, 1.39439}, + {1000.0, 20, 361.05, 1.39353}, + {1000.0, 30, 361.05, 1.39224}, + {1000.0, 0 , 404.41, 1.34431}, + {1000.0, 10, 404.41, 1.34404}, + {1000.0, 20, 404.41, 1.34329}, + {1000.0, 30, 404.41, 1.34218}, + {1000.0, 0 , 589.00, 1.33447}, + {1000.0, 10, 589.00, 1.33422}, + {1000.0, 20, 589.00, 1.33350}, + {1000.0, 30, 589.00, 1.33243}, + {1000.0, 0 , 632.80, 1.33321}, + {1000.0, 10, 632.80, 1.33296}, + {1000.0, 20, 632.80, 1.33224}, + {1000.0, 30, 632.80, 1.33118}, + {1000.0, 0 , 1013.98, 1.32626}, + {1000.0, 10, 1013.98, 1.32604}, + {1000.0, 20, 1013.98, 1.32537}, + {1000.0, 30, 1013.98, 1.32437}, + {1000.0, 0 , 2325.42, 1.27663}, + {1000.0, 10, 2325.42, 1.27663}, + {1000.0, 20, 2325.42, 1.27627}, + {1000.0, 30, 2325.42, 1.27563}, +}; /* Table of the values of solid angle for various values of r0/r and L/r. * @@ -66,6 +106,26 @@ int isclose(double a, double b, double rel_tol, double abs_tol) return fabs(a-b) <= fmax(rel_tol*fmax(fabs(a),fabs(b)),abs_tol); } +int test_refractive_index(char *err) +{ + /* Tests the get_index() function. */ + int i; + double n; + struct refractive_index_results result; + + for (i = 0; i < sizeof(refractive_index_results)/sizeof(struct refractive_index_results); i++) { + result = refractive_index_results[i]; + n = get_index(result.p, result.wavelength, result.T); + + if (!isclose(n, result.n, 1e-4, 0)) { + sprintf(err, "n = %.5f, but expected %.5f", n, result.n); + return 1; + } + } + + return 0; +} + int test_solid_angle(char *err) { /* Tests the get_solid_angle() function. */ @@ -101,12 +161,21 @@ int test_solid_angle(char *err) int main(int argc, char **argv) { char err[256]; + int retval = 0; if (!test_solid_angle(err)) { printf("[\033[92mok\033[0m] test_solid_angle\n"); - return 0; } else { printf("[\033[91mfail\033[0m] test_solid_angle: %s\n", err); - return 1; + retval = 1; } + + if (!test_refractive_index(err)) { + printf("[\033[92mok\033[0m] test_refractive_index\n"); + } else { + printf("[\033[91mfail\033[0m] test_refractive_index: %s\n", err); + retval = 1; + } + + return retval; } |