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

#define N 3

#define Ca  f[0]
#define IP3 f[1]
#define z   f[2]

double f[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.5;
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 d_Ca = 0.001;
//double d_IP3 = 0.12;

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 UllahJung(int i, double f[N])
{
	switch (i)
	{
	case 0: // Ca
		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) ) - ( v_3 * pow(Ca, 2) / (pow(k_3, 2) + pow(Ca, 2)) ) + ( c_1 * v_2 * ( c_0 / c_1 - (1 + 1 / c_1) * Ca ) ) + ( v_5 + v_6 * pow(IP3, 2) / (pow(k_2, 2) + pow(IP3, 2)) ) - k_1 * Ca;

	case 1: // IP3
		return ( IP3_star - IP3 ) / tau_IP3 + v_4 * (Ca + (1 - alpha) * k_4) / (Ca + k_4);

	case 2: // z
		return a_2 * ( d_2 * ((IP3 + d_1) / (IP3 + d_3)) * (1 - z) - Ca * z );
	}
  return 0;
}

void RungeKutta(double dt, double f[N], double f_next[N])
{
	double k[N][4];

	// k1
	for (int i = 0; i < N; i++)
		k[i][0] = UllahJung(i, f) * dt;

	double phi_k1[N];
	for (int i = 0; i < N; i++)
		phi_k1[i] = f[i] + k[i][0] / 2;

	// k2
	for (int i = 0; i < N; i++)
		k[i][1] = UllahJung(i, phi_k1) * dt;

	double phi_k2[N];
	for (int i = 0; i < N; i++)
		phi_k2[i] = f[i] + k[i][1] / 2;

	// k3
	for (int i = 0; i < N; i++)
		k[i][2] = UllahJung(i, phi_k2) * dt;

	double phi_k3[N];
	for (int i = 0; i < N; i++)
		phi_k3[i] = f[i] + k[i][2] / 2;

	// k4
	for (int i = 0; i < N; i++)
		k[i][3] = UllahJung(i, phi_k3) * dt;

	for (int i = 0; i < N; i++)
		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;
	
	//Initial values at t = 0
	f[0] = Ca_0;
	f[1] = IP3_0;
	f[2] = z_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   = 180; // 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_Ca_IP3_z.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);
}