/* Copyright (c) 2019, Anthony Latorre <tlatorre at uchicago>
 *
 * This program is free software: you can redistribute it and/or modify it
 * under the terms of the GNU General Public License as published by the Free
 * Software Foundation, either version 3 of the License, or (at your option)
 * any later version.

 * This program is distributed in the hope that it will be useful, but WITHOUT
 * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
 * FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for
 * more details.

 * You should have received a copy of the GNU General Public License along with
 * this program. If not, see <https://www.gnu.org/licenses/>.
 */

#include <stdio.h>
#include <gsl/gsl_spline.h>
#include <stdlib.h> /* for size_t, strtod() */
#include <errno.h> /* for errno */
#include <string.h> /* for strerror(), strtok() */
#include "sno.h"
#include "misc.h"

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

static gsl_interp_accel *acc_dEdx_rad;
static gsl_spline *spline_dEdx_rad;

static gsl_interp_accel *acc_dEdx;
static gsl_spline *spline_dEdx;

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

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

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

        if (!len) continue;
        else if (line[0] == '#') continue;

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

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

        n += 1;
    }

    x = malloc(sizeof(double)*n);
    dEdx_rad = malloc(sizeof(double)*n);
    dEdx = 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 8 lines since it's just a header. */
        if (i <= 8) continue;

        if (!len) continue;
        else if (line[0] == '#') continue;

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

        j = 0;
        while (str) {
            value = strtod(str, NULL);
            switch (j) {
            case 0:
                x[n] = value;
                break;
            case 2:
                dEdx_rad[n] = value;
                break;
            case 3:
                dEdx[n] = value;
                break;
            }
            j += 1;
            str = strtok(NULL," \n");
        }

        n += 1;
    }

    fclose(f);

    acc_dEdx_rad = gsl_interp_accel_alloc();
    spline_dEdx_rad = gsl_spline_alloc(gsl_interp_linear, size);
    gsl_spline_init(spline_dEdx_rad, x, dEdx_rad, size);

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

    return 0;

err:
    fclose(f);

    return -1;
}

/* Returns the approximate dE/dx for a electron in water with kinetic energy
 * `T`. If `T` is outside the bounds of the range table, extrapolate linearly
 * using the last two points.
 *
 * `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. */
double electron_get_dEdx(double T, double rho)
{
    if (T < spline_dEdx->x[0]) {
        return spline_dEdx->y[0];
    } else if (T < spline_dEdx->x[size-1]) {
        return gsl_spline_eval(spline_dEdx, T, acc_dEdx)*rho;
    }

    /* We extrapolate using the last two points. */
    return spline_dEdx->y[size-1] + (spline_dEdx->y[size-1]-spline_dEdx->y[size-2])*(T-spline_dEdx->x[size-1])/(spline_dEdx->x[size-1]-spline_dEdx->x[size-2]);
}

/* Returns the approximate radiative dE/dx for a electron in water with kinetic
 * energy `T`. If `T` is outside the bounds of the range table, extrapolate
 * linearly using the last two points.
 *
 * `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. */
double electron_get_dEdx_rad(double T, double rho)
{
    if (T < spline_dEdx_rad->x[0]) {
        return spline_dEdx_rad->y[0];
    } else if (T < spline_dEdx_rad->x[size-1]) {
        return gsl_spline_eval(spline_dEdx_rad, T, acc_dEdx_rad)*rho;
    }

    /* We extrapolate using the last two points. */
    return spline_dEdx_rad->y[size-1] + (spline_dEdx_rad->y[size-1]-spline_dEdx_rad->y[size-2])*(T-spline_dEdx_rad->x[size-1])/(spline_dEdx_rad->x[size-1]-spline_dEdx_rad->x[size-2]);
}

double electron_get_range(double T0)
{
    double T, dx, range;

    T = T0;

    dx = 1e-5;

    range = 0.0;
    while (T > 0) {
        double dEdx2 = electron_get_dEdx(T, WATER_DENSITY);
        T -= dEdx2*dx;
        range += dx;
    }

    return range;
}

int main(int argc, char **argv)
{
    int i;
    double T0;
    double Ts[18] = {
        2.000E+04,
        3.000E+04,
        4.000E+04,
        5.000E+04,
        6.000E+04,
        7.000E+04,
        8.000E+04,
        9.000E+04,
        1.000E+05,
        2.000E+05,
        3.000E+05,
        4.000E+05,
        5.000E+05,
        6.000E+05,
        7.000E+05,
        8.000E+05,
        9.000E+05,
        1.000E+06};

    init();

    for (i = 0; i < LEN(Ts); i++) {
        T0 = Ts[i];

        printf("%.3E - %.3E %.3E %.3E - -\n", T0, electron_get_dEdx_rad(T0,WATER_DENSITY), electron_get_dEdx(T0,WATER_DENSITY), electron_get_range(T0));
    }

    return 0;
}