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

#define N 1000
#define K 5

const int w_rand_glial_min = 3;
const int w_rand_glial_max = 4;

const int w_rand_neuron_min = 5;
const int w_rand_neuron_max = 6;

const int theta_rand_min = 0;
const int theta_rand_max = 2 * M_PI;

double A[N][N];
double w[N];
double theta[N];

double kuramoto(int i, double theta[N])
{
	double sum = 0;

	for (int j = 0; j < N; j++)
		sum += A[i][j] * sin(theta[j] - theta[i]);

	return w[i] + sum * (double)K / N;
}

void RungeKutta(double dt, double theta[N], double theta_next[N])
{
	double k[N][4];

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

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

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

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

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

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

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

	for (int i = 0; i < N; i++)
		theta_next[i] = theta[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];
}

void main()
{
	srand(time(NULL));

	for (int i = 0; i < N/2; i++)
	{
		w[i] = ((double)rand() / RAND_MAX) * (w_rand_glial_max - w_rand_glial_min) + w_rand_glial_min;
	}

	for (int i = N/2 +1 ; i < N; i++)
	{
		w[i] = ((double)rand() / RAND_MAX) * (w_rand_neuron_max - w_rand_neuron_min) + w_rand_neuron_min;
	}

	/*w[0] = 5;
	w[1] = 5;
	w[2] = 5;
	w[3] = 5;
	w[4] = 5;*/

	for (int i = 0; i < N; i++)
	{
		theta[i] = ((double)rand() / RAND_MAX) * (theta_rand_max - theta_rand_min) + theta_rand_min;
	}

	/*theta[0] = 0;
	theta[1] = M_PI / 2;
	theta[2] = M_PI;
	theta[3] = ((double)3 / 2) * M_PI;
	theta[4] = 2 * M_PI;*/

	/*{
	A[0][0] = 0;
	A[0][1] = 1;
	A[0][2] = 0;
	A[0][3] = 0;
	A[0][4] = 1;

	A[1][0] = 1;
	A[1][1] = 0;
	A[1][2] = 1;
	A[1][3] = 0;
	A[1][4] = 0;

	A[2][0] = 0;
	A[2][1] = 1;
	A[2][2] = 0;
	A[2][3] = 1;
	A[2][4] = 0;

	A[3][0] = 0;
	A[3][1] = 0;
	A[3][2] = 1;
	A[3][3] = 0;
	A[3][4] = 1;

	A[4][0] = 1;
	A[4][1] = 0;
	A[4][2] = 0;
	A[4][3] = 1;
	A[4][4] = 0;
	}*/

	/*{
		A[0][0] = 1;
		A[0][1] = 1;
		A[0][2] = 1;
		A[0][3] = 1;
		A[0][4] = 1;

		A[1][0] = 1;
		A[1][1] = 1;
		A[1][2] = 1;
		A[1][3] = 1;
		A[1][4] = 1;

		A[2][0] = 1;
		A[2][1] = 1;
		A[2][2] = 1;
		A[2][3] = 1;
		A[2][4] = 1;

		A[3][0] = 1;
		A[3][1] = 1;
		A[3][2] = 1;
		A[3][3] = 1;
		A[3][4] = 1;

		A[4][0] = 1;
		A[4][1] = 1;
		A[4][2] = 1;
		A[4][3] = 1;
		A[4][4] = 1;
	}*/

	for (int i = 0; i < N; i++)
	{
		for (int j = 0; j < N; j++)
		{
			A[i][j] = 1;
		}
	}


	const double t_start = 0;
	const double t_max = 10;
	const double dt = 0.1;

	double t = t_start;

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

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

		double theta_next[N];

		RungeKutta(dt, theta, theta_next);
		CopyArray(theta_next, theta);

		t += dt;
		printf("Progress: %d%%\n", (int)(100 * t / t_max));
	}

	fclose(fp);
}