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

#define Node_count 2*1
#define Node_wire_width (int)sqrt(Node_count)

#define Equations_per_node 8 // !!! Don't change !!!
#define Equations_count Node_count * Equations_per_node

double f[Equations_count];
double f_diff[Equations_count];

double c_0 = 2; // uM
double c_1 = 0.185;
double v_1 = 6; // s^-1
double v_2 = 0.11; // s^-1
double v_3 = 2.2; // uM/s
double v_4[Node_count]; // uM/s - Controling parameter // 0.5 //double v_4[Node_count]; // uM/s - Controling parameter //0.495
double v_5 = 0.025; // uM/s
double v_6 = 0.2;  // uM/s
double k_1 = 0.5; // s^-1
double k_2 = 1; // uM
double k_3 = 0.1;
double k_4 = 1.1; // uM/s
double a_2 = 0.14; // uM/s
double d_1 = 0.13; // uM
double d_2 = 1.049; // uM
double d_3 = 0.9434; // uM
double d_5 = 0.082; // uM
double alpha = 0.8;
double tau_IP3 = 7.143; // s
double IP3_star = 0.16; // uM
double d_Ca = 0.001; // 0.001
double d_IP3 = 0.2; // 0.12
double alpha_Glu = 0; // 2
double g_astro = 0; // 5

double C_m = 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;

double g_syn;// 0.18 // 0.04 // 0.2 // 1.6
double k_syn = 0.2; // 0.2
double E_syn[Node_count];

double alpha_G = 25; //s^-1
double beta_G = 500; //s^-1

double I_app[Node_count];

double** A_A;
double** B_A;
double*  C_A;

double** A_N;
double** B_N;
double*  C_N;

//double tau[Node_count][Node_count];

//#define tau_min  2 // ms
//#define tau_max 12 // ms

//#define ms_to_step 200 // (0.001 / dt) !!! Don't forget !!!

//#define Max_delay tau_max * ms_to_step
//double V_old_array[Node_count][Max_delay];

const double Freq = 500; // Hz
const double Min_magintude = -0.13; // pA
const double Max_magintude =  8.0; // pA
const double Duration = 0.001; // sec

double Meander_start_from_zero[Node_count];
double Meander_width[Node_count];
double Meander_height[Node_count];
double Meander_interval[Node_count];
double last_meander_end[Node_count];

double I_stim(int i, double t)
{
	if (i == 1)
		return 0;

	if (t < Meander_start_from_zero[i])
		return 0;

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

	return t < Meander_width[i] ? Meander_height[i] : 0;
}

double Ca(int i)
{
	return f[i * Equations_per_node];
}

void SetCa(int i, double value)
{
	f[i * Equations_per_node] = value;
}

double IP3(int i)
{
	return f[i * Equations_per_node + 1];
}

void SetIP3(int i, double value)
{
	f[i * Equations_per_node + 1] = value;
}

double z(int i)
{
	return f[i * Equations_per_node + 2];
}

void Setz(int i, double value)
{
	f[i * Equations_per_node + 2] = value;
}

double G(int i)
{
	return f[i * Equations_per_node + 3];
}

void SetG(int i, double value)
{
	f[i * Equations_per_node + 3] = value;
}

double V(int i)
{
	return f[i * Equations_per_node + 4];
}

void SetV(int i, double value)
{
	f[i * Equations_per_node + 4] = value;
}

double m(int i)
{
	return f[i * Equations_per_node + 5];
}

void Setm(int i, double value)
{
	f[i * Equations_per_node + 5] = value;
}

double n(int i)
{
	return f[i * Equations_per_node + 6];
}

void Setn(int i, double value)
{
	f[i * Equations_per_node + 6] = value;
}

double h(int i)
{
	return f[i * Equations_per_node + 7];
}

void Seth(int i, double value)
{
	f[i * Equations_per_node + 7] = value;
}

/*double V_old(int i, int delay)
{
	return V_old_array[i][Max_delay - 1 - delay];
}*/

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 g_syn_eff(double* f, int i)
{
	if (Ca(i) >= 0.3)
	{
		return g_syn * (1 + g_astro * Ca(i));
	}

	return g_syn;
}

double I_syn(double* f, int i)
{
	double sum_3 = 0;

	for (int j = 0; j < Node_count; j++)
	{
		sum_3 += A_N[i][j] * g_syn_eff(f, i) * (E_syn[i] - V(i)) / (1 + exp(-(V(j)) / k_syn));
	}

	return sum_3;
}

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

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

