public class LorenzAttractor {

	// Parameters of the Lorenz equation
	static final double S = 10.0;   // Sigma
	static final double R = 100.0;  // Rho
	static final double B = 5.0;    // Beta

	// Parameters to solve the equation numerically
	static final double DT = 1.0E-5;                 // Time-step for each cycle
	static final double TIME_MAX = 8.0;              // The rough end-time to quit the computation
	static final int LOOPS = (int)(TIME_MAX / DT);   // End count of loops

	// Parameters for outputting data
	static final int OUTPUT_POINTS = 2400;   // The rough number of points to plot
	static final int OUTPUT_INTERVAL = LOOPS / OUTPUT_POINTS; // Interval to output data in the loop 

	// Note:
	//   For keeping code simpler,
	//   the actual end-time and the actual number of points on the graph
	//   may be different with the above specified values in some extent.
	//   If you want to control them exactly, please customize code yourself.

	// Initial values of x/y/z coordinates
	static final double X_INIT = 0.1;
	static final double Y_INIT = 0.1;
	static final double Z_INIT = 0.1;

	// The function giving the time-derivative of the x coordinate: dx/dt
	static double dxdt(double x, double y, double z){
		return -S*x + S*y;
	}

	// The function giving the time-derivative of the y coordinate: dy/dt
	static double dydt(double x, double y, double z){
		return -x*z + R*x - y;
	}

	// The function giving the time-derivative of the z coordinate: dz/dt
	static double dzdt(double x, double y, double z){
		return x*y - B*z;
	}


	public static void main(String[] args) {

		// Coordinate variables
		double x = X_INIT;
		double y = Y_INIT;
		double z = Z_INIT;

		// ================================================================================
		// Solve the differential equation by using Runge-Kutta 4th order (RK4) method
		// ================================================================================

		for( int i=0; i<LOOPS; i++ ){

			// Output data for plotting graph
			if( i % OUTPUT_INTERVAL == 0 ){
				System.out.println(x + "\t" + y + "\t" + z);
			}

			double k1_x = dxdt(x, y, z);
			double k1_y = dydt(x, y, z);
			double k1_z = dzdt(x, y, z);

			double k2_x = dxdt(x + k1_x*DT/2.0, y + k1_y*DT/2.0, z + k1_z*DT/2.0);
			double k2_y = dydt(x + k1_x*DT/2.0, y + k1_y*DT/2.0, z + k1_z*DT/2.0);
			double k2_z = dzdt(x + k1_x*DT/2.0, y + k1_y*DT/2.0, z + k1_z*DT/2.0);

			double k3_x = dxdt(x + k2_x*DT/2.0, y + k2_y*DT/2.0, z + k2_z*DT/2.0);
			double k3_y = dydt(x + k2_x*DT/2.0, y + k2_y*DT/2.0, z + k2_z*DT/2.0);
			double k3_z = dzdt(x + k2_x*DT/2.0, y + k2_y*DT/2.0, z + k2_z*DT/2.0);

			double k4_x = dxdt(x + k3_x*DT, y + k3_y*DT, z + k3_z*DT);
			double k4_y = dydt(x + k3_x*DT, y + k3_y*DT, z + k3_z*DT);
			double k4_z = dzdt(x + k3_x*DT, y + k3_y*DT, z + k3_z*DT);

			x += (k1_x + 2.0*k2_x + 2.0*k3_x + k4_x) * DT / 6.0;
			y += (k1_y + 2.0*k2_y + 2.0*k3_y + k4_y) * DT / 6.0;
			z += (k1_z + 2.0*k2_z + 2.0*k3_z + k4_z) * DT / 6.0;
		}

		// Output last data
		System.out.println(x + "\t" + y + "\t" + z);
	}
}