//  Searching Program
//  Lucas Finco
//  Date:  Dec. 1, 1998

#include <iostream.h>
#include <fstream.h>
#include <stdlib.h>
#include <math.h>
#include <time.h>

double a=1.0, b=1.0, c=1.0, d=1.0, e=1.0, f=1.0;

void rk4(double& x, double& y, double& z, double h); 	//Runge-Kutta Method
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(){
	double x=0.0, y=0.0, z=0.0, h=0.01, x2=x+1e-8, y2=y, z2=z, t=0.0;
	double d1=1.0, df, rs, lsum=0, l1, ox, oy, oz;
	int cont=1;
	long i=100000;

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

//output text for chaotic cases
	ofstream fout;
	fout.open("S2output.txt");
	if (fout.fail()){ cout << "File failed to open!"; exit(0); }

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

//loop on cases, ie, run one million systems
	for (long d=0; d<1000000; d++){

//re-seed randomizer every now and then
		if (d%1000==0){
			srand(time(NULL));
		}

//load random variables into coeff. and starting positions
		t=0.0;
		d1=1.0;
		lsum=0;
		a=10.0-(double(rand()%10000)/500.0);
		b=10.0-(double(rand()%10000)/500.0);
		c=10.0-(double(rand()%10000)/500.0);
		e=10.0-(double(rand()%10000)/500.0);
		f=10.0-(double(rand()%10000)/500.0);
		x=1.0-(double(rand()%10000)/5000.0);
		y=1.0-(double(rand()%10000)/5000.0);
		z=1.0-(double(rand()%10000)/5000.0);
		x2=x+1e-8;
		y2=y;
		z2=z;

//save starting coords.
		ox=x;
		oy=y;
		oz=z;

//keep up-to-date on status
		if (d%10000==0){
			cout << d << endl;
		}

//loop on a system
		for (long k=1; (k<i); k++){
			rk4(x, y, z, h);
			rk4(x2,y2, z2, h);
			d1=((x2-x)*(x2-x))+((y2-y)*(y2-y))+((z2-z)*(z2-z));
			if (d1>0.0){
				df=d1*1e16;
				rs=1.0/sqrt(df);
				x2=x+rs*(x2-x);
				y2=y+rs*(y2-y);
				z2=z+rs*(z2-z);
				lsum=lsum+log(df);
				l1=0.5*lsum/k/h;
			}

//test early, save time
			if (k>1000){
				if (((absl(y)+absl(x)+absl(z))>500.0)
					||((absl(y)+absl(x)+absl(z))<1e-6)
					||(l1<1e-4)){
					break;
				}
			}
			t+=h;
		}

//test for output
		if ((l1>0.01)&&((absl(y)+absl(x)+absl(z))<500.0)
				&&((absl(y)+absl(x)+absl(z))>1e-6)){
			fout << "a=" << a << ";b=" << b << ";c=" << c;
			fout << ";d=" << d << ";e=" << e << ";f=" << f;
			fout << ";\nx=" << ox << ";y=" << oy << ";z=" << oz;
			fout << ";\nHighest LE : " << l1 << endl << endl;
			cout << "___________________________________________________\n";
			cout << "\tHighest LE : " << l1 << endl;
		}
	d++;
	}
	cout << "Done";
	return 0;
}

void rk4(double& x, double& y, double& z, double h){
//Runge-Kutta Method
	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 (a*y)-(b*z);
}

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

double fz(double x, double y, double z){
	return e-(f*x);
}

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