double J_leak(double* f, int i)
{
	return c_1 * v_2 * (c_0 / c_1 - (1 + 1 / c_1) * Ca(i));
}

double J_pump(double* f, int i)
{
	return v_3 * pow(Ca(i), 2) / (pow(k_3, 2) + pow(Ca(i), 2));
}

double J_in(double* f, int i)
{
	return v_5 + v_6 * pow(IP3(i), 2) / (pow(k_2, 2) + pow(IP3(i), 2));
}

double J_out(double* f, int i)
{
	return k_1 * Ca(i);
}

double J_Glu(double* f, int i)
{
  if (E_syn[i] == 0)
  {
		//printf("J_Glu = %f\n", alpha_Glu / (1 + exp(-(G(i) - 0.4) / 0.01)));
    return alpha_Glu / (1 + exp(-(G(i) - 0.22) / 0.01));
  }

  return 0;
}

double alpha_m(double* f, int i)
{
	return 0.182 * (V(i) + 35) / (1 - exp(-(V(i) + 35) / 9));
}

double beta_m(double* f, int i)
{
	return -0.124 * (V(i) + 35) / (1 - exp((V(i) + 35) / 9));
}

double alpha_n(double* f, int i)
{
	return 0.02 * (V(i) - 25) / (1 - exp(-(V(i) - 25) / 9));
}

double beta_n(double* f, int i)
{
	return -0.002 * (V(i) - 25) / (1 - exp((V(i) - 25) / 9));
}

double alpha_h(double* f, int i)
{
	return 0.25 * exp(-(V(i) + 90) / 12);
}

double beta_h(double* f, int i)
{
	return 0.25 * exp((V(i) + 62) / 6) / exp((V(i) + 90) / 12);
}

