#include <iostream>


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

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

// 値の出力用のパラメータ
const long long int OUTPUT_POINTS = 2400LL;   // 出力する座標点数 (の目安、ある程度のずれあり)
const long long 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)法で解曲線の座標値を計算
// ================================================================================

int main() {

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

		// グラフプロット周期ごとに座標値を書き出す
		if( i % OUTPUT_INTERVAL == 0 ){
			std::cout << x << "\t" << y << "\t" << z << std::endl;
		}

		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;
	}

	// 最後の座標値を書き出しておく
	std::cout << x << "\t" << y << "\t" << z << std::endl;

	return 0;
}