add a function to compute the refractive index of water as a function of wavelength

4 files changed, 108 insertions, 3 deletions

diff --git a/Makefile b/Makefile
index 98f700d..bfdfb22 100644
--- a/Makefile
+++ b/Makefile
@@ -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
diff --git a/test.c b/test.c
index 3d5708d..a823711 100644
--- a/test.c
+++ b/test.c
@@ -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;
 }