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

#define N 8

#define Ca  f[0]
#define IP3 f[1]
#define z   f[2]
#define V   f[3]
#define m   f[4]
#define n   f[5]
#define h   f[6]
#define G   f[7]

double f[N];

double k[N][4];
double phi_k1[N];
double phi_k2[N];
double phi_k3[N];

double c_0 = 2;
double c_1 = 0.185;
double v_1 = 6;
double v_2 = 0.11;
double v_3 = 2.2;
double v_4 = 0.0;
double v_5 = 0.025;
double v_6 = 0.2;
double k_1 = 0.5;
double k_2 = 1;
double k_3 = 0.1;
double k_4 = 1.1;
double a_2 = 0.14;
double d_1 = 0.13;
double d_2 = 1.049;
double d_3 = 0.9434;
double d_5 = 0.082;
double alpha_G = 25;
double beta_G = 500;
double alpha = 0.8;
double tau_IP3 = 7.143;
double IP3_star = 0.16;
double alpha_Glu = 2;
//double d_Ca = 0.001;
//double d_IP3 = 0.12;

double C = 1; // muF/cm^2
double g_K  =  35; // mS/cm^2
double g_Na =  40; // mS/cm^2
double g_L  = 0.3; // mS/cm^2
double E_K  = -77; // mV
double E_Na =  55; // mV
double E_L  = -65; // mV
double I_app = 1.2;

int RandomI(int min, int max)
{
	return ((double)rand() / (RAND_MAX - 1)) * (max - min) + min;
}

double RandomD(double min, double max)
{
	return ((double)rand() / RAND_MAX) * (max - min) + min;
}

double J_channel(double f[N])
{
	return c_1 * v_1 * pow(IP3, 3) * pow(Ca, 3) * pow(z, 3) * (c_0 / c_1 - (1 + 1 / c_1) * Ca) / pow((IP3 + d_1) * (Ca + d_5), 3);
}

double J_pump(double f[N])
{
	return v_3 * pow(Ca, 2) / (pow(k_3, 2) + pow(Ca, 2));
}

double J_leak(double f[N])
{
	return c_1 * v_2 * (c_0 / c_1 - (1 + 1 / c_1) * Ca);
}

double J_in(double f[N])
{
	return v_5 + v_6 * pow(IP3, 2) / (pow(k_2, 2) + pow(IP3, 2));
}

double J_out(double f[N])
{
	return k_1 * Ca;
}

double J_PLC(double f[N])
{
	return v_4 * (Ca + (1 - alpha) * k_4) / (Ca + k_4);
}

double J_Glu(double f[N])
{
	return alpha_Glu / (1 + exp(-(G - 0.25) / 0.01));
}

double alpha_n(double f[N])
{
	return 0.02 * (V - 25) / (1 - exp(-(V - 25) / 9));
}

double beta_n(double f[N])
{
	return -0.002 * (V - 25) / (1 - exp((V - 25) / 9));
}

double alpha_m(double f[N])
{
	return 0.182 * (V + 35) / (1 - exp(-(V + 35) / 9));
}

double beta_m(double f[N])
{
	return -0.124 * (V + 35) / (1 - exp((V + 35) / 9));
}

double alpha_h(double f[N])
{
	return 0.25 * exp(-(V + 90) / 12);
}

double beta_h(double f[N])
{
	return 0.25 * exp((V + 62) / 6) / exp((V + 90) / 12);
}

double UllahJung_HodgkinHuxley(int i, double f[N])
{
	switch (i)
	{
	case 0: // Ca
		return J_channel(f) - J_pump(f) + J_leak(f) + J_in(f) - J_out(f);

	case 1: // IP3
		return (IP3_star - IP3) / tau_IP3 + J_PLC(f) + J_Glu(f);

	case 2: // z
		return a_2 * ( d_2 * ((IP3 + d_1) / (IP3 + d_3)) * (1 - z) - Ca * z );
	
	case 3: // V
		return 1000 * ((g_Na * pow(m, 3) * h * (E_Na - V) + g_K * n * (E_K - V) + g_L * (E_L - V) + I_app) / C);

	case 4: // m
		return 1000 * (alpha_m(f) * (1 - m) - beta_m(f) * m);

	case 5: // n
		return 1000 * (alpha_n(f) * (1 - n) - beta_n(f) * n);

	case 6: // h
		return 1000 * (alpha_h(f) * (1 - h) - beta_h(f) * h);
		
	case 7: // G
		return -alpha_G * G + beta_G * (1 / (1 + exp(-V / 0.5)));
	}

    return 0;
}

