波の屈折のシミュレーション

coding UTF-8;

import Math;
import Graphics;
import Graphics3D;
import graphics3d.Graphics3DFramework;


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


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


// シミュレーション終了時刻（単位はフレーム数）
const int END_TIME_COUNT = 500;

// 頂点座標を格納する配列（最右次元は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;

// 頂点の質量
float mass[ Y_N ][ X_N ];

// 面密度分布の画像ファイルのパス
const string DENSITY_FILE_PATH = "density.png";

// 面密度のスケール値
const float DENSITY_SCALE = 10.0;


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

// 張力（単位幅のリボン形状を引っ張った際の張力）
const float TENSION = 3.0;

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

// パルス波の幅
const float PULSE_WIDTH = 0.1;

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

// サイン波の角速度
const float SIN_SPEED = 20.0;

// サイン波の高さ
const float SIN_HEIGHT = 0.2;

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

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

// シミュレーションの時間刻み
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;

// 波のタイプ（実行時に選択）
const int WAVE_TYPE_SIN = 1; // 正弦波
const int WAVE_TYPE_PULSE = 2; // パルス波（孤立波）
int waveType = WAVE_TYPE_SIN;



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


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

	// 波のタイプを選択
	selectWaveType();

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

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


	// 波座標系に座標軸モデルを生成して配置
	mountModel(
		newAxisModel(3.0,3.0,3.0),
		renderer, waveCoordinate
	);


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

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

	// 波モデルの初期形状を設定
	if(waveType == WAVE_TYPE_PULSE){
		setGaussianWave(PULSE_WIDTH, PULSE_HEIGHT);
	}else{
		setFlatWave();
	}

	// 面密度分布ファイルを読み込む
	loadDensityFile(DENSITY_FILE_PATH);
}


/**
 * 波のタイプを選択するダイアログを提示します。
 */
void selectWaveType(){
	alert("波のタイプを選択して下さい：");

	string sin = "Sin Wave - 正弦波（サイン波）";
	string pulse = "Pulse Wave - 孤立波（パルス波）";
	string options[] = {sin, pulse};

	string selected = select(options);

	if(selected == sin){
		waveType = WAVE_TYPE_SIN;
	}else if(selected == pulse){
		waveType = WAVE_TYPE_PULSE;
	}else{
		error("波のタイプが不正です：" + selected);
		exit();
	}
}


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

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

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


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

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

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

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


/**
 * 頂点配列をガウス波束の形状に設定します（パルス波の場合用）。
 */
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;

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

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


/**
 * 頂点配列をフラットな形状に設定します（サイン波の場合用）。
 */
