// 軌道の座標値
	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);