#include <stdio.h>
#include <errno.h>
#include <string.h>
#include <stdlib.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_spline.h>
#include <math.h>

static int initialized = 0;

static double *x, *dEdx, *csda_range;
static size_t size;

static gsl_interp_accel *acc_dEdx;
static gsl_spline *spline_dEdx;

static gsl_interp_accel *acc_range;
static gsl_spline *spline_range;

static const double MUON_CRITICAL_ENERGY = 1.029e6;

static int init()
{
    int i, j;
    char line[256];
    char *str;
    double value;
    int n;

    FILE *f = fopen("muE_water_liquid.txt", "r");

    if (!f) {
        fprintf(stderr, "failed to open muE_water_liquid.txt: %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 10 lines since it's just a header. */
        if (i <= 10) continue;

        if (!len) continue;
        else if (line[0] == '#') continue;
        else if (strstr(line, "Minimum ionization")) continue;
        else if (strstr(line, "Muon critical energy")) continue;

        str = strtok(line," \n");

        while (str) {
            value = strtod(str, NULL);
            str = strtok(NULL," \n");
        }

        n += 1;
    }

    x = malloc(sizeof(double)*n);
    dEdx = malloc(sizeof(double)*n);
    csda_range = malloc(sizeof(double)*n);
    size = n;

    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 10 lines since it's just a header. */
        if (i <= 10) continue;

        if (!len) continue;
        else if (line[0] == '#') continue;
        else if (strstr(line, "Minimum ionization")) continue;
        else if (strstr(line, "Muon critical energy")) continue;

        str = strtok(line," \n");

        j = 0;
        while (str) {
            value = strtod(str, NULL);
            /* According to the file, the values are stored for wavelengths
             * between 230 and 700 in 1 nm increments. */
            switch (j) {
            case 0:
                x[n] = value;
                break;
            case 7:
                dEdx[n] = value;
                break;
            case 8:
                csda_range[n] = value;
                break;
            }
            j += 1;
            str = strtok(NULL," \n");
        }

        n += 1;
    }

    fclose(f);

    acc_dEdx = gsl_interp_accel_alloc();
    spline_dEdx = gsl_spline_alloc(gsl_interp_linear, size);
    gsl_spline_init(spline_dEdx, x, dEdx, size);

    acc_range = gsl_interp_accel_alloc();
    spline_range = gsl_spline_alloc(gsl_interp_linear, size);
    gsl_spline_init(spline_range, x, csda_range, size);

    initialized = 1;

    return 0;

err:
    fclose(f);

    return -1;
}

double get_range(double T, double rho)
{
    /* Returns the approximate range a muon with kinetic energy `T` will travel
     * in water before losing all of its energy. This range is interpolated
     * based on data from the PDG which uses the continuous slowing down
     * approximation.
     *
     * `T` should be in MeV, and `rho` should be in g/cm^3.
     *
     * Return value is in cm.
     *
     * See http://pdg.lbl.gov/2018/AtomicNuclearProperties/adndt.pdf. */
    if (!initialized) {
        if (init()) {
            exit(1);
        }
    }

    return gsl_spline_eval(spline_range, T, acc_range)/rho;
}

double get_E(double T, double x, double rho)
{
    /* Returns the approximate energy of a muon in water after travelling `x`
     * cm with an initial kinetic energy `T`.
     *
     * `T` should be in MeV, `x` in cm, and `rho` in g/cm^3.
     *
     * Return value is in MeV.
     *
     * See http://pdg.lbl.gov/2018/AtomicNuclearProperties/adndt.pdf. */
    double a, b, range, E;

    if (!initialized) {
        if (init()) {
            exit(1);
        }
    }

    range = gsl_spline_eval(spline_range, T, acc_range)/rho;
    /* FIXME: Need to double check if it's kosher to be using kinetic energies
     * here instead of the total energy. Equation 1 in the document uses the
     * total energy, but here I'm using the critical energy in kinetic energy,
     * so I should check to see if I need to convert both. */
    b = log(1 + T/MUON_CRITICAL_ENERGY)/range;
    a = MUON_CRITICAL_ENERGY*b;

    E = T + a*(1-exp(b*x))/b;

    if (E < 0) return 0;

    return E;
}

double get_dEdx(double T, double rho)
{
    /* Returns the approximate dE/dx for a muon in water with kinetic energy
     * `T`.
     *
     * `T` should be in MeV and `rho` in g/cm^3.
     *
     * Return value is in MeV/cm.
     *
     * See http://pdg.lbl.gov/2018/AtomicNuclearProperties/adndt.pdf. */
    if (!initialized) {
        if (init()) {
            exit(1);
        }
    }

    return gsl_spline_eval(spline_dEdx, T, acc_dEdx)/rho;
}