void setFlatWave(){
	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;

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


/**
 * 密度分布が描かれた画像ファイルを読み込み、設定します。
 */
void loadDensityFile(string filePath){

	// 画像ファイルからグラフィックスデータを生成
	int graphics = newGraphics(filePath);

	// ピクセルの色成分を取得（ [ Y ][ X ][ R/G/B ] ）
	int pixel[][][] = getGraphicsPixel(graphics);

	// グラフィックスデータを破棄
	deleteGraphics(graphics);

	// 色配列の要素数を確認
	int width = length(pixel, 1);
	int height = length(pixel, 0);
	if(width!=X_N || height!=Y_N){
		alert("密度分布のX/Yサイズと、メッシュのX/Yサイズが一致していません。");
		exit();
	}

	// 面密度スケールを頂点質量のスケールに変換
	float massScale = DENSITY_SCALE
		* X_LENGTH * Y_LENGTH
		/ ( X_N * Y_N );

	// グローバルの頂点質量配列 mass に値を設定
	mass = getMassFromPixel(pixel, massScale);

	// 波モデルに色を塗る
	setModelColorFromPixel(waveModel, pixel);
}


/**
 * ピクセルの色配列の明るさを密度分布と見なし、
 * 頂点質量の配列に変換してを返します。
 *
 * @param pixel ピクセルの色配列（ [ X ][ Y ][ R/G/B ] ）
 * @param densityScale 質量のスケール値（係数）
 */
float[][] getMassFromPixel(int pixel[][][], float massScale){

	const int R = 0; // 赤色のインデックス
	const int G = 1; // 青色のインデックス
	const int B = 2; // 緑色のインデックス

	// ピクセルの明るさから頂点質量を計算、設定
	for(int i=0; i<Y_N; i++){
		for(int j=0; j<X_N; j++){

			// 2D画像内の、3D頂点位置に対応する座標
			int x2d = j;
			int y2d = Y_N-1-i; //2D画像は左上頂点なのでYを反転

			// 明るさ
			float brightness = (
				pixel[ y2d ][ x2d ][ R ] / 255.0
				+
				pixel[ y2d ][ x2d ][ G ] / 255.0
				+
				pixel[ y2d ][ x2d ][ B ] / 255.0
			) / 3.0;

			// 明るいほうが軽い、暗いほうが重い
			mass[i][j] = (1.0-brightness) * massScale;
		}
	}

	return mass;
}


/**
 * 2D画像のピクセル色配列から、
 * 格子メッシュ形式のモデルに色を塗ります。
 *
 * ピクセル色配列は頂点色と対応しますが、
 * 現状はポリゴンをベタ塗り（フラットシェーディング）するので、
 * 画像右端と下端のラインは使用されません。
 *
 * ピクセル色配列は、左上頂点である事を前提としています。
 *
 * @param model 四角形格子メッシュ形式(QUADRANGLE_GRID)のモデル
 * @param pixel 2D画像のピクセル色配列[ Y ][ X ][ R/G/B/A ]
 */
void setModelColorFromPixel(int model, int pixel[][][]){

	for(int i=0; i<Y_N-1; i++){
		for(int j=0; j<X_N-1; j++){

			// モデル内のポリゴンのインデックス
			int index = i*(X_N-1) + j;

			// 2D画像内の、ポリゴン位置に対応する座標
			int x2d = j;
			int y2d = Y_N-1-i; //2D画像は左上頂点なのでYを反転

			// ポリゴンに対応する位置にあるピクセルの色を取得
			int polygonColor[4];
			polygonColor[0] = pixel[ y2d ][ x2d ][0]; //RED
			polygonColor[1] = pixel[ y2d ][ x2d ][1]; //GREEN
			polygonColor[2] = pixel[ y2d ][ x2d ][2]; //BLUE
			polygonColor[3] = pixel[ y2d ][ x2d ][3]; //ALPHA

			// ポリゴンの色を設定
			setModelPolygonColor(model, index, polygonColor);
		}
	}
}





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


/** 時刻変数 */
int timeCount = 0;


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

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

		timeCount++;
		if(timeCount == END_TIME_COUNT){
			alert("終了時刻です。より長い時間のシミュレーションを行うには、END_TIME_COUNT の値を大きく設定してください。");
		} else if(timeCount > END_TIME_COUNT){
			// 終了時刻以降は何しない。
			return;
		}

		// X / Y方向の線張力（※TENSIONは単位幅あたりの張力）
		float tensionX = TENSION * Y_LENGTH / Y_N;
		float tensionY = TENSION * X_LENGTH / X_N;

		// 波のタイプがサイン波の場合は、境界をsin関数で振動させる
		if(waveType == WAVE_TYPE_SIN){
			for(int i=0; i<Y_N; i++){
				vertex[i][0][Z] = SIN_HEIGHT * sin(timeCount * DT * SIN_SPEED);
			}
		}

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

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

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

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

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

	}
}


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


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

	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], tensionX, DX);
			force += getForce(vertex[i][j-1][Z], vertex[i][j][Z], tensionX, DX);

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

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


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

	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], tensionX, DX);
			}
			if(j != X_N-1){
				force += getForce(vertex[i][j+1][Z], vertex[i][j][Z], tensionX, DX);
			}

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

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

...


// 頂点の質量
float mass[ Y_N ][ X_N ];

// 面密度分布の画像ファイルのパス
const string DENSITY_FILE_PATH = "density.png";

// 面密度のスケール値
const float DENSITY_SCALE = 10.0;


...

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

	// 波のタイプを選択
	selectWaveType();

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

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


	// 波座標系に座標軸モデルを生成して配置
	mountModel(
		newAxisModel(3.0,3.0,3.0),
		renderer, waveCoordinate
	);


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

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

	// 波モデルの初期形状を設定
	if(waveType == WAVE_TYPE_PULSE){
		setGaussianWave(PULSE_WIDTH, PULSE_HEIGHT);
	}else{
		setFlatWave();
	}

	// 面密度分布ファイルを読み込む
	loadDensityFile(DENSITY_FILE_PATH);
}

/**
 * 密度分布が描かれた画像ファイルを読み込み、設定します。
 */
