// Graphing Program
// Author:  Lucas Finco
// Date:  Dec. 1, 1998

#include <iostream.h>
#include <graphics.h>
#include <stdlib.h>
#include <math.h>
#include <time.h>
#include <stdio.h>
#include <conio.h>
#include <dos.h>

double b=1.0, a=1.0;

void rk4(double& x, double& y, double& z, double h);
//Runge-Kutta Method--provided by J. C. Sprott

double fx(double x, double y, double z);		//x'=fx(x, y, z)
double fy(double x, double y, double z);		//y'=fy(x, y, z)
double fz(double x, double y, double z);		//z'=fz(x, y, z)

double absl(double a);				//absolute value function

int main(){
	int cont=1;
	long r=0;
	double scale=85.0;
	double x=0.001, y=0.001, z=0.001, h=0.01, x2=x+1e-8, y2=y, z2=z, t=0.0;
	double d1=1.0, df, rs, lsum=0, l1=0.0;
	double lx1, lx2, ly1, ly2, gx, gy, gz, xt, yt, gx2, gy2, gz2;

//seed randomizer
	srand(time(NULL));

//set output format
	cout.setf(ios::showpoint);
	cout.precision(10);

//set graphics
	// request autodetection
	int gdriver = DETECT, gmode, errorcode;
	// initialize graphics mode
	initgraph(&gdriver, &gmode, "");
	// read result of initialization
	errorcode = graphresult();
	if (errorcode != grOk)    // an error occurred
	{
		cout << "Graphics error: " << grapherrormsg(errorcode);
		cout << "Press any key to halt:";
		getch();
		exit(1);               // return with error code
	}

//possible loop for re-draw,
	while (cont && (r<100)){
//init
		long i=100000;
		x=-0.1422;y=-0.5594;z=-0.6496;
		x2=x+1e-8; y2=y; z2=z; t=0.0; lsum=0;
		a=-2.252;b=9.71;

//draw axes
		cleardevice();
		line(320, 0, 320, 480);
		line(0, 240, 640, 240);

//loop on system, output graph
		for (long k=1; (k<i); k++){

//call Runge Kutta Method
			rk4(x, y, z, h);
			rk4(x2,y2, z2, h);

//calculate LE
			d1=((x2-x)*(x2-x))+((y2-y)*(y2-y))+((z2-z)*(z2-z));
			if (d1>0){
				df=d1*1e16;
				rs=1.0/sqrt(df);
				x2=x+rs*(x2-x);
				y2=y+rs*(y2-y);
				z2=z+rs*(z2-z);
				if (k>1000){
					lsum=lsum+log(df);
					l1=0.5*lsum/k/h;
				}
			}

//graphical output
			gx=(z*scale);
			gy=(x*scale);
//			putpixel(t*scale, (240-gx), 2);
//			putpixel((320-gz), (240-gy), (int(x)%40));
//			putpixel((320-gz), (240-gx), (int(y)%40));
//			putpixel((320-gx), (240-gy), 4);
			putpixel((320+gx), (380+gy), 4);
			t+=h;

//look for escapers
			if ((absl(y)+absl(x)+absl(z))>500){
				break;
			}
		}
		char a;

//is chaotic?
		if (((absl(y)+absl(x)+absl(z))<500)&&(l1>0.001)){
//			cout << "T : " << t << endl;
//			cout << "X : " << x << endl;
			cout << "Highest LE : " << l1 << endl;
		}

//set loop variable 'cont'
		cout << "Again? (y/n) : ";
		cin >> a;
//a='y';
		if (a!='y' && a!='Y'){ cont=0;r++; }
		else {cont=1;r++;}
	}
	return 0;
}

void rk4(double& x, double& y, double& z, double h){
//Runge-Kutta Method
//provided by Prof. J. C. Sprott
	double k1x, k2x, k3x, k4x, k1y, k2y, k3y, k4y, k1z, k2z, k3z, k4z;
	k1x=h*fx(x, y, z);
	k1y=h*fy(x, y, z);
	k1z=h*fz(x, y, z);
	k2x=h*fx(x+(0.5*k1x), y+(0.5*k1y), z+(0.5*k1z));
	k2y=h*fy(x+(0.5*k1x), y+(0.5*k1y), z+(0.5*k1z));
	k2z=h*fz(x+(0.5*k1x), y+(0.5*k1y), z+(0.5*k1z));
	k3x=h*fx(x+(0.5*k2x), y+(0.5*k2y), z+(0.5*k2z));
	k3y=h*fy(x+(0.5*k2x), y+(0.5*k2y), z+(0.5*k2z));
	k3z=h*fz(x+(0.5*k2x), y+(0.5*k2y), z+(0.5*k2z));
	k4x=h*fx(x+k3x, y+k3y, z+k3z);
	k4y=h*fy(x+k3x, y+k3y, z+k3z);
	k4z=h*fz(x+k3x, y+k3y, z+k3z);
	x+=(k1x+(2.0*(k2x+k3x))+k4x)/6.0;
	y+=(k1y+(2.0*(k2y+k3y))+k4y)/6.0;
	z+=(k1z+(2.0*(k2z+k3z))+k4z)/6.0;
	return;
}

double fx(double x, double y, double z){
	return y+z;
}

double fy(double x, double y, double z){
	return (a*(absl(x)))-y;
}

double fz(double x, double y, double z){
	return 1.0-x;
}

double absl(double a){
	if (a>0.0){ return a; }
	else { return -a; }
}