coding Shift_JIS;

import Math;
import tool.Graph2D;


// ラインの頂点数
const int N = 1024;

// ラインの頂点の座標
float vertexX[N];
float vertexY[N];

// ラインの頂点の速度
float vertexVY[N];

// ラインの頂点に働く力
float vertexFY[N];


// 媒質の長さ
const float LENGTH = 10.0;

// 単位長さあたりの質量密度
const float DENSITY = 100.0;

// 張力
const float TENSION = 10.0;

// 減衰力係数
const float FRICTION = 0.0;

// 初期パルスの幅
const float PULSE_WIDTH = 0.3;

// 初期パルスの高さ
const float PULSE_HEIGHT = 1.0;

// 自由端かどうか
const bool FREE_END = true;

// シミュレーションの時間刻み
const float DT = 0.01;

// 1フレームのアニメーションウェイト
const int WAIT = 20;

// グラフのIDを格納する変数
int graph;





/**
 * main関数、最初にシステムから呼び出されます。
 */
void main(){

	alert("シミュレーションを開始します");

	// 初期化処理
	initialize();

	// メインループへ
	mainLoop();

}



/**
 * メインループ部分です。
 */
void mainLoop(){

	// 描画を一定周期ごとに行うためのカウンタ
	int plotTiming = 50;
	int plotCounter = 0;

	// メインループ
	while(true){

		// 状態の更新処理（力学シミュレーション）
		update();

		// plotTiming回に一回、グラフを描画する
		if(plotCounter == plotTiming){

			// 状態をグラフに描画
			setGraph2DData(graph, vertexX, vertexY);

			// 少し待機（アニメーションウェイト）
			sleep(WAIT);

			// 描画タイミングのカウンタをリセット
			plotCounter = 0;
		}

		// 描画タイミングのカウンタを加算
		plotCounter++;
	}
}



/**
 * 初期値の設定を行います。
 */
void initialize(){

	// 頂点Xの初期値設定
	for(int i=0; i<N; i++){
		float x = i * LENGTH / N;
		vertexX[i] = x;
	}

	// 頂点Yの初期値設定（ガウス波束）
	float centerX = LENGTH / 2;
	for(int i=0; i<N; i++){

		float x = i * LENGTH / N;

		// ガウス関数
		vertexY[i] = PULSE_HEIGHT * 
			exp(
				- (x-centerX)*(x-centerX) 
				/ (2*PULSE_WIDTH*PULSE_WIDTH)
			);

		// 初期速度は 0（静止状態）
		vertexVY[i] = 0.0;
	}

	// グラフを起動
	graph = newGraph2D();

	// 範囲の自動調整機能をOFFに
	setGraph2DAutoRange(graph, false, false);

	// 範囲を設定
	setGraph2DRange(graph, 0.0, LENGTH, PULSE_HEIGHT, -PULSE_HEIGHT);

	// 初期値をプロット
	setGraph2DData(graph, vertexX, vertexY);
}



/**
 * 1ステップの時間発展を行います。
 */
void update(){

	// 質点間の距離
	float dx = LENGTH / N;

	// 質点1個あたりの質量
	float m = DENSITY * LENGTH / (N+1);

	// 端以外の頂点にかかる力の計算
	for(int i=1; i<N-1; i++){

		vertexFY[i] =
			(vertexY[i+1]+vertexY[i-1]-2.0*vertexY[i])
			* TENSION/dx - FRICTION*vertexVY[i];
	}

	// 両端にかかる力の計算
	if(FREE_END){ // 自由端

		// 隣接する1個の質点からのみ力を受ける
		vertexFY[0] = (vertexY[1]-vertexY[0])
			* TENSION/dx - FRICTION*vertexVY[0];

		vertexFY[N-1] = (vertexY[N-2]-vertexY[N-1])
			* TENSION/dx - FRICTION*vertexVY[N-1];

	}else{ //固定端

		// 壁から固定する抗力を受けるので、合力は0で動かない
		vertexFY[0] = 0.0;
		vertexFY[N-1] = 0.0;
	}

	// 位置と速度の計算
	for(int i=0; i<N; i++){

		// 頂点の速度変化の計算
		vertexVY[i] += vertexFY[i] / m * DT;

		// 頂点の座標変化の計算
		vertexY[i] += vertexVY[i] * DT;
	}
}