void loadDensityFile(string filePath){

	// 画像ファイルからグラフィックスデータを生成
	int graphics = newGraphics(filePath);

	// ピクセルの色成分を取得（ [ Y ][ X ][ R/G/B ] ）
	int pixel[][][] = getGraphicsPixel(graphics);

	// グラフィックスデータを破棄
	deleteGraphics(graphics);

	// 色配列の要素数を確認
	int width = length(pixel, 1);
	int height = length(pixel, 0);
	if(width!=X_N || height!=Y_N){
		alert("密度分布のX/Yサイズと、メッシュのX/Yサイズが一致していません。");
		exit();
	}

	// 面密度スケールを頂点質量のスケールに変換
	float massScale = DENSITY_SCALE
		* X_LENGTH * Y_LENGTH
		/ ( X_N * Y_N );

	// グローバルの頂点質量配列 mass に値を設定
	mass = getMassFromPixel(pixel, massScale);

	// 波モデルに色を塗る
	setModelColorFromPixel(waveModel, pixel);
}

/**
 * ピクセルの色配列の明るさを密度分布と見なし、
 * 頂点質量の配列に変換してを返します。
 *
 * @param pixel ピクセルの色配列（ [ X ][ Y ][ R/G/B ] ）
 * @param densityScale 質量のスケール値（係数）
 */
float[][] getMassFromPixel(int pixel[][][], float massScale){

	const int R = 0; // 赤色のインデックス
	const int G = 1; // 青色のインデックス
	const int B = 2; // 緑色のインデックス

	// ピクセルの明るさから頂点質量を計算、設定
	for(int i=0; i<Y_N; i++){
		for(int j=0; j<X_N; j++){

			// 2D画像内の、3D頂点位置に対応する座標
			int x2d = j;
			int y2d = Y_N-1-i; //2D画像は左上原点なのでYを反転

			// 明るさ
			float brightness = (
				pixel[ y2d ][ x2d ][ R ] / 255.0
				+
				pixel[ y2d ][ x2d ][ G ] / 255.0
				+
				pixel[ y2d ][ x2d ][ B ] / 255.0
			) / 3.0;

			// 明るいほうが軽い、暗いほうが重い
			mass[i][j] = (1.0-brightness) * massScale;
		}
	}

	return mass;
}

/**
 * 2D画像のピクセル色配列から、
 * 格子メッシュ形式のモデルに色を塗ります。
 *
 * ピクセル色配列は頂点色と対応しますが、
 * 現状はポリゴンをベタ塗り（フラットシェーディング）するので、
 * 画像右端と下端のラインは使用されません。
 *
 * ピクセル色配列は、左上頂点である事を前提としています。
 *
 * @param model 四角形格子メッシュ形式(QUADRANGLE_GRID)のモデル
 * @param pixel 2D画像のピクセル色配列[ Y ][ X ][ R/G/B/A ]
 */
void setModelColorFromPixel(int model, int pixel[ ][ ][ ]){

	for(int i=0; i < Y_N-1; i++){
		for(int j=0; j < X_N-1; j++){

			// モデル内のポリゴンのインデックス
			int index = i*(X_N-1) + j;

			// 2D画像内の、ポリゴン位置に対応する座標
			int x2d = j;
			int y2d = Y_N-1-i; //2D画像は左上原点なのでYを反転

			// ポリゴンに対応する位置にあるピクセルの色を取得
			int polygonColor[4];
			polygonColor[0] = pixel[ y2d ][ x2d ][0]; //RED
			polygonColor[1] = pixel[ y2d ][ x2d ][1]; //GREEN
			polygonColor[2] = pixel[ y2d ][ x2d ][2]; //BLUE
			polygonColor[3] = pixel[ y2d ][ x2d ][3]; //ALPHA

			// ポリゴンの色を設定
			setModelPolygonColor(model, index, polygonColor);
		}
	}
}

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

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

		timeCount++;
		if(timeCount == END_TIME_COUNT){
			alert("終了時刻です。より長い時間のシミュレーションを行うには、END_TIME_COUNT の値を大きく設定してください。");
		} else if(timeCount > END_TIME_COUNT){
			// 終了時刻以降は何しない。
			return;
		}

		// X / Y方向の線張力（※TENSIONは単位幅あたりの張力）
		float tensionX = TENSION * Y_LENGTH / Y_N;
		float tensionY = TENSION * X_LENGTH / X_N;

		// 波のタイプがサイン波の場合は、境界をsin関数で振動させる
		if(waveType == WAVE_TYPE_SIN){
			for(int i=0; i<Y_N; i++){
				vertex[i][0][Z] = SIN_HEIGHT * sin(timeCount * DT * SIN_SPEED);
			}
		}

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

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

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

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

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

	}
}