#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);
}