coding Shift_JIS;

import Math;
import tool.Graph3D;

const string FILE_NAME = "lorenz.dat3d"; // 出力ファイル名
const double S = 10.0; // ローレンツ方程式のパラメータ
const double B = 4.7; // ローレンツ方程式のパラメータ
const double R = 100.7; // ローレンツ方程式のパラメータ
const double TIME_MAX = 8.0; // 終了時刻
const double DT = 1.0E-4;  // 1サイクルの時間進行
const int GRAPH_POINTS = 2400; // グラフプロット点数
const int LOOPS = TIME_MAX / DT; // 必要ループ回数

void main(){

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

	// RK4法で必要な変数
	double k1_x, k1_y, k1_z;
	double k2_x, k2_y, k2_z;
	double k3_x, k3_y, k3_z;
	double k4_x, k4_y, k4_z;

	// TSV書き込みモードでファイル「map.dat2d」を開く
	int fp = open(FILE_NAME, "wtsv");

	println("計算を開始します...");

	// RK4（4次ルンゲ=クッタ）法で座標値を計算
	for( int i=0; i<LOOPS; i++ ){

		k1_x = fX(x, y, z, S, B, R);
		k1_y = fY(x, y, z, S, B, R);
		k1_z = fZ(x, y, z, S, B, R);

		k2_x = fX(x+k1_x*DT*0.5, y+k1_y*DT*0.5, z+k1_z*DT*0.5, S, B, R);
		k2_y = fY(x+k1_x*DT*0.5, y+k1_y*DT*0.5, z+k1_z*DT*0.5, S, B, R);
		k2_z = fZ(x+k1_x*DT*0.5, y+k1_y*DT*0.5, z+k1_z*DT*0.5, S, B, R);

		k3_x = fX(x+k2_x*DT*0.5, y+k2_y*DT*0.5, z+k2_z*DT*0.5, S, B, R);
		k3_y = fY(x+k2_x*DT*0.5, y+k2_y*DT*0.5, z+k2_z*DT*0.5, S, B, R);
		k3_z = fZ(x+k2_x*DT*0.5, y+k2_y*DT*0.5, z+k2_z*DT*0.5, S, B, R);

		k4_x = fX(x+k3_x*DT, y+k3_y*DT, z+k3_z*DT, S, B, R);
		k4_y = fY(x+k3_x*DT, y+k3_y*DT, z+k3_z*DT, S, B, R);
		k4_z = fZ(x+k3_x*DT, y+k3_y*DT, z+k3_z*DT, S, B, R);

		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;

		// ファイルに座標値を書き出し
		if( i % (LOOPS/GRAPH_POINTS) == 0 ){
			writeln(fp, x, y, z);
		}
	}

	// ファイルを閉じる
	close(fp);

	println("計算が終了しました。グラフにプロットします。");

	// ファイルをグラフにプロット
	int graph = newGraph3D();
	setGraph3DFile(graph, FILE_NAME);
}


/**
 * X方向の変化率を与える関数
 */
double fX(double x, double y, double z, double s, double b, double r){
	return -s*x + s*y;
}

/**
 * Y方向の変化率を与える関数
 */
double fY(double x, double y, double z, double s, double b, double r){
	return -x*z + r*x - y;
}

/**
 * Z方向の変化率を与える関数
 */
double fZ(double x, double y, double z, double s, double b, double r){
	return x*y - b*z;
}