» 詳しい使用方法や、エラーで展開できない際の対応方法などはこちら
乱雑な密度分布における波のシミュレーション
このVCSSLプログラムは、乱雑な密度分布の媒質中を、波が伝播する様子をシミュレーションします。 密度境界で波が複雑に反射・屈折し、先端部が急速にエネルギーを失いながら進行する様子が見られます。
スポンサーリンク
使用方法
ダウンロードと展開(解凍)
まず、PC(スマホは未対応)で上の画面の「 ダウンロード 」ボタンを押してください。 するとZIP形式で圧縮されたファイルがダウンロードされます。
その後、ZIPファイルを右クリックして「すべて展開」や「ここに展開」などで展開(解凍)してください。 展開が成功すると、ZIPファイルと同じ名前のフォルダができ、その中にZIPファイルの中身が入っています。
» 展開がエラーで止まってしまう場合や、ファイル名が文字化けしてしまう場合は…
プログラムの起動
Windows をご使用の場合
上記でZIPファイルを展開したフォルダ内にある、以下のバッチファイルをダブルクリック実行してください:
もしプログラムを書き変えながら使いたい場合は、代わりに「 VCSSL_Editor__プログラム編集はこちら.bat 」を実行してください。
正常に起動できると、初回のみ、Java実行環境を入手するか等を尋ねられるので、適時答えて済ませると、プログラムが起動します。 2回目以降はすぐに起動します。
Linux 等をご使用の場合
ZIPファイルを展開したフォルダ内へコマンドライン端末で cd して、以下の通り入力して実行してください:
(プログラムの内容を書き変えながら使いたい場合は、代わりに VCSSL_Editor.jar を実行)
» javaコマンドが使用できない等のエラーが表示される場合は…
起動後の操作方法
起動すると、ウィンドウ上に画面が表示され、乱雑な高密度媒質を通過する様子がアニメーションで見られます。
画面ウィンドウ上では、下記のマウス操作を行えます。
- 左ドラッグ … 視点の回転
- 右ドラッグ … 視点の平行移動
- ホイールスクロール … 拡大/縮小
密度分布の変更
密度分布は、画像ファイル 「 density.png 」 として描いたものを読み込むようになっています。実際にこのシミュレーションで用いているのは以下の画像です:
この通り、乱雑な形状が描かれています。この画像ファイルの内容を描きかえる事により、自由な密度分布に設定できます。 画像ファイルのピクセルごとの明るさが、その点の媒質密度に対応します。 具体的には、暗いほど密度が重く、明るいほど軽くなります。
パラメータの変更
波に関する各種パラメータは、プログラムのコードの先頭領域にまとめて記載されています。 それらの値を書き換えると、色々と条件を変えてシミュレーションを楽しめます。
シミュレーションの題材解説
このプログラムでは、乱雑な密度分布の媒質を、波が通過する様子をシミュレーションしています。 ここでは、その際に生じる現象について、簡単に考察してみましょう。
乱雑な密度分布
上でも述べましたが、このシミュレーションで用いた媒質の密度分布は、以下のPNG形式画像ファイルとして読み込んでいます。

上記の図の、ピクセルごとの明度情報が、その点の密度分布に対応します。 暗い部分ほど密度が重く、明るい部分ほど軽くなります。
上図の通り、このシミュレーションの媒質は、 乱雑な(でたらめな)形状の細かい領域が組み合わさったような模様となっており、各領域の密度も乱雑です。 このような分布は、ミクロスケールのものでは多結晶、マクロスケールのものでは岩石などに見られます。
このシミュレーションを実行した様子が下図です。

