#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include "numcip.h"
#include "random.h"


int seed;
long lseed;
double pi,pi2;

//system variables
double *xx;//phase variable
double *w;//natural freq
long N;//total number of oscillators
double K;//coupling value
double TotalT;//total time
double del,w0;//parameters for lorenztian

//global variables for RK4
double *k1,*k2,*k3,*k4;
double *y4;
double *df;
double h;

void kuramoto(double *df,double *x);
double ranlorenz(void);
double meanfield(double *);
double norm(double *x,int dim);

void rk4(void);

main()
{
	long i;
	int n;
	double totaltime;
	double r;
	double kstep,endk;
	int M;

	FILE *outfile;

	seed=-1;
	lseed=-1;
	pi=4.0*atan(1.0);
	pi2=2.0*pi;

	//parameters
	del=1.0;//width of lorentzian
	w0=2.0;//center of lorentzian
	h=1e-1;//integration step

	printf("Kuramoto Model\n");
	printf("Number of Oscillators: [10000] [<32000]\n");
	scanf("%d",&N);
	printf("Starting coupling value for K: [0.5]\n");
	scanf("%lg",&K);
	printf("Ending coupling value for K: [5.5]\n");
	scanf("%lg",&endk);
	printf("Number of points in [Kmin,Kmax] [100]\n");
	scanf("%d",&M);
	printf("Total time for transient: [50]\n");
	scanf("%lg",&TotalT);

	kstep=(endk-K)/M;//step in K

	//memory allocations
	xx=bigvector(1,N);
	w=bigvector(1,N);
	y4=bigvector(1,N);
	df=bigvector(1,N);
	k1=bigvector(1,N);
	k2=bigvector(1,N);
	k3=bigvector(1,N);
	k4=bigvector(1,N);


	outfile=fopen("rasy.dat","w");
	fprintf(outfile,"kval\trasy\n");

	for(n=1;n<=M;n++)
	{
		//initialize
		for(i=1;i<=N;i++)
			{
			xx[i]=pi2*ran1(&seed);
			w[i]=ranlorenz();
			}

		//preiterate system to take out transient
		totaltime=0;
		do
			{
			rk4();
			totaltime+=h;
			}
		while(totaltime<TotalT);

		//estimate r_asy by evaluating its time average
		totaltime=0;
		r=0.0;
		do
			{
			rk4();
			r+=(meanfield(xx)*h);
			totaltime+=h;
			}
		while(totaltime<100);
		r/=totaltime;
		fprintf(outfile,"%lg\t%lg\n",K,r);
		fflush(outfile);

		K+=kstep;//increment K
	}

	return 0;
}

double meanfield(double *x)
{
	double cx,cy;
	long i;

	cx=cy=0.0;
	for(i=1;i<=N;i++)
		{
		cx+=cos(x[i]);
		cy+=sin(x[i]);
		}
	cx/=N;
	cy/=N;

	return sqrt(cx*cx+cy*cy);
}

double ranlorenz(void)
{
	return del*tan(pi*ran11(&lseed))+w0;
}

void kuramoto(double *df,double *x)
{
	double cx,cy;//phase in complex plane
	double r,phi;
	double coup;
	long i;

	//calculating the complex mean field
	//adding exp(itheta) in rect coord
	cx=cy=0.0;
	for(i=1;i<=N;i++)
		{
		cx+=cos(x[i]);
		cy+=sin(x[i]);
		}
	cx/=N;
	cy/=N;

	//remap x,y to r,phi[0,2pi]
	r=sqrt(cx*cx+cy*cy);
	if(cx>0)
		{
		if(cy>0)
			phi=atan(cy/cx);
		else
			phi=pi2+atan(cy/cx);
		}
	else
		{
		phi=pi+atan(cy/cx);
		}

	coup=r*K;

	for(i=1;i<=N;i++)
		df[i]=w[i]+coup*sin(phi-x[i]);

}


void rk4(void)
{
	long i;

	kuramoto(df,xx);
	for(i=1;i<=N;i++)
		{
		k1[i]=h*df[i];
		y4[i]=xx[i]+k1[i]/2.0;
		}

	kuramoto(df,y4);
	for(i=1;i<=N;i++)
		{
		k2[i]=h*df[i];
		y4[i]=xx[i]+k2[i]/2.0;
		}

	kuramoto(df,y4);
	for(i=1;i<=N;i++)
		{
		k3[i]=h*df[i];
		y4[i]=xx[i]+k3[i];
		}
	
	kuramoto(df,y4);
	for(i=1;i<=N;i++)
		k4[i]=h*df[i];

	for(i=1;i<=N;i++)
		{
		xx[i]+=(k1[i]+2*k2[i]+2*k3[i]+k4[i])/6.0;
		xx[i]=fmod(xx[i],pi2);
		}
}