#define _USE_MATH_DEFINES
#include "math.h"
#include "malloc.h"

class collisionmodel
{
private:
	// 

public:

	// Model parameters
	double a00f, ai0f, eps0f, b0f, bMf, qf, qepsf, q2f;
	double g1a00f, g1ai0f, g1eps0f, g1b0f, g1bMf, g1qf, g1qepsf, g1q2f;
	double chi_global_e0, chi_global_bmin, chi_q1, chi_q2, chi_bM, chi_qe;
	double chi_global_coefs[6];


	void LoadModelFromFile(char *modelfile)
	{
		
		FILE *mf;
		mf=fopen(modelfile,"r");
		
		fscanf(mf,"%lf %lf %lf %lf %lf %lf %lf %lf",&a00f, &ai0f, &eps0f, &b0f, &bMf, &qf, &qepsf, &q2f);
		fscanf(mf,"%lf %lf %lf %lf %lf %lf %lf %lf",&g1a00f, &g1ai0f, &g1eps0f, &g1b0f, &g1bMf, &g1qf, &g1qepsf, &g1q2f);
		fscanf(mf,"%lf %lf %lf %lf %lf %lf",&chi_global_e0, &chi_global_bmin, &chi_q1, &chi_q2, &chi_bM, &chi_qe);
		fscanf(mf,"%lf %lf %lf %lf %lf %lf",&(chi_global_coefs[0]),&(chi_global_coefs[1]),&(chi_global_coefs[2]),&(chi_global_coefs[3]),&(chi_global_coefs[4]),&(chi_global_coefs[5]));

		fclose(mf);
	}

	void ComputeCollision(double b, double er, double ein[2], double *er_after, double *ein_after, double *chi)
	{
		double eps = er + ein[0] + ein[1];
		double gr = er/eps;
		double g1 = ein[0]/(ein[0]+ein[1]);
		double gr_after, g1_after;

		// Compute new gamma_r, gamma_1 and chi
		CollModel(eps, b, gr, g1, &gr_after, &g1_after, chi);

		// Compute new relative and internal energies
		*er_after = gr_after*eps;
		ein_after[0] = g1_after * (1-gr_after)*eps;
		ein_after[1] = (1-g1_after) * (1-gr_after)*eps;
	}

	void CollModel(double eps, double b, double g_r, double g_1, double *g_r_after, double *g_1_after, double *chi)
	{
		if (b < bMf) {
			double gr_alpha = MyAlpha(eps, b, a00f, ai0f, eps0f, b0f, bMf, qf, qepsf, q2f);
			double g1_alpha = MyAlpha(eps, b, g1a00f, g1ai0f, g1eps0f, g1b0f, g1bMf, g1qf, g1qepsf, g1q2f);
			*chi = ChiFunction(eps, b, g_r, chi_global_e0, chi_global_bmin, chi_global_coefs, chi_q1, chi_q2, chi_bM, chi_qe);

			if (RandUniform(0, 1) > gr_alpha) { *g_r_after = g_r; }
			else { *g_r_after = equilibrium_gr(); }
			if (RandUniform(0, 1) > g1_alpha) { *g_1_after = g_1; }
			else { *g_1_after = equilibrium_g1(); }
		}
		else {
			*g_r_after = g_r;
			*g_1_after = g_1;
			*chi = 0;
		}
	}

    double MyAlpha(double eps, double b, double a00, double ai0, double eps0, double b0, double bM, double q, double qeps, double q2)
    {

		double g = exp(-pow(eps / eps0, qeps));
		double f0 = 0;
		if (b <= bM) f0 = pow(1 - pow(b / bM, q2), 1.0 / q2);
		double fi = exp(-pow(b / b0, q));

		return a00 * g * f0 + ai0 * (1 - g) * fi;
	}

	double MyChi_spline(double b, double bM, double k1, double k2, double b1, double q1, double q2)
    {
        double res;
        double rb = b/bM;

        if (b <= b1) res = M_PI * (1 + k1 * pow(rb,q1) + k2 * pow(rb,q2));
        else
        {
            double sb = bM + b1;
            double db = bM - b1;
            double rb1 = b1 / bM;

            double f1 = M_PI * (1 + k1 * pow(rb1, q1) + k2 * pow(rb1, q2));
            double d1 = M_PI * (k1*q1/bM * pow(rb1, q1-1) + k2 * q2 / bM * pow(rb1, q2-1));

            double a[4] = { 0, 0, 0, 0 };
            a[3] = (d1 * db + 2 * f1) / pow(db, 3);
            a[2] = -(d1 / db + 3 * sb * a[3]) / 2.0;
            a[1] = -(f1 / db + a[3] * (bM * bM + b1 * b1 + b1 * bM) + a[2] * sb);
            a[0] = -((a[3] * bM  + a[2]) * bM + a[1]) * bM;

            res = a[0] + (a[1] + (a[2] + a[3] * b) * b) * b;
        }


        return res;
    }

	double Myb1(double er, double e0, double bM, double bmin, double qe)
    {
		return (bM - bmin) * exp(-pow(er / e0, qe)) + bmin;
    }


    double ChiFunction(double eps, double b, double gr, double e0, double bmin, double k[6], double q1, double q2, double bM, double qe)
    {
        double k1 = k[0] + k[1] * gr + k[2] * sqrt(eps);
        double k2 = k[3] + k[4] * gr + k[5] * sqrt(eps);
        double b1 = Myb1(eps * gr, e0, bM, bmin, qe);

        return MyChi_spline(b, bM, k1, k2, b1, q1, q2);
    }


	// Equilibrium distributions of collision parameters

	double equilibrium_gr()
	{
		return -cos((acos(2*RandUniform(0,1)-1)-2.0*M_PI)/3.0)+0.5;
	}

	double equilibrium_g1()
	{
		return RandUniform(0,1);
	}


	// Statistical functions

	double RandUniform(double From, double To)
	{
		return (From + (To - From) * mRandNorm());
	}

	double mRandNorm()
	{
		double res = 0;

		while (res < 0.0000000001 || res > 0.9999999999)
		{
			res = ((double)rand()) / (double)(RAND_MAX);
		}

		return res;
	}

};