public class LorenzAttractor {

	// ローレンツ方程式のパラメータ
	static final double S = 10.0;   // σ (Sigma)
	static final double R = 100.0;  // ρ (Rho)
	static final double B = 5.0;    // β (Beta)

	// 微分方程式を数値的に解くパラメータ
	static final double DT = 1.0E-5;                 // 時間刻み (1サイクル計算ごとの時間進行)
	static final double TIME_MAX = 8.0;              // 終了時刻 (の目安、ある程度のずれあり)
	static final int LOOPS = (int)(TIME_MAX / DT);   // [自動設定] 必要な計算ループ回数

	// 値の出力用のパラメータ
	static final int OUTPUT_POINTS = 2400;   // 出力する座標点数 (の目安、ある程度のずれあり)
	static final int OUTPUT_INTERVAL = LOOPS / OUTPUT_POINTS; // [自動設定] ループ内での出力周期

	// x, y, z 座標の初期値
	static final double X_INIT = 0.1;
	static final double Y_INIT = 0.1;
	static final double Z_INIT = 0.1;


	// X方向の時間変化率を与える関数
	static double dxdt(double x, double y, double z){
		return -S*x + S*y;
	}

	// Y方向の時間変化率を与える関数
	static double dydt(double x, double y, double z){
		return -x*z + R*x - y;
	}

	// Z方向の時間変化率を与える関数
	static double dzdt(double x, double y, double z){
		return x*y - B*z;
	}


	public static void main(String[] args) {

		// 軌道の座標値変数
		double x = X_INIT;
		double y = Y_INIT;
		double z = Z_INIT;

		// ================================================================================
		// 以下、4次ルンゲ=クッタ(RK4)法で解曲線の座標値を計算
		// ================================================================================

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

			// グラフプロット周期ごとに座標値を書き出す
			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;
		}

		// 最後の座標値を書き出しておく
		System.out.println(x + "\t" + y + "\t" + z);
	}
}