coding Shift_JIS;

import Math;
import Graphics3D;
import graphics3d.Graphics3DFramework;



// グリッドのX/Y分割数
const int X_N = 64;
const int Y_N = 64;

// 座標値配列で、X/Y/Z座標を表すインデックス
const int X = 0;
const int Y = 1;
const int Z = 2;

// 頂点座標を格納する配列（最右次元は0=X,1=Y,2=Z）
float vertex[ Y_N ][ X_N ][ 3 ];

// 頂点のZ方向速度
float vertexVZ[ Y_N ][ X_N ] = 0.0;

// 頂点のZ方向の力
float vertexFZ[ Y_N ][ X_N ] = 0.0;



// 領域の辺の長さ
const float X_LENGTH = 4.0;
const float Y_LENGTH = 4.0;

// 単位面積あたりの質量密度
const float DENSITY = 10.0;

// 張力
const float TENSION = 3.0;

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

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

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

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



/**
 * 1フレームあたり、力学計算（時間発展）を行う回数です。
 *
 * 波を速くしたい場合、DTを荒くしたり、張力を上げると、
 * 精度限界を超えて運動が発散してしまう場合があります。
 * その回避策で、1フレーム内に細かいDTで繰り返し計算します。
 */
const int UPDATE_PER_FRAME = 3;

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

// 格子点間距離（ポリゴンの辺の長さ）
const float DX = X_LENGTH / (X_N-1);
const float DY = Y_LENGTH / (Y_N-1);



// 波モデルを配置する座標系
int waveCoordinate;

// モデルのIDを格納する
int waveModel;





/* --------------------------------------------------
 * ここから、初期化処理に関する部分
 * （モデルの生成や初期形状設定など）
 */


/*
 * ここはフレームワークから、
 * プログラム起動後に1度だけコールされます。
 */
void initialize(int renderer){

	// 波座標系を生成
	initializeCoordinate();

	// 波座標系を配置
	mountCoordinate(waveCoordinate, renderer);


	// 波モデルを生成して配置
	initializeModel();

	// 波モデルを波座標系に配置
	mountModel(waveModel, renderer, waveCoordinate);


	// 波モデルの初期形状をガウス波束で設定
	setGaussianWave(PULSE_WIDTH, PULSE_HEIGHT);
}


/**
 * 台座となる座標系を生成します。
 */
void initializeCoordinate(){

	// 波モデルを配置する座標系を生成
	waveCoordinate = newCoordinate();

	// 原点を 領域長/2 だけ平行移動
	moveCoordinate(
		waveCoordinate,
		-X_LENGTH/2.0, -Y_LENGTH/2.0, 0.0
	);

	// 角度を調整
	rotZCoordinate(waveCoordinate, PI/3.0);
	rotXCoordinate(waveCoordinate, -PI/4.0);
}


/**
 * モデルを生成します。
 */
void initializeModel(){

	// 四角形グリッド形式で波のモデルを生成
	waveModel = newModel(QUADRANGLE_GRID, vertex);

	// 波モデルの色設定（赤、緑、青、α値）
	setModelColor(waveModel, 0, 0, 255, 255);

	// 裏面のカリング（描画省略）を無効化
	setModelCull(waveModel, false, false);

	// 描画を面にするかワイヤーにするか、ユーザーに選択を求める
	bool fill = confirm(
		"立体的な面による描画を行いますか ?"
		+ lf()
		+ "(「いいえ」でワイヤーによる描画 )"
	);

	// モデル構成ポリゴンのフィル設定(falseでワイヤー描画)
	setModelFill(waveModel, fill);
}


/**
 * 頂点配列をガウス波束の形状に設定します。
 */