void RungeKutta(double dt, double f[N], double f_next[N])
{
	#pragma omp parallel for
	for (int i = 0; i < N; i++)
	{
		k[i][0] = UllahJung_HodgkinHuxley(i, f) * dt;
		phi_k1[i] = f[i] + k[i][0] / 2;
		k[i][1] = UllahJung_HodgkinHuxley(i, phi_k1) * dt;
		phi_k2[i] = f[i] + k[i][1] / 2;
		k[i][2] = UllahJung_HodgkinHuxley(i, phi_k2) * dt;
		phi_k3[i] = f[i] + k[i][2] / 2;
		k[i][3] = UllahJung_HodgkinHuxley(i, phi_k3) * dt;

		f_next[i] = f[i] + (k[i][0] + 2 * k[i][1] + 2 * k[i][2] + k[i][3]) / 6;
	}
}

void CopyArray(double source[N], double target[N])
{
	for (int i = 0; i < N; i++)
		target[i] = source[i];
}

bool Approximately(double a, double b)
{
	if (a < 0)
		a = -a;

	if (b < 0)
		b = -b;

	return a - b <= 0.000001;
}

int main() 
//main(int argc, char *argv[])
{
	//sscanf(argv[1], "%lf", &v_4);

	FILE *fp0;
	srand(time(NULL));

	double Ca_0  = 0.07;
	double IP3_0 = 0.16;
	double z_0   = 0.67;
	
	double V_0 = -58.7085;
	double m_0 = 0.0953;
	double n_0 = 0.000913;
	double h_0 = 0.3662;
	
	double G_0 = 0;
	
	//Initial values at t = 0
	f[0] = Ca_0;
	f[1] = IP3_0;
	f[2] = z_0;
	
	f[3] = V_0;
	f[4] = m_0;
	f[5] = n_0;
	f[6] = h_0;
	
	f[7] = G_0;
		
	/*fp0 = fopen("last_values.txt", "r");
	for (int i = 0; i < N; i++)
	{
		fscanf(fp0, "%lf", &f[i]);
	}
	fclose(fp0);*/

	const double t_start = 0;
	const double t_max   = 120; // sec
	const double dt      = 0.00001;

	double t = t_start;

	//fp0 = fopen("v_4.txt", "w+");
	//fprintf(fp0, "%f\n", v_4);
	//fclose(fp0);

	fp0 = fopen("results.txt", "w+");
	//setlocale(LC_NUMERIC, "French_Canada.1252");

	clock_t start_rk4, end_rk4;
	start_rk4 = clock();
	int lastPercent = -1;

	while (t < t_max || Approximately(t, t_max))
	{
		fprintf(fp0, "%f\t", t);
		for (int i = 0; i < N; i++)
		{
			fprintf(fp0, i == N - 1 ? "%f" : "%f\t", f[i]);
		}
		fprintf(fp0, "\n");

		double f_next[N];

		RungeKutta(dt, f, f_next);
		CopyArray(f_next, f);

		t += dt;

		int percent = (int)(100 * (t - t_start) / (t_max - t_start));
		if (percent != lastPercent)
		{
			printf("Progress: %d%%\n", percent);
			lastPercent = percent;
		}
	}

	end_rk4 = clock();
	double extime_rk4 = (double)(end_rk4 - start_rk4) / CLOCKS_PER_SEC;
	int minutes = (int)extime_rk4 / 60;
	int seconds = (int)extime_rk4 % 60;
	printf("\nExecution time is: %d minutes %d seconds\n ", minutes, seconds);

	fclose(fp0);

	fp0 = fopen("time_exec.txt", "w+");
	fprintf(fp0, "%f\n", extime_rk4);
	fclose(fp0);

	/*fp0 = fopen("last_values.txt", "w+");
	for (int i = 0; i < N; i++)
	{
		fprintf(fp0, i == N - 1 ? "%f" : "%f\t", f[i]);
	}
	fclose(fp0);*/
}