/*****************************************************************************************************************	
							This program was written by Kostas Chlouverakis 
							-----------------------------------------------

	It calculates the Lyapunov exponents spectrum for a DDE.
	It is done using the Euler integration, following Farmer's work and ideas:
	Farmer, J. D. "Chaotic Attractors of an Infinite-Dimensional Dynamical System." Physica D 4 (1982): 366-393.

	Variable N should be as high as possible. Try 100-400. 
	If N>100 then it will need many-many hours...
	The integration step should be equal to: h=tau/N.
	Time-step h should have usually a maximum value at 0.1-0.2. 
	So, if tau is high then N should be high as well in order the calculation to be accurate.

	Discussions with Prof. J.C. Sprott, Dr. Colet and Dr. P. Grassberger are greatly acknowledged!!!!!!!!
*******************************************************************************************************************/

#include <iostream.h>
#include <fstream.h>
#include <stdlib.h>
#include <math.h>
#include <time.h>
#include <stdio.h>
#include <cstdlib>
#include <ctime>
#define	pi 3.1415926535

#define	value_zero	0.005 
/*	If an LE is really zero, numerically it will be: LE=|value_zero| and never LE=0.00...0
	It is up to you to set it as low as you want to.
*/

#define	N 100
#define	Na (N+1)*(N+1)+N+2	/* Don't mess with this!! */ 
#define	NN (N+1)*(N+1)+N	/* Don't mess with this too!! */			

double	tau = 5. ;
double	b = 5.;
double	fi= -pi/4.;

double	imax=pow(10,5);
int		num=0;

double	X[Na], Xnew[Na] ;
double	CUM[N+1], Znorm[N+1], GSC[N+1], lyp_exp[N+1] ;

void main(void)
{
	double h=tau/N, t=0; 
	int i, k, j, l;
    
	ofstream fout;
	ofstream fout1;
	ofstream fout2;
	ofstream fout3;
	
	fout.open("xOF.txt");
	fout2.open("yOF.txt");
	fout1.open("tOF.txt");
	fout3.open("LEs.txt");

	/*  Init. Condts  */
	for (int mal=1; mal<=N; mal++)
		X[mal]=1. ;
	
	/*  Init. Condts for linear system (orthonormal frame)  */
    for (i = N+1; i <= NN; i++)
       X[i] = 0.0;
    for (i = 1; i <= N; i++)
    {
		X[(N+1)*i] = 1.0;
		CUM[i] = 0.0;
    }
	
	/* Start iterating... */
	for (int d=1; d<imax; d++)
	{	
		/* Ikeda Equations */
		Xnew[1] = X[1] + h*( b*pow(sin(X[N]+fi),2) - X[1] );
		for (k=1; k<=N; k++)
		{
			Xnew[k+1]=X[k];
			X[k]=Xnew[k];
		}

		/* Linearized Equations */	
		for (i = 0; i < N; i++)
			Xnew[(1)*N+1+i] = X[(1)*N+1+i] 
			+ h*( X[(N)*N+1+i]*2*b*sin(X[N]+fi)*cos(X[N]+fi) - X[(1)*N+1+i] );
		for (k=1; k<=N; k++)
		{
			for (i = 0; i < N; i++)
			{				
				Xnew[(k+1)*N+1+i]=X[(k)*N+1+i];
				X[(k)*N+1+i]=Xnew[(k)*N+1+i];
			}
		}

		t+=h;
		
		/*  Write output files for the orbit for i.e. num=20000 points */
		num=num+1;
		if( num<20000 )
		{
			fout<<X[1]<<endl;
			fout2<<X[N]<<endl;
			fout1<<t<<endl;
		}

		/*  Normalize first vector  */
		Znorm[1] = 0.0;
		for (j = 1; j <= N; j++)
			Znorm[1] = Znorm[1] + X[N*j+1] * X[N*j+1];
		Znorm[1] = sqrt(Znorm[1]);

		for (j = 1; j <= N; j++)
			X[N*j+1] = X[N*j+1] / Znorm[1];

		/*  Generate new orthonormal set of vectors  */
		for (j = 2; j <= N; j++)
		{
			/*  Generate J-1 GSC coefficients  */
			for (k = 1; k <= j-1; k++)
			{
				GSC[k] = 0.0;
				for (l = 1; l <= N; l++)
					GSC[k] = GSC[k] + X[N*l+j] * X[N*l+k];
			} 

			/*  Construct a new vector  */
			for (k = 1;  k <= N; k++)
				for (l = 1; l <= j-1; l++)
					X[N*k+j] = X[N*k+j] - GSC[l] * X[N*k+l];

			/*  Calculate the vector's norm  */
			Znorm[j] = 0.0;
			for (k = 1; k <= N; k++)
				Znorm[j] = Znorm[j] + X[N*k+j] * X[N*k+j];
			Znorm[j] = sqrt(Znorm[j]);

			/*  Normalize the new vector  */
			for (k = 1; k <= N; k++)
				X[N*k+j] = X[N*k+j] / Znorm[j];	
		}
	
		/*  Update running vectors magnitudes  */
		for (k = 1; k <= N; k++)
			CUM[k] = CUM[k] + log(Znorm[k]) / log(2.0);

		if (d % 10 == 0)
		{
		  for (j = 1; j <= N; j++)
               lyp_exp[j] = log(2.0)*CUM[j] / t;
		}
	  
		/*  Calculate the Lyapunov dimension Dky  */
		double sum=0, lsum, Dky, LEsum=0.;
		int	 kmax;

		for (int j1=1; j1<=N; j1++)
		{
		   sum = sum + lyp_exp[j1] ;
		   if( sum>=0 )
		   {
			   kmax=j1;
			   lsum=sum;
		   }
		}
		Dky = kmax - lsum/lyp_exp[kmax+1];

		/*  The value of the Jacobian trace, that is equal to the sum of all LEs  */
		for (int j2=1; j2<=N; j2++)
			LEsum = LEsum + lyp_exp[j2] ;

		/*  Check positive, zero and negative LEs  */
		int pos=0 ,zero=0, negative=0, jl;
		
		for (jl=1; jl<=N; jl++)
		{
		  if(lyp_exp[jl]>value_zero)
			  pos=pos+1;
		  if(fabs(lyp_exp[jl])<value_zero)
			  zero=zero+1;
		  if(lyp_exp[jl]<-value_zero)
			  negative=negative+1;
		}

		if( d%200==0)
			cout<<d/1000000.<<"\t"<<lyp_exp[1]<<"\t"<<pos<<"\t"<<zero<<"\t"
			<<negative<<"\t"<<LEsum<<"\t"<<Dky<<endl;
		  
	}
	/* End iterating! */

	/* Write output file for the LEs spectrum*/
	for (int jk=1; jk<=N; jk++)
			fout3<<lyp_exp[jk]<<endl;	
}