void setGaussianWave(float width, float height){
	for(int i=0; i<Y_N; i++){
		for(int j=0; j<X_N; j++){

			// 頂点のX/Y座標
			float x = j * DX;
			float y = i * DY;

			// 中心から頂点までの距離の2乗

			float r2 =
				(x - X_LENGTH/2.0)*(x - X_LENGTH/2.0)
				+ (y - Y_LENGTH/2.0)*(y - Y_LENGTH/2.0);

			// 頂点のZ座標を、中心からの距離とガウス関数で計算
			float z = height * exp( -r2 / (2*width*width) );

			// 頂点配列に座標値をセット
			vertex[i][j][X] = x;
			vertex[i][j][Y] = y;
			vertex[i][j][Z] = z;
		}
	}
}





/* --------------------------------------------------
 * ここから、更新処理に関する部分
 * （力学演算による形状の時間発展など）
 */


/*
 * ここはフレームワークから、
 * ここは1秒間に数十回、繰り返しコールされます。
 */
void update(int renderer){

	// 質点1個あたりの質量
	float m = DENSITY
		* X_LENGTH * Y_LENGTH
		/ ( (X_N+1) * (Y_N+1) );


	// UPDATE_PER_FRAME回だけ繰り返し力学演算（時間発展）
	for(int cycle=0; cycle<UPDATE_PER_FRAME; cycle++){

		// 質点にはたらく力を計算
		if(FREE_END){
			updateForceFreeEnd();
		}else{
			updateForceFixedEnd();
		}

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

				// 頂点の速度変化の計算
				vertexVZ[i][j] += vertexFZ[i][j] / m * DT;

				// 頂点の座標変化の計算
				vertex[i][j][Z] += vertexVZ[i][j] * DT;
			}
		}

		// モデルの頂点配列を更新
		setModelVertex(waveModel, QUADRANGLE_GRID, vertex);

	}
}


/**
 * 1個の隣接質点から受ける力を返します。引数 delta には DX または DY を指定します。
 * （高速化目的でマクロ関数にしているため、実行前にインライン展開されます）
 */
macro getForce(float zSide, float zCenter, float tension, float delta){
	return ( (zSide - zCenter) * tension / delta );
}


/**
 * 固定端の場合、質点にはたらく力を計算します。
 */
void updateForceFixedEnd(){

	float force;

	// 端以外の場合のみ、質点にかかる力を計算する
	for(int i=1; i<Y_N-1; i++){
		for(int j=1; j<X_N-1; j++){

			// バネからかかる力
			force = 0.0;

			// X方向に隣接する質点とのバネからかかる力
			force += getForce(vertex[i][j+1][Z], vertex[i][j][Z], TENSION, DX);
			force += getForce(vertex[i][j-1][Z], vertex[i][j][Z], TENSION, DX);

			// Y方向に隣接する質点とのバネからかかる力
			force += getForce(vertex[i+1][j][Z], vertex[i][j][Z], TENSION, DY);
			force += getForce(vertex[i-1][j][Z], vertex[i][j][Z], TENSION, DY);

			// バネからかかる力 + 減衰力 = 質点にかかる力
			vertexFZ[i][j] = force - FRICTION * vertexVZ[i][j];
		}
	}
}


/**
 * 自由端の場合、質点にはたらく力を計算します。
 */
void updateForceFreeEnd(){

	float force;

	// 質点にかかる力の計算
	for(int i=0; i<Y_N; i++){
		for(int j=0; j<X_N; j++){

			// バネからかかる力
			force = 0.0;

			// X方向に隣接する質点とのバネからかかる力
			if(j != 0){
				force += getForce(vertex[i][j-1][Z], vertex[i][j][Z], TENSION, DX);
			}
			if(j != X_N-1){
				force += getForce(vertex[i][j+1][Z], vertex[i][j][Z], TENSION, DX);
			}

			// Y方向に隣接する質点とのバネからかかる力
			if(i != 0){
				force += getForce(vertex[i-1][j][Z], vertex[i][j][Z], TENSION, DY);
			}
			if(i != Y_N-1){
				force += getForce(vertex[i+1][j][Z], vertex[i][j][Z], TENSION, DY);
			}

			// バネからかかる力 + 減衰力 = 質点にかかる力
			vertexFZ[i][j] = force - FRICTION * vertexVZ[i][j];
		}
	}
}