/* integrands used in quadrature integration with scipy's LowLevelCallables */
#include <math.h>

const double KB = 8.617333262145179e-05;

/* FFHS functions */
double ffhsSD(double x, void *user_data) {
	double *c = (double *)user_data;

	double omega = c[0];
	double tau = c[1];

	double u = x*x;
	double res = u*u * tau;

	res /= 81. + 9.*u - 2.*u*u + u*u*u;
	res /= u*u + omega*omega * tau*tau;
	return res;
}

/* log-gaussian functions */
double logNormalDist(double tau, double tau0, double sigma) {
    return exp(- pow((log(tau/tau0) / sigma), 2) / 2.) / sqrt(2*M_PI)/sigma;
}

double logGaussian_imag_high(double u, void *user_data) {
	double *c = (double *)user_data;

	double omega = c[0];
	double tau = c[1];
	double sigma = c[2];

    double uu = exp(-u);
    double dist = logNormalDist(1./uu, tau, sigma);

    return dist * omega * uu / (pow(uu, 2) + pow(omega, 2));
}

double logGaussian_imag_low(double u, void *user_data) {
	double *c = (double *)user_data;

	double omega = c[0];
	double tau = c[1];
	double sigma = c[2];

    double uu = exp(u);

    double dist = logNormalDist(uu, tau, sigma);

    return dist * omega * uu / (1. + pow(omega*uu, 2));
}

double logGaussian_real_high(double u, void *user_data) {
	double *c = (double *)user_data;

	double omega = c[0];
	double tau = c[1];
	double sigma = c[2];

    double uu = exp(-2.*u);
    double dist = logNormalDist(exp(uu), tau, sigma);

    return dist * uu / (uu + pow(omega, 2));
}

double logGaussian_real_low(double u, void *user_data) {
	double *c = (double *)user_data;

	double omega = c[0];
	double tau = c[1];
	double sigma = c[2];

    double uu = exp(u);

    double dist = logNormalDist(uu, tau, sigma);

    return dist / (1. + pow(omega*uu, 2));
}

double logGaussianCorrelation(double x, void *user_data) {
    double *c = (double *)user_data;

    double t = c[0];
	double tau = c[1];
	double sigma = c[2];

    double uu = exp(x);

    double dist = logNormalDist(uu, tau, sigma);

    return dist * exp(-t/uu);
}

// functions for distribution of energy
double normalDist(double x, double x0, double sigma) {
    return exp(- pow((x-x0) / sigma, 2) / 2.) / sqrt(2 * M_PI) / sigma;
}

double rate(double tau0, double ea, double t) {
    return exp(-ea / t / KB) / tau0;
}

double energyDist_SD(double x, void *user_data) {
    double *c = (double *)user_data;

    double omega = c[0];
    double tau0 = c[1];
    double e_m = c[2];
    double e_b = c[3];
    double temp = c[4];

    double r = rate(tau0, x, temp);

    return r/(pow(r, 2) + pow(omega, 2)) * normalDist(x, e_m, e_b);
}

double energyDistSuscReal(double x, void *user_data) {
    double *c = (double *)user_data;

    double omega = c[0];
    double tau0 = c[1];
    double e_m = c[2];
    double e_b = c[3];
    double temp = c[4];

    double r = rate(tau0, x, temp);

    return 1 / (pow(r, 2) + pow(omega, 2)) * normalDist(x, e_m, e_b);
}

double energyDistSuscImag(double x, void *user_data) {
    double *c = (double *)user_data;

    double omega = c[0];
    double tau0 = c[1];
    double e_m = c[2];
    double e_b = c[3];
    double temp = c[4];

    double r = rate(tau0, x, temp);

    return omega * r / (pow(r, 2) + pow(omega, 2)) * normalDist(x, e_m, e_b);
}

double energyDistCorrelation(double x, void *user_data) {
    double *c = (double *)user_data;

    double t = c[0];
    double tau0 = c[1];
    double e_m = c[2];
    double e_b = c[3];
    double temp = c[4];

    double r = rate(tau0, x, temp);

    return normalDist(x, e_m, e_b) * exp(-t * r);
}

// Generalised Gamma Function
double genGammaAlphaDist(double x, void *user_data) {
    double *c = (double *)user_data;

    double omega = c[0];
    double tau0 = c[1];
    double alpha = c[2];
    double beta = c[3];

    double b_to_a = beta / alpha;
    double tau_to_tau0 = tau / tau0;

    double norm = exp(-lgamma(b_to_a) + b_to_a * log(b_to_a)) * alpha;

    return norm * exp(-b_to_a * pow(tau_to_tau0, alpha)) * pow(tau_to_tau0, beta);
}

