move everything to src directory

1 files changed, 125 insertions, 0 deletions

diff --git a/src/quantum_efficiency.c b/src/quantum_efficiency.c
new file mode 100644
index 0000000..f8ae972
--- /dev/null
+++ b/src/quantum_efficiency.c
@@ -0,0 +1,125 @@
+#include <stdio.h>
+#include <errno.h>
+#include <string.h>
+#include <stdlib.h>
+#include <gsl/gsl_errno.h>
+#include <gsl/gsl_spline.h>
+
+static int initialized = 0;
+
+static double *x, *y;
+static size_t size;
+
+gsl_interp_accel *acc;
+gsl_spline *spline;
+
+static int init()
+{
+    int i;
+    char line[256];
+    char *str;
+    double value;
+    int n;
+
+    FILE *f = fopen("pmt_pcath_response.dat", "r");
+
+    if (!f) {
+        fprintf(stderr, "failed to open pmt_pcath_response.dat: %s", strerror(errno));
+        return -1;
+    }
+
+    i = 0;
+    n = 0;
+    /* For the first pass, we just count how many values there are. */
+    while (fgets(line, sizeof(line), f)) {
+        size_t len = strlen(line);
+        if (len && (line[len-1] != '\n')) {
+            fprintf(stderr, "got incomplete line on line %i: '%s'\n", i, line);
+            goto err;
+        }
+
+        i += 1;
+
+        /* Skip the first 32 lines since it's just a header. */
+        if (i < 32) continue;
+
+        if (!len) continue;
+        else if (line[0] == '#') continue;
+
+        str = strtok(line," \n");
+
+        while (str) {
+            value = strtod(str, NULL);
+            n += 1;
+            str = strtok(NULL," \n");
+        }
+    }
+
+    x = malloc(sizeof(double)*(n+2));
+    y = malloc(sizeof(double)*(n+2));
+    size = n + 2;
+
+    /* Make sure we extrapolate to 0. */
+    x[0] = 0;
+    y[0] = 0;
+    x[n+1] = 1000;
+    y[n+1] = 0;
+
+    i = 0;
+    n = 0;
+    /* Now, we actually store the values. */
+    rewind(f);
+    while (fgets(line, sizeof(line), f)) {
+        size_t len = strlen(line);
+        if (len && (line[len-1] != '\n')) {
+            fprintf(stderr, "got incomplete line on line %i: '%s'\n", i, line);
+            goto err;
+        }
+
+        i += 1;
+
+        /* Skip the first 32 lines since it's just a header. */
+        if (i < 32) continue;
+
+        if (!len) continue;
+        else if (line[0] == '#') continue;
+
+        str = strtok(line," \n");
+
+        while (str) {
+            value = strtod(str, NULL);
+            /* According to the file, the values are stored for wavelengths
+             * between 230 and 700 in 1 nm increments. */
+            x[n+1] = 230 + n*1.0;
+            y[n+1] = value;
+            n += 1;
+            str = strtok(NULL," \n");
+        }
+    }
+
+    fclose(f);
+
+    acc = gsl_interp_accel_alloc();
+    spline = gsl_spline_alloc(gsl_interp_linear, size);
+    gsl_spline_init(spline, x, y, size);
+
+    initialized = 1;
+
+    return 0;
+
+err:
+    fclose(f);
+
+    return -1;
+}
+
+double get_quantum_efficiency(double wavelength)
+{
+    if (!initialized) {
+        if (init()) {
+            exit(1);
+        }
+    }
+
+    return gsl_spline_eval(spline, wavelength, acc);
+}