double UllahJung_HodgkinHuxley(int i, double* f, double t)
{
	int in = i / Equations_per_node;
	int il = i % Equations_per_node;

	switch (il)
	{
	  case 0: // Ca
	  {
		  double sum_1 = 0;

		  for (int j = 0; j < Node_count; j++)
		  {
		    sum_1 += d_Ca * (Ca(j) - Ca(in));
		  }

		 /*for (int j = 0; j < C_A[in]; j++)
		 {
			 sum_1 += d_Ca * (Ca((int)B_A[in][j]) - Ca(in));
		 }*/

		  return J_channel(f, in) - J_pump(f, in) + J_leak(f, in) + J_in(f, in) - J_out(f, in) + sum_1;
	  }

	case 1: // IP3
	{
		double sum_2 = 0;

		for (int j = 0; j < Node_count; j++)
		{
		  sum_2 += d_IP3 * (IP3(j) - IP3(in));
		}

		/*for (int j = 0; j < C_A[in]; j++)
		{
			sum_2 += d_IP3 * (IP3((int)B_A[in][j]) - IP3(in));
		}*/

		return (IP3_star - IP3(in)) / tau_IP3 + J_PLC(f, in) + sum_2 + J_Glu(f, in);
	}

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

	case 3: // G
	{
		return -alpha_G * G(in) + beta_G * (1 / (1 + exp(-V(in) / 0.5)));
	}

	case 4: // V
	{

		//double sum = 0;
    //for (int j = 0; j < Node_count; j++)
		//{
		//sum += A[in][j] * g_syn * (V(in) - V_old(j, tau[in][j]));
		//sum += A[in][j] * g_syn * (V(j) - V(in));
		//sum += A[in][j] * g_syn * (V(in) - E_syn[in]) / (1 + exp(-V_old(j, tau[in][j]) / k_syn));
		//sum += 1 / (0.2 * Node_count_half) * A[in][j] * g_syn * (V(in) - E_syn[in]) / (1 + exp(-V(j) / k_syn)); // i up, j down
		//	sum += /*1 / (0.2 * Node_count) */ A_N[in][j] * g_syn * (1 + g_astro * Ca(in)) * (E_syn[in] - V(in)) / (1 + exp(-(V(j)) / k_syn)); // i up, j down
		//printf("in = %d\t Node_count = %d\t A[in][j] = %f\t V(in) = %f\t E_syn[in] = %f\t sum = %f\n", in, j, A_N[in][j], V(in), E_syn[in], sum);
		//}

		/*for (int j = 0; j < C[in]; j++)
		{
		sum += sigma[in][(int)B[in][j]] * (V((int)B[in][j]) - V(in));
		}*/

		//for (int j = 0; j < C_N[in]; j++)
		//{
			//sum += g_syn * (V((int)B[in][j]) - V(in)); // устаревшая часть, нужна для проверки разностной схемы
			//sum += A[in][j] * g_syn * (V(in) - E_syn[in]) / (1 + exp(-V_old((int)B[in][j]) / k_syn)); // устаревшая часть
			//sum += A[in][(int)B[in][j]] * g_syn * (V(in) - E_syn[in]) / (1 + exp(-V_old((int)B[in][j], tau[in][(int)B[in][j]]) / k_syn));
			//sum += A[in][(int)B[in][j]] * g_syn * (V((int)B[in][j]) - E_syn[(int)B[in][j]]) / (1 + exp(-V_old(in, tau[in][(int)B[in][j]]) / k_syn));
			//sum += 1 / (0.2 * Node_count_half) * A[in][(int)B[in][j]] * g_syn * (V((int)B[in][j]) - E_syn[(int)B[in][j]]) / (1 + exp(-V(in) / k_syn)); // j up, i down
			//sum += 1 / (0.2 * Node_count) * /*(int)A_N[in][(int)B_N[in][j]] * */ g_syn * (1 + g_astro * Ca(in)) * (V(in) - E_syn[in]) / (1 + exp(-V((int)B_N[in][j]) / k_syn)); // i up, j down
			//sum += 1 / (0.2 * Node_count) * /*(int)A_N[in][(int)B_N[in][j]] * */ g_syn * (1 + g_astro * Ca(in)) * (E_syn[in] - V(in)) / (1 + exp(-V((int)B_N[in][j]) / k_syn)); // i up, j down
			//printf("i = %d\t j = %d\t A[i, j] = %d\n", in, (int)B_N[in][j], (int)A_N[in][(int)B_N[in][j]]);
      //printf("i = %d\t V_old = %f\t exp = %f\n", in, Vold, ee);
		//}
		//printf("i = %d\t sum = %f\n", in, sum);

		return 1000 * ((g_Na * pow(m(in), 3) * h(in) * (E_Na - V(in)) + g_K * n(in) * (E_K - V(in)) + g_L * (E_L - V(in)) + I_app[in] + /*I_stim(in, t) +*/ I_syn(f, in)) / C_m); // V
	}

	case 5: // m
	{
		return 1000 * (alpha_m(f, in) * (1 - m(in)) - beta_m(f, in) * m(in)); // m
	}

	case 6: // n
	{
		return 1000 * (alpha_n(f, in) * (1 - n(in)) - beta_n(f, in) * n(in)); // n
	}

	case 7: // h
	{
		return 1000 * (alpha_h(f, in) * (1 - h(in)) - beta_h(f, in) * h(in)); // h
	}
	}

	return 0;
}

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

	// k1
	for (int i = 0; i < Equations_count; i++)
		k[i][0] = UllahJung_HodgkinHuxley(i, f, t) * dt;

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

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

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

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

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

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

	for (int i = 0; i < Equations_count; 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, double* target, int 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;
}

bool CheckSameLine(int i, int j)
{
	return i / Node_wire_width == j / Node_wire_width;
}

bool IsWireNeighbors(int i, int j)
{
	if (CheckSameLine(i, j) && (i == j - 1 || i == j + 1))
		return true;

	if (i == j - Node_wire_width || i == j + Node_wire_width)
		return true;

	return false;
}

// http://preshing.com/20111007/how-to-generate-random-timings-for-a-poisson-process/
double nextTime(double rateParameter)
{
	return -log(1.0 - (double)rand() / (RAND_MAX)) / rateParameter;
}

double ext_freq; // Hz
void GenerateRandomMeander(int i, double min_start_time)
{
	double offset;

	if (ext_freq == 0)
	{
		offset = 1000;
	}
	else
	{
		offset = 1 / ext_freq; // 0.027-0.001; //nextTime(Freq); // secs
	}

	if (offset < 0)
	{
		int a = 0;
	}

	Meander_start_from_zero[i] = min_start_time + offset;
	Meander_width[i] = Duration;
	//Meander_height[i] = RandomD(Min_magintude, Max_magintude);
	Meander_height[i] = Max_magintude;
}

