#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>
#include "optics.h"
#include "quantum_efficiency.h"
#include "solid_angle.h"
#include "pdg.h"
#include "vector.h"
#include "muon.h"
#include "sno.h"
#include "scattering.h"
#include "pmt_response.h"
#include "misc.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;

/* Muon critical energy in H2O and D2O. These values are used in computing the
 * kinetic energy of the muon as a function of distance in get_T().
 *
 * These values come from
 * http://pdgprod.lbl.gov/~deg/AtomicNuclearProperties/HTML/deuterium_oxide_liquid.html for D2O and
 * http://pdg.lbl.gov/2018/AtomicNuclearProperties/HTML/water_liquid.html for H2O.
 */
static const double MUON_CRITICAL_ENERGY_H2O = 1029.0e6;
static const double MUON_CRITICAL_ENERGY_D2O = 967.0e3;

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

    FILE *f = fopen("muE_deuterium_oxide_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);
            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;
}

muon_energy *muon_init_energy(double T0, double rho, size_t n)
{
    /* Returns a struct with arrays of the muon position and kinetic energy.
     * This struct can then be passed to muon_get_energy() to interpolate the
     * muon's kinetic energy at any point along the track. For example:
     *
     *     muon_energy *m = muon_init_energy(1000.0, HEAVY_WATER_DENSITY, 100);
     *     double T = muon_get_energy(x, m);
     *
     * To compute the kinetic energy as a function of distance we need to solve
     * the following integral equation:
     *
     *     T(x) = T0 - \int_0^x dT(T(x))/dx
     *
     * which we solve by dividing the track up into `n` segments and then
     * numerically computing the energy at each step. */
    size_t i;
    double range, dEdx;

    /* Get the range of the muon. */
    range = get_range(T0, rho);

    muon_energy *m = malloc(sizeof(muon_energy));
    m->x = malloc(sizeof(double)*n);
    m->T = malloc(sizeof(double)*n);
    m->n = n;

    m->x[0] = 0;
    m->T[0] = T0;

    /* Now we loop over the points along the track and assume dE/dx is constant
     * between points. */
    for (i = 1; i < n; i++) {
        m->x[i] = range*i/(n-1);
        dEdx = get_dEdx(m->T[i-1], rho);
        m->T[i] = m->T[i-1] - dEdx*(m->x[i]-m->x[i-1]);
    }

    return m;
}

double muon_get_energy(double x, muon_energy *m)
{
    /* Returns the approximate kinetic energy of a muon in water after
     * travelling `x` cm with an initial kinetic energy `T`.
     *
     * Return value is in MeV. */
    double T;

    T = interp1d(x,m->x,m->T,m->n);

    if (T < 0) return 0;

    return T;
}

void muon_free_energy(muon_energy *m)
{
    free(m->x);
    free(m->T);
    free(m);
}

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

    if (T < spline_dEdx->x[0]) return spline_dEdx->y[0];

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

double get_expected_charge(double x, double T, double theta0, double *pos, double *dir, double *pmt_pos, double *pmt_normal, double r)
{
    double pmt_dir[3], cos_theta, n, wavelength0, omega, E, p, beta, z, R, f, cos_theta_pmt, absorption_length_h2o, absorption_length_d2o, l_h2o, l_d2o, distance;

    z = 1.0;

    SUB(pmt_dir,pmt_pos,pos);

    distance = NORM(pmt_dir);

    normalize(pmt_dir);

    cos_theta_pmt = DOT(pmt_dir,pmt_normal);

    if (cos_theta_pmt > 0) return 0;

    /* Calculate the cosine of the angle between the track direction and the
     * vector to the PMT. */
    cos_theta = DOT(dir,pmt_dir);

    /* Calculate total energy */
    E = T + MUON_MASS;
    p = sqrt(E*E - MUON_MASS*MUON_MASS);
    beta = p/E;

    omega = get_solid_angle_approx(pos,pmt_pos,pmt_normal,r);

    R = NORM(pos);

    /* FIXME: I just calculate delta assuming 400 nm light. */
    wavelength0 = 400.0;

    if (R <= AV_RADIUS) {
        n = get_index_snoman_d2o(wavelength0);
    } else {
        n = get_index_snoman_h2o(wavelength0);
    }

    if (beta < 1/n) return 0;

    f = get_weighted_pmt_response(acos(-cos_theta_pmt));

    absorption_length_d2o = get_weighted_absorption_length_snoman_d2o();
    absorption_length_h2o = get_weighted_absorption_length_snoman_h2o();

    get_path_length(pos,pmt_pos,AV_RADIUS,&l_d2o,&l_h2o);

    return f*exp(-l_d2o/absorption_length_d2o-l_h2o/absorption_length_h2o)*omega*FINE_STRUCTURE_CONSTANT*z*z*(1-(1/(beta*beta*n*n)))*get_probability(beta, cos_theta, theta0);
}