#define _USE_MATH_DEFINES
#include "math.h"
#include <stdlib.h>
#include <stdio.h>
#include <locale.h>
#include <time.h>
#include <stdbool.h>

#define N 4

#define V f[0]
#define m f[1]
#define n f[2]
#define h f[3]

double f[N];

double C = 1;

double g_K  = 36;
double g_Na = 120;
double g_L  = 0.3;

double E_K  = -77;   // -77   // -12
double E_Na =  55;   //  55   // 115
double E_L  = -54.4; // -54.4 // 10

double I_app = 5.2;

const double Meander_start_from_zero = 10;
const double Meander_width = 90;
const double Meander_height = 5.1;
const double Meander_interval = 10;

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 alpha_n(double f[N])
{
	return 0.01 * (V + 55) / (1 - exp(-(V + 55) / 10));
  //return (10 - V) / (100*(exp((10-V)/10)) - 1);
}

double beta_n(double f[N])
{
	return 0.125 * exp(-(V + 65) / 80);
  //return 0.125 * exp(-V/80);
}

double alpha_m(double f[N])
{
	return 0.1 * (V + 40) / (1 - exp(-(V + 40) / 10));
  //return (25 - V) / (10*(exp((25-V)/10) - 1));
}

double beta_m(double f[N])
{
	return 4 * exp(-(V + 65) / 18);
  //return 4 * exp(-V/18);
}

double alpha_h(double f[N])
{
	return 0.07 * exp(-(V + 65) / 20);
  //return 0.07 * exp(-V/20);
}

double beta_h(double f[N])
{
	return 1 / (exp(-(V + 35) / 10) + 1);
  //return 1 / (exp((30-V)/10) + 1);
}

double I_stim(double t)
{
    if (t < Meander_start_from_zero)
        return 0;

    t -= Meander_start_from_zero;
    t = fmod(t, Meander_width + Meander_interval);

    return t < Meander_width ? Meander_height : 0;
}

double HodgkinHuxley(int i, double f[N], double t)
{
	switch (i)
	{
	case 0:
		return (g_Na * m * m * m * h * (E_Na - V) + g_K * n * n * n * n * (E_K - V) + g_L * (E_L - V) + I_stim(t)) / C;

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

	case 2:
		return alpha_n(f) * (1 - n) - beta_n(f) * n;

	case 3:
		return alpha_h(f) * (1 - h) - beta_h(f) * h;
	}
  return 0;
}

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

	// k1
	for (int i = 0; i < N; i++)
		k[i][0] = HodgkinHuxley(i, f, t) * 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] = HodgkinHuxley(i, phi_k1, t) * 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] = HodgkinHuxley(i, phi_k2, t) * 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] = HodgkinHuxley(i, phi_k3, t) * 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(int argc, char *argv[])
{
  //sscanf(argv[1], "%lf", &Meander_height);

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

	for (int i = 0; i < N; i++)
		f[i] = 0;

	const double t_start = 0;
	const double t_max = 100;
	const double dt = 0.0001;

	double t = t_start;

  fp0 = fopen("I_stim_height.txt", "a");
	fprintf(fp0, "%f\t", Meander_height);
	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);
    fprintf(fp0, "%f\t", I_stim(t));
		for (int i = 0; i < N; i++)
		{
			fprintf(fp0, i == N - 1 ? "%f" : "%f\t", f[i]);
		}
		fprintf(fp0, "\n");

		double phi_next[N];

		RungeKutta(t, dt, f, phi_next);
		CopyArray(phi_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);
}