coding UTF-8;   // 文字コード宣言 (文字化け予防)

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

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

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

// 軌道の座標値変数
double x = 0.1;
double y = 0.1;
double z = 0.1;


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

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

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


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

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

	// グラフプロット周期ごとに座標値を書き出す
	if( i % OUTPUT_INTERVAL == 0 ){
		println(x, y, 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;
}

// 最後の座標値を書き出しておく
println(x, y, z);