上図において、波は左から右へ入射しています。最初はまとまった(パルス)波が入射したにもかかわらず、ばらけながら進んでいる様子が見て取れます。
もし密度分布が一様であれば、理想的な(線形な)波は、入射時の形のまま、ばらけずに進行します。
しかし乱雑な密度分布の中では、波は密度境界で複雑に反射・屈折を繰り返し、ばらけながら進行します。 結果、波の先端部は急速にエネルギーを失いながら進みます。
コード解説
このプログラムのコードは、プログラミング言語VCSSLで記述されています。
実は、コード内容は以下のページのプログラムと全く同じで、 密度分布を設定する画像ファイル( density.png )が異なるだけです。 モデルやシミュレーション方法、コードなどについての詳しい解説は、以下のページをご参照下さい。
ここでは、コード全体のみを掲載しておきます。
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[] = {pulse, sin};
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];
}
}
}
RandomDensityWave.vcssl
ライセンス
このVCSSL/Vnanoコード( 拡張子が「.vcssl」や「.vnano」のファイル )は実質的な著作権フリー(パブリックドメイン) である CC0 の状態で公開しています※。 記事中にC言語/C++/Java言語などでのサンプルコードが掲載されいてる場合は、それらについても同様です。 そのままでのご利用はもちろん、改造や流用などもご自由に行ってください。
※ ただし、このコードの配布フォルダ内には、ダウンロード後すぐに実行できるように、 VCSSLの実行環境も同梱されており、そのライセンス文書は「 License 」フォルダ内に同梱されています (要約すると、商用・非商用問わず自由に使用できますが、使用の結果に対して開発元は一切の責任を負いません、といった具合の内容です)。 配布フォルダ内の各構成物の一覧やライセンスについては「 ReadMe_使用方法_必ずお読みください.txt 」をご参照ください。
※ Vnano の実行環境については、別途スクリプトエンジンのソースコードも一般公開しており、 何らかのソフトウェア内に組み込んでご利用いただく事も可能です。詳細はこちらをご参照ください。
この記事中の商標などについて
- OracleとJavaは、Oracle Corporation 及びその子会社、関連会社の米国及びその他の国における登録商標です。文中の社名、商品名等は各社の商標または登録商標である場合があります。
- Windows は、米国 Microsoft Corporation の米国およびその他の国における登録商標です。この記事は独立著作物であり、Microsoft Corporation と関連のある、もしくはスポンサーを受けるものではありません。
- Linux は、Linus Torvalds 氏の米国およびその他の国における商標または登録商標です。
- その他、文中に使用されている商標は、その商標を保持する各社の各国における商標または登録商標です。
Vnano版 | ローレンツ方程式を数値的に解くプログラム |
|
![]() |
ローレンツ方程式を4次ルンゲ=クッタ法によって解き、グラフ描画用のデータを出力するプログラムです。 |
波の干渉(面上の円形波)のアニメーション表示 |
|
![]() |
面上の円形波が干渉する様子を、パラメータを操作しながらアニメーションで見られるプログラムです。 |
円形波のアニメーション表示 |
|
![]() |
振幅・波長・周期をスライダ―で操作しながら、円形波のグラフをアニメーションで見られるプログラムです。 |
波の干渉(線上の正弦波)のアニメーション表示 |
|
![]() |
線上(1次元の)の正弦波が干渉する様子を、パラメータを操作しながらアニメーションで見られるプログラムです。 |
正弦波のアニメーション表示 |
|
![]() |
振幅・波長・周期をスライダ―で操作しながら、正弦波のグラフをアニメーションで見られるプログラムです。 |
凹レンズを通過する波のシミュレーション |
|
![]() |
凹レンズ形状の高密度媒質を通過する、波のシミュレーションです。 |
凸レンズを通過する波のシミュレーション |
|
![]() |
凸レンズ形状の高密度媒質を通過する、波のシミュレーションです。 |
乱雑な密度分布における波のシミュレーション |
|
![]() |
密度分布が乱雑な媒質中における、波の伝播のシミュレーションです。 |
ローレンツアトラクタ(ファイル出力版) |
|
![]() |
4次精度ルンゲ=クッタ法により、ローレンツアトラクタを求めるプログラムです。 |
波の屈折のシミュレーション |
|
![]() |
密度の異なる領域を、波が屈折しながら通過するシミュレーションです。 |
力学アルゴリズムによる波のシミュレーション(面上の波) |
|
![]() |
媒質をバネと格子点で近似し、力学的なアルゴリズムで動かす事による、波のシミュレーションです。 |
手動で波を発生させるシミュレーション |
|
![]() |
スライダーをマウスで動かす事により、波を発生させるシミュレーションです。 |
力学アルゴリズムによる波のシミュレーション(線上の波) |
|
![]() |
媒質をバネと格子点で近似し、力学的なアルゴリズムで動かす事による、波のシミュレーションです。 |
二重振り子のシミュレーション |
|
![]() |
ラグランジュ方程式を用いた、二重振り子のシミュレーションです。 |
ローレンツアトラクタ(GUI版) |
|
![]() |
4次精度ルンゲ=クッタ法により、ローレンツアトラクタを求めるプログラムです。 |