/* 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://ww
/* 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 "scattering.h"
#include "quantum_efficiency.h"
#include <gsl/gsl_sf_bessel.h>
#include <gsl/gsl_integration.h>
#include "optics.h"
#include "sno.h"
#include <stddef.h> /* for size_t */
#include <stdlib.h> /* for exit() */
#include <stdio.h> /* for fprintf() */
#include <gsl/gsl_errno.h> /* for gsl_strerror() */
#include "pdg.h"
#include "misc.h"

static int initialized = 0;

static double xlo = -1.0;
static double xhi = 1.0;
static size_t nx = 1000;
static double ylo = 0.0;
static double yhi = MAX_THETA0;
static size_t ny = 1000;

static double *x;
static double *y;
static double *z;

static double x2lo = 0.0;
static double x2hi = 1.0;
static size_t nx2 = 1000;

static double *x2;
static double *y2;

/* Quantum efficiency weighted by the Cerenkov spectrum. */
static double weighted_qe;

static double prob_scatter(double wavelength, void *params)
{
    /* Calculate the number of photons emitted per unit solid angle per cm at
     * an angle `theta` for a particle travelling with velocity `beta` with an
     * angular distribution of width `theta0`. */
    double qe, delta, index;

    double beta_cos_theta = ((double *) params)[0];
    double beta_sin_theta_theta0 = ((double *) params)[1];

    qe = get_quantum_efficiency(wavelength);

    index = get_index_snoman_d2o(wavelength);

    delta = (1.0/index - beta_cos_theta)/beta_sin_theta_theta0;

    return FINE_STRUCTURE_CONSTANT*qe*exp(-pow(delta,2)/2.0)/pow(wavelength,2)*1e7/sqrt(2*M_PI)/beta_sin_theta_theta0;
}

static double prob_scatter2(double wavelength, void *params)
{
    /* Calculate the number of photons emitted per per cm for a particle
     * travelling with velocity `beta`. */
    double qe, index;

    double beta = ((double *) params)[0];

    qe = get_quantum_efficiency(wavelength);

    index = get_index_snoman_d2o(wavelength);

    return 2*M_PI*FINE_STRUCTURE_CONSTANT*(1-(1/(beta*beta*index*index)))*qe/pow(wavelength,2)*1e7;
}

double get_weighted_quantum_efficiency(void)
{
    /* Returns the quantum efficiency weighted by the Cerenkov wavelength
     * distribution, i.e. the probability that a photon randomly sampled from
     * the Cerenkov wavelenght distribution between 200 and 800 nm is detected. */
    if (!initialized) {
        fprintf(stderr, "quantum efficiency hasn't been initialized!\n");
        exit(1);
    }

    return weighted_qe;
}

static double gsl_quantum_efficiency(double wavelength, void *params)
{
    /* Returns the quantum efficiency times the Cerenkov wavelength
     * distribution. */
    double qe;

    qe = get_quantum_efficiency(wavelength);

    return qe/pow(wavelength,2);
}

void init_interpolation(void)
{
    size_t i, j;
    double params[2];
    double result, error;
    size_t nevals;
    int status;
    gsl_integration_cquad_workspace *w;
    gsl_function F;

    x = malloc(nx*sizeof(double));
    y = malloc(ny*sizeof(double));
    z = malloc(nx*ny*sizeof(double));

    for (i = 0; i < nx; i++) {
        x[i] = xlo + (xhi-xlo)*i/(nx-1);
    }

    for (i = 0; i < ny; i++) {
        y[i] = ylo + (yhi-ylo)*i/(ny-1);
    }

    w = gsl_integration_cquad_workspace_alloc(100);

    F.function = &prob_scatter;
    F.params = params;

    for (i = 0; i < nx; i++) {
        for (j = 0; j < ny; j++) {
            params[0] = x[i];
            params[1] = y[j];

            status = gsl_integration_cquad(&F, 200, 800, 0, 1e-2, w, &result, &error, &nevals);

            if (status) {
                fprintf(stderr, "error integrating photon angular distribution: %s\n", gsl_strerror(status));
                exit(1);
            }

            z[i*ny + j] = result;
        }
    }

    x2 = malloc(nx2*sizeof(double));
    y2 = malloc(nx2*sizeof(double));

    for (i = 0; i < nx2; i++) {
        x2[i] = x2lo + (x2hi-x2lo)*i/(nx2-1);
    }

    F.function = &prob_scatter2;
    F.params = params;

    for (i = 0; i < nx2; i++) {
        params[0] = x2[i];

        status = gsl_integration_cquad(&F, 200, 800, 0, 1e-2, w, &result, &error, &nevals);

        if (status) {
            fprintf(stderr, "error integrating photon angular distribution: %s\n", gsl_strerror(status));
            exit(1);
        }

        y2[i] = result;
    }

    F.function = &gsl_quantum_efficiency;
    F.params = params;

    status = gsl_integration_cquad(&F, 200, 800, 0, 1e-2, w, &result, &error, &nevals);

    if (status) {
        fprintf(stderr, "error integrating quantum efficiency distribution: %s\n", gsl_strerror(status));
        exit(1);
    }

    weighted_qe = result*800/3.0;

    initialized = 1;

    gsl_integration_cquad_workspace_free(w);
}

double get_probability(double beta, double cos_theta, double sin_theta, double theta0)
{
    /* Make sure theta0 is less than MAX_THETA0, otherwise it's possible that
     * interp2d() will quit. */
    if (theta0 > MAX_THETA0) theta0 = MAX_THETA0;

    return beta*interp2d(beta*cos_theta, beta*sin_theta*theta0, x, y, z, nx, ny);
}

double get_probability2(double beta)
{
    return interp1d(beta, x2, y2, nx2);
}

void free_interpolation(void)
{
    free(x);
    free(y);
    free(z);

    free(x2);
    free(y2);
}