void FillA_N() // 1 - post synaptic neuron marker
{
	A_N[0][0] = 0;
	A_N[0][1] = 0;
	A_N[1][0] = 1;
	A_N[1][1] = 0;
}

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

  FILE *fp_g_syn;
	fp_g_syn = fopen("g_syn.txt", "w");
  fprintf(fp_g_syn, "%f\t", g_syn);
	fclose(fp_g_syn);

	sscanf(argv[2], "%lf", &I_app_1);
	FILE *fp_I_app;
	fp_I_app = fopen("I_app.txt", "w");
  fprintf(fp_I_app, "%f\t", I_app_1);
	fclose(fp_I_app);

	double g_syn_real;
  g_syn_real = /*1 / (0.2 * Node_count) */ g_syn;
  printf("g_syn_real = %f\n", g_syn_real);

	FILE *fp0;
  FILE *fp_I_stim;
	FILE *fp_I_syn;
	FILE *fp_Ca;
	FILE *fp_IP3;
	//FILE *fp_z;
	FILE *fp_G;
	FILE *fp_V;
	//FILE *fp_m;
	//FILE *fp_n;
	//FILE *fp_h;
	//FILE *fp_V_spikes;
	FILE *fp_Esyn;
	srand(time(NULL));

	//for (int i = 0; i < Node_count; i++)
	//	V_old_length[i] = 0;

	A_A = malloc(Node_count * sizeof(double));
	for (int i = 0; i < Node_count; i++)
		A_A[i] = malloc(Node_count * sizeof(double));

	B_A = malloc(Node_count * sizeof(double));
	for (int i = 0; i < Node_count; i++)
		B_A[i] = malloc(Node_count * sizeof(double));

	C_A = malloc(Node_count * sizeof(double));

	A_N = malloc(Node_count * sizeof(double));
	for (int i = 0; i < Node_count; i++)
		A_N[i] = malloc(Node_count * sizeof(double));

	B_N = malloc(Node_count * sizeof(double));
	for (int i = 0; i < Node_count; i++)
		B_N[i] = malloc(Node_count * sizeof(double));

	C_N = malloc(Node_count * sizeof(double));

	FillA_N();

	/*fp0 = fopen("tau.txt", "w+");
	for (int i = 0; i < Node_count; i++)
	{
		for (int j = 0; j < Node_count; j++)
		{
			fprintf(fp0, "%f\t", tau[i][j] / ms_to_step);
		}
		fprintf(fp0, "\n");
	}
	fclose(fp0);*/

	//setlocale(LC_NUMERIC, "French_Canada.1252");
	fp0 = fopen("test_Poisson.txt", "w+");
	for (int i = 0; i < 1000; i++)
		fprintf(fp0, "%f\n", nextTime(Freq));
	fclose(fp0);

	/*for (int i = 0; i < 6; i++)
	{
		v_4[i] = 0.6;
	}*/
	for (int i = 0; i < Node_count; i++)
	{
		v_4[i] = 0.4; // 0.4
	}

	// Initial values
	/*for (int i = 0; i < Equations_count; i++)
	{
	f[i] = 0;
	}

	for (int i = 0; i < Equations_count; i++)
	{
	I_app[i] = RandomD(9, 40);
	}*/

	for (int i = 0; i < Node_count; i++) // init array for all nodes
	{
		SetG(i, 0); // G
	}

	//double percent_stable_state = 0.4; // 0.40
	double eps_persent = 0.05; //0.05

	double Ca0  = 0.07;
	double IP30 = 0.16;
	double z0   = 0.67;

	// Initial values at t = 0
	/*for (int i = 0; i < Node_count; i++)
	{
		SetCa(i, Ca0); // Ca
		SetIP3(i, IP30); // IP3
		Setz(i, z0); // z
	}*/
	for (int i = 0; i < Node_count; i++)
	{
		SetCa(i, Ca0 + RandomD(-Ca0 * eps_persent, Ca0 * eps_persent)); // Ca
		SetIP3(i, IP30 + RandomD(-IP30 * eps_persent, IP30 * eps_persent)); // IP3
		Setz(i, z0 + RandomD(-z0 * eps_persent, z0 * eps_persent)); // z
	}

	// Stable Focus
	double V0 = -58.7085;
	double m0 = 0.0953;
	double n0 = 0.000913;
	double h0 = 0.3662;

	// Limit Cicle
	double V1 = 14.8409;
	double m1 = 0.9174;
	double n1 = 0.0140;
	double h1 = 0.0539;

	/*for (int i = 0; i < Node_count; i++) // init array for all nodes
	{
		SetV(i, V0); // V
		Setm(i, m0); // m
		Setn(i, n0); // n
		Seth(i, h0); // h
	}*/

	SetV(0, V1); // V
	Setm(0, m1); // m
	Setn(0, n1); // n
	Seth(0, h1); // h

	SetV(1, V0); // V
	Setm(1, m0); // m
	Setn(1, n0); // n
	Seth(1, h0); // h

	/*for (int i = 0; i < Node_count; i++) // init only for neurons
	{
	  double random = RandomD(0, 1);

	  SetV(i, random < percent_stable_state ? V0 + RandomD(-V0 * eps_persent, V0 * eps_persent) : V1 + RandomD(-V1 * eps_persent, V1 * eps_persent)); // V
	  Setm(i, random < percent_stable_state ? m0 + RandomD(-m0 * eps_persent, m0 * eps_persent) : m1 + RandomD(-m1 * eps_persent, m1 * eps_persent)); // m
	  Setn(i, random < percent_stable_state ? n0 + RandomD(-n0 * eps_persent, n0 * eps_persent) : n1 + RandomD(-n1 * eps_persent, n1 * eps_persent)); // n
	  Seth(i, random < percent_stable_state ? h0 + RandomD(-h0 * eps_persent, h0 * eps_persent) : h1 + RandomD(-h1 * eps_persent, h1 * eps_persent)); // h
	}*/

	/*for (int i = 0; i < Node_count; i++) // init only for neurons
	{
		SetV(i, RandomD(-80, 20)); // V
		Setm(i, RandomD(0, 1)); // m
		Setn(i, RandomD(0, 1)); // n
		Seth(i, RandomD(0, 1)); // h
	}*/

	/*for (int i = 0; i < Node_count; i++)
	{
		I_app[i] = 0.82; // init for neurons; Bifurcation point: I_app = 0.82; I_app_up = 1.04 // 0.81 or 0.93 or 1.05
	}*/

	I_app[0] = I_app_1;
	I_app[1] = 0.8;

	/*double percent_excitable = 0.8; // 0.8

	double E_syn0 =   0;
	double E_syn1 = -90;

  fp_Esyn = fopen("results_Esyn.txt", "w+");
	for (int i = 0; i < Node_count; i++)
	{
		double random = RandomD(0, 1);

		E_syn[i] = random < percent_excitable ? E_syn0 : E_syn1;
		fprintf(fp_Esyn, "%f\n", E_syn[i]);
		//printf("i = %d\t E_syn = %f\n", i, E_syn[i]);
	}
	fclose(fp_Esyn);*/

	E_syn[0] = 0; // 1st neuron
	E_syn[1] = 0; // 2nd neurom

	for (int i = 0; i < Node_count; i++)
	{
		GenerateRandomMeander(i, 0);
		last_meander_end[i] = Meander_start_from_zero[i] + Duration;
	}

	const double t_start = 0;
	const double t_max = 240.0; // 100 msec = 0.1 sec // 240
	const double dt = 0.000025; // 0.01 msec = 0.00001 sec; 0.1 msec = 0.0001 sec; 1 msec = 0.001 sec

	double t = t_start;

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

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

	//FillVOldFromCurrent();

	fp_I_stim = fopen("results_I_stim.txt", "w+");
	fp_I_syn = fopen("results_I_syn.txt", "w+");
	fp_Ca  = fopen("results_Ca.txt", "w+");
	fp_IP3 = fopen("results_IP3.txt", "w+");
	//fp_z   = fopen("results_z.txt", "w+");
	fp_G = fopen("results_G.txt", "w+");
	fp_V = fopen("results_V.txt", "w+");
	//fp_m = fopen("results_m.txt", "w+");
	//fp_n = fopen("results_n.txt", "w+");
	//fp_h = fopen("results_h.txt", "w+");
	//fp_V_spikes = fopen("results_V_spikes.txt", "w+");

	while (t < t_max || Approximately(t, t_max))
	{
		fprintf(fp_I_stim, "%f\t", t);
		fprintf(fp_I_syn, "%f\t", t);
		fprintf(fp_Ca, "%f\t", t);
		fprintf(fp_IP3, "%f\t", t);
		//fprintf(fp_z, "%f\t", t);
		fprintf(fp_G, "%f\t", t);
		fprintf(fp_V, "%f\t", t);
		//fprintf(fp_m, "%f\t", t);
		//fprintf(fp_n, "%f\t", t);
		//fprintf(fp_h, "%f\t", t);
		//fprintf(fp_V_spikes, "%f\t", t);

		for (int i = 0; i < Node_count; i++)
		{
			if (t > last_meander_end[i])
			{
				GenerateRandomMeander(i, t);
				last_meander_end[i] = Meander_start_from_zero[i] + Duration;
			}

			fprintf(fp_I_stim, "%f\t", I_stim(i, t));
		}
		fprintf(fp_I_stim, "\n");

		for (int i = 0; i < Node_count; i++)
		{
			fprintf(fp_I_syn, "%f\t", I_syn(f, i));
		}
		fprintf(fp_I_syn, "\n");

		for (int i = 0; i < Equations_count; i += Equations_per_node)
			fprintf(fp_Ca, i == Equations_count - 1 ? "%f" : "%f\t", f[i]); // Ca

		for (int i = 1; i < Equations_count; i += Equations_per_node)
			fprintf(fp_IP3, i == Equations_count - 1 ? "%f" : "%f\t", f[i]); // IP3

		//for (int i = 2; i < Equations_count; i += Equations_per_node)
		//	fprintf(fp_z, i == Equations_count - 1 ? "%f" : "%f\t", f[i]); // z

		for (int i = 3; i < Equations_count; i += Equations_per_node)
			fprintf(fp_G, i == Equations_count - 1 ? "%f" : "%f\t", f[i]); // G

		for (int i = 4; i < Equations_count; i += Equations_per_node)
			fprintf(fp_V, i == Equations_count - 1 ? "%f" : "%f\t", f[i]); // V

		//for (int i = 5; i < Equations_count; i += Equations_per_node)
		//	fprintf(fp_m, i == Equations_count - 1 ? "%f" : "%f\t", f[i]); // m

		//for (int i = 6; i < Equations_count; i += Equations_per_node)
		//	fprintf(fp_n, i == Equations_count - 1 ? "%f" : "%f\t", f[i]); // n

		//for (int i = 7; i < Equations_count; i += Equations_per_node)
		//	fprintf(fp_h, i == Equations_count - 1 ? "%f" : "%f\t", f[i]); // h

    fprintf(fp_Ca, "\n");
    fprintf(fp_IP3, "\n");
    //fprintf(fp_z, "\n");
		fprintf(fp_G, "\n");
		fprintf(fp_V, "\n");
		//fprintf(fp_m, "\n");
		//fprintf(fp_n, "\n");
		//fprintf(fp_h, "\n");

		double f_next[Equations_count];

		RungeKutta(t, dt, f, f_next);

		/*for (int i = 0; i < Equations_count; i += Equations_per_node)
		{
			double diff = f_next[i] - f[i];

			fprintf(fp_V_spikes, i == Equations_count - 1 ? "%d" : "%d\t", diff < 0 && f_diff[i] > 0 && f[i] > -10 ? 1 : 0);

			f_diff[i] = diff;
		}*/

		//fprintf(fp_V_spikes, "\n");

		CopyArray(f_next, f, Equations_count);

		t += dt;

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

		//printf("V(24) = %f\t V_old(24) = %f\n", f[24*4], V_old(24));
		//UpdateVOld();
	}

	fclose(fp_I_stim);
	fclose(fp_I_syn);
	fclose(fp_Ca);
	fclose(fp_IP3);
	//fclose(fp_z);
	fclose(fp_G);
	fclose(fp_V);
	//fclose(fp_m);
	//fclose(fp_n);
	//fclose(fp_h);
	//fclose(fp_V_spikes);

	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);

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

	for (int i = 0; i < Node_count; i++)
		free(A_A[i]);
	free(A_A);

	for (int i = 0; i < Node_count; i++)
		free(B_A[i]);
	free(B_A);

	free(C_A);

	for (int i = 0; i < Node_count; i++)
		free(A_N[i]);
	free(A_N);

	for (int i = 0; i < Node_count; i++)
		free(B_N[i]);
	free(B_N);

	free(C_N);
}