【本サイトのご利用指針】
当サイトはSPAM対策等のためJavaScriptを使っています。
JavaScriptの実行を可能な状態にしてご利用下さい。
Please enable the execution of JavaScript!





































































【Excel VBA】Excel VBAで音響FDTDシミュレーション【ソース公開】
　 【リンクフリー】　私設研究所ネオテックラボ Neo-Tech-Lab.co.uk
【記載者】 【私設研究所Neo-Tech-Lab】　上田 智章
【掲載日】2010年12月19日
 
ここにチェックボックス型外部コンテンツ・メニューが入ります。



【ブラウザはGoogle Chrome】

●最新版ダウンロードはこちら

【お知らせ】

現在、特にお知らせする事項はありません。

【インターフェース2014年7月号】

インターフェース2014年5月号のAE関連の捕捉記事として書いた記事が、インターフェース2014年7月号に掲載されましたので、下記ダウンロードリンクを公開いたします。
簡易型対消滅壁も実装しています。

　　●順方向問題Execel VBAプログラムダウンロード

音響逆方向問題に関しては、『ポイントクラウドから音源座標を求める』を参考にしてください。
 今月号の特集とは全く無関係な記事です。2014年5月号のAEセンシング関連になります。

【図１.a】１次元音響FDTDシミュレーションの様子

【図１.b】２次元音響FDTDシミュレーションの様子

【図１.c】Excel VBAグラフィック表示用ライブラリの使い方

【はじめに】

　YouTubeの公開動画【Tidal Wave Simulation by Excel VBA -- 津波のシミュレーション(Excel VBA使用)】を見た方から本ツールを使いたいとのリクエストがあったので、複雑過ぎると手に負えないので、機能を限定して切り出したライブラリです。（最後に改訂したのが１年前の2009年12月、最初の版は４年程前2005/7/20頃のコーディングです。YouTube動画は厳密には物理演算の式が異なります。津波シミュレーションを行うには海底の３次元形状(深さ)や波高を考慮した２．５次元の演算を行う必要があります。）

　FDTD(Finite-difference time-domain method)とは本来は電磁場解析シミュレーションの一手法で、日本語では時間領域差分法、有限差分時間領域法などと訳されている手法のことです。FDTD法は偏微分方程式を離散化することによって得られる簡単な演算式を交互に使って複雑な波動現象をシミュレートすることができます。有限要素法(FEM)や境界要素法(BEM)等のプログラムと違ってFDTDは簡単な式の反復計算だけで複雑な現象のシミュレーションを行うことができる特徴があります。
　このページではMicrosoft Excelに標準実装されたVBA(Visual Basic for Application)を使って１次元/２次元の音響FDTDシミュレーションを行う方法をご紹介いたします。2005年7月1日に佐賀大学理工学部の村松和弘准教授（所属　佐賀大学 理工学部 電気電子工学科 電子システム工学講座）に音響FDTDの原理をご教授いただき、その後自力でExcel VBAでコーディングしたものです。Excel VBAにはグラフィック描画を行う機能が実装されていませんので、RGB24bitのBitMap形式ファイルを介してImageオブジェクトを使って表示するツールを制作しました。まずはこのExcelソースをダウンロードして実行してみてください。
【注意事項】このプログラムはExcelのマクロ(VBA)を使用していますのでマクロを実行可能にするためにExcelのセキュリティー設定を変更する必要があります。ご注意ください。

【２次元音響シミュレーションの実施例】パイプ中の音響伝搬現象
　丸い地球の表面が平面にしか見えないのと同様、パイプ中を伝搬して音源から離れれば離れるほど波面は平面に近づいてゆきます。パイプから出ると回折波も確認できます。

【追加情報】

【ブラウザで３次元CGを楽しもう！　WebGL】

　追加情報を示します。現在のところ、利用できるブラウザはGoogle Chromeだけのようですが、JavaScriptでかなり高機能のリアルタイム３次元グラフィックスを楽しめます。
　推奨する理由は、最新Google Chromeさえあれば、Windowsアクセサリのメモ帳だけで、JavaScript, HTMLの編集を行って３次元グラフィックスが楽しめるので、Excelを購入しなくてもいいからです。
　描画速度もDirectX9なのでかなり早いです。Excelより100倍速いどころではないでしょう。O3Dと違いプラグインも必要ありません。

【ダウンロード・リンク】

　●【Excel VBA版】音響FDTDシミュレーション・プログラム ≪Ver. 0.01≫ /2010/12/19/
　●【Excel VBA版】音響FDTDシミュレーション・プログラム ≪Ver. 0.02≫ /2010/12/22/
　　左右上下端を簡易吸音壁とし、デモにフェーズドアレイを追加したバージョンです。

【追記事項】
　もしかして、Excel VBAの３次元グラフィックス・ライブラリに興味をお持ちでしたら、こちらへどうぞ。ライブラリ公開中です。
　【関連する技術に関するサイト内リンク】 　『MikuMikuDanceのPMDモデルをExcel VBAで表示するぞ！』
　　●ポリゴンフィルとテクスチャーマッピング←【ライブラリ・ダウンロード】
　　●透視変換について
　　●光源計算について
　　●３次元CG【参考】

【事前準備①】

【事前準備②】

【事前準備③】

【事前準備④】

【事前準備】Excel VBAのマクロの実行を許可する

Microsoft Excelに標準実装されているVBA(Visual Basic for Application)を使うには以下の設定を行う必要があります。Excel 2007で説明します。
①Excel左上にある基本設定ボタンをクリックして、プルダウンメニューを表示し、最下段にある『Excelのオプション(I)』ボタンをクリックします。
②セキュリティーセンターの設定を選択し、『セキュリティーセンターの設定(T)』ボタンをクリックします。
③セキュリティーセンターのマクロの設定を選択し、オプションボタンを『すべてのマクロを有効にする』を選択し、『VBAプロジェクトオブジェクトモデルへのアクセスを信頼する(V)』にもチェックを入れて下さい。OKボタンをクリックしてオプション設定に戻ります。
④Excelのオプションで基本設定を選択し、[開発]タブをリボンに表示する(D)にチェックを入れ、OKボタンをクリックします。
⑤一旦、Excelを終了し、再度起動してください。リボンに開発タブが表示されたはずです。

【動作確認】

１次元音響FDTDの『音源駆動関数の表示』ボタンをクリックして、灰色の部分にグラフが表示されたらOKです。


【説明】

　ボタンをクリックするだけでプログラムがスタートします。
　『１次元音響FDTD』『２次元音響FDTD』『BitMapライブラリ』の３シートで構成されています。

【１次元音響FDTD】

　１次元音響FDTDシミュレーションで音波の固定端反射と自由端反射の違いを説明するデモプログラムになっています。

【１次元音響FDTDの原理】

　１次元音響FDTDはオイラーの式と連続の式を離散化して得られる反復式です。ここでｐが音圧、ｕ,ｖ,ｗは粒子速度、ρは媒質密度、ｋは体積弾性率(ρｃ^2)、ｃは音速です。
１次元では音圧ｐと粒子速度ｕに関する２つの方程式が得られます。音圧ｐと粒子速度ｕの位置は図のように半グリッドずれて交互に配置されているところがミソで、交互に方程式を演算すること(リープフロッグという。)で精度のよい演算を行うことができます。


【１次元音響FDTDのデモ】

　１次元音響FDTDシミュレーションで音波の固定端反射と自由端反射の違いを説明する以外に、在来技術であるアクティブ・ノイズ・キャンセラー(ANC、最後はControlという場合もある。)についてもデモを用意しています。
　ANC(アクティブ・ノイズ・キャンセラー)は制御点に逆位相音圧を人工的に加算する静音化方法ですが、P(x,t)に対して-P(x,t)を加えるということは、P(x,t)＝P(x,t)－P(x,t)＝０となり、制御点での音圧を常に０になるように維持する制御を意味します。
　３つめのボタン『逆位相音圧制御』は中央部分でANCを行うデモになっています。制御点で、逆位相の音を加算したので左から右に向かって伝搬してきたノイズは相殺されて制御点より左にはノイズが伝搬するのを阻止することができました。しかし、制御点が加算した逆位相音圧の音源位置となるので左に向かって逆位相音が放射される結果となります。このように在来技術のANCは自分は雑音を阻止することはできるのですが、逆に周囲の騒音は増加する結果となってしまう迷惑な静音化方法なのです。

　プログラムの実行後に、プログラムを格納していたフォルダーの中を見て下さい。
『..........な、なんじゃ、こりゃあ！？』
山ほどbmp形式の画像ファイルが作られていると思います。実はこのプログラムはWindowsに標準添付のWindowsムービーメーカーで動画を作るために演算結果をbmpファイル形式で残す仕様になっています。数字はシミュレーションの時刻ステップの番号を示しています。
なので、十分にディスクスペースのあるパソコンをお使い下さい。(^_^;)

【デモに関連する知識】

　自由端反射と固定端反射、ANC(Active Noize Canceller)について少しだけ図を作りました。ANC自体はもう40年以上前の技術(当時はアナログ技術)ではないかと思います。DSPで演算を実行したりしていますが、静音化できるのは逆位相音を出すスピーカー周辺に限定され、拘束性があるのが欠点です。



【VBAソースの閲覧方法】

　リボンの右端にある[開発]タブをクリックし、左上端のVisual Basicをクリックします。
　手元に動作可能状態のExcel 2003はもはや存在しないので、すっかり忘れてしまいましたがツールメニューのマクロからVisual Basicに入れたはずです。(=_=)　≪忘却≫
　Win32APIという名前のモジュールもありますが、これは自作ハードウェアによる計測でUSB経由でExcell Sheetにデータを読み込んだりするときに使うもので、削除しても問題ないはずです。




２次元音響FDTDシミュレーション

　用意した２次元FDTDのプログラムは、24ビット ビットマップ形式の画像ファイルで定義された任意の２次元空間の音響伝搬現象をシミュレートすることができます。デフォルトで４パターン用意しました。４つのうちどれかをクリックすると、シミュレーション対象を変更することができます。

【２次元音響FDTDの原理】

　２次元FDTDは先ほどと同様オイラーの式と連続の式を離散化して得られます。音圧と粒子速度の格子の接続関係は下のようになっています。


【別のモデルの２次元音響FDTDシミュレーションを行うには】

　２次元FDTDの演算モデル空間は、Windowsアクセサリのペイントを使って作成し、24ビット ビットマップ形式(.bmp)で保存したものが使えるようにしています。デモには４パターン用意しています。
　プログラムを実行させるとわかりますが、『白色』の空間を媒質として取り扱います。他の色が付いている部分は『密度の濃い物体』が存在し、Version 0.01では境界面で固定端反射が起こるようになっています。このソフトは機能限定版なので端も固定端反射を起こします。
　(Version 0.02では私独自の簡易吸音壁[＝一次の対消滅壁]に変更しました。)
　デフォルトのソースでは２次元FDTDの演算モデル空間の大きさはNx=200, Ny=200としています。Sheet2(２次元音響FDTD)のソースのPrivate Constを変更すれば変更できます。
　音源の座標もSx,Syで与えていますので、適宜駆動音源の位置を変更することができます。
　自分で作成した２次元FDTDの演算モデル空間を読込むには以下のようにします。
　　■リボンの[開発]タブを選択し、『デザインモード』に変更します。
　　■次に『プロパティー』を表示します。
　　■モデルに使用する画像データを格納するImageオブジェクトをクリックして選択します。
　　■プロパティーのPictureの部分(ビットマップ)をクリックします。
　　■右端に『･･･』と書かれたボタンが表示されるので、これをクリックするとファイル『ピクチャの読込』ダイアログが表示されるので目的のファイルを選んで『開く(O)』ボタンをクリックすればOKです。

シミュレーション結果の動画の作り方

　Windows標準添付のムービーメーカーを使えば簡単です。
　まず、ムービーメーカーを起動したらメニューから『ツール(T)』－『オプション(O)』を実行します。オプション・ダイアログの詳細タブを選び、規定の再生時間を最小の0.125秒に設定して『OK』ボタンをクリックします。
　ムービーメーカーに編集対象の画像ファイルを読込むには『メディアの読込』ボタンをクリックして読込むこともできますが、対象ファイルを全部選択してムービーメーカー内にドラッグしてくれば読込めます。
　読込んだ画像ファイルを再び全部選択してビデオのタイムライン内にドラッグすれば番号順に読込まれます。
　特殊効果やタイトル、音楽の編集の後で、『ムービーの発行』をクリックしてコンピュータ内に.WMV形式のムービーを作成することができます。


作成した動画の例

【回折波】デモから

【パイプ中の音波伝搬現象】デモから

【フェーズドアレイ】

　わずかにソースをいじって制作した動画です。
　７素子のトランスデューサを駆動するタイミングをわずかにずらせてやると音波の進行方向をフレキシブルに変更してやることができます。これをPhased Arrayと言います。メカ的にトランスデューサの方向を変えることなく、音波の方向を変えることができるので医療用超音波断層撮影装置などに使われている技術です。

【フェーズドアレイ２】

　さらにソースをいじって上下左右端を私オリジナルの『簡易吸音壁』にしてみた例です。フェーズドアレイの素子数を６１に増やしました。

【センチネルリンパ節を磁性微粒子の励磁音響効果で検出する】

　ちょっとソースをいじれば下のような感じでシミュレーション結果を動画化することができます。

【Active Wave Barrier】≪NTLの真面目な研究内容のひとつ≫

騒音と逆位相の音を出して騒音をキャンセルするActive Noise Controlが知られています。
しかし、原理的に加算した逆位相音が騒音の入射方向に拡散するため、透過波は抑制できますが、
防音壁として利用しようとしても大きな反響音が出てしまう欠点がありました。
この欠点を克服するため、開発した技術が対消滅壁です。
自由端反射と固定端反射を同一箇所で実装することにより、透過波の抑制は勿論、反射波も抑制することに成功しました。
この技術を採用することで、音響、電磁波、ガンマ線などの波動現象におけるバリア、あるいはステルス人工素材が実現します。
アクティブ防音壁、人工静寂空間、防音室といった音響系技術だけでなく、電磁波、光波領域でも利用が期待できます。

対消滅壁には幾つかの実現方法がありますが、
ソース公開している【Excel VBA版】音響FDTDシミュレーション・プログラム ≪Ver. 0.02≫にも簡易版を実装しています。
通常の音響FDTD(【Excel VBA版】音響FDTDシミュレーション・プログラム ≪Ver. 0.01≫)では
上下左右端で大きな反射現象が発生してしまう欠点がありますが、Ver.0.02では反射が抑制されていることがわかると思います。
実はいわゆる多層減衰壁ではなく、私オリジナルの『対消滅壁』という技術なのです。

対消滅といってもエネルギーが消えるわけではありません。
エネルギー保存則は成立します。
自由端反射、固定端反射のそれぞれによって波は拡散してゆきますが、
近接して逆位相(音圧と粒子速度の両方)の反射波が存在するために
あたかも消失したかのように見えるだけなのです。
これはビッグバンによって生じた物質と反物質が対消滅によってあたかも存在しなくなったかのように見えながら
その実、宇宙にエネルギーとして残存しているのと同じことなのです。
音響FDTD Ver.0.02でも上下左右端からの反射は激減したかのように見えますが、エネルギーは残留し、伝播を続けています。
この対消滅によって拡散させた波動を逆拡散するプロセスを確立することができれば
人類はエネルギー問題から開放されることでしょう。

　●【Excel VBA版】音響FDTDシミュレーション・プログラム ≪Ver. 0.01≫ /2010/12/19/
　●【Excel VBA版】音響FDTDシミュレーション・プログラム ≪Ver. 0.02≫ /2010/12/22/


【３次元音響シールドへの適用】

【対消滅壁の原理説明部分の動画】

【標準モジュール：NTL_Lib_2DCG_VBA001のソース】

点と直線、ボックスしか描画できない極めて限定された機能のみの２次元グラフィックスライブラリです。
基本色や指定色(RGB24bit)による描画の他、シミュレーション結果の数値を色に変換するのに便利なカラーインデックス方式での描画が可能です。
'%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
'% 私設研究所 Neo-Tech-Lab.com Excel VBA ２次元ＣＧライブラリ Ver. 0.01
'% 【作成者】東京工業大学ソリューション研究機構　上田智章
'% 【作成日】2010年12月17日
'% 【目　的】2次元FDTDシミュレーション用ライブラリ
'% 【仕　様】上記目的での利用に限定した点と直線描画のみ行う簡易版(basic)
'% 【ご利用条件】
'% 　●本ライブラリによる演算結果は非商用、商用を問わず自由にご利用になれます。
'% 　●このライブラリを使うことによって発生した損害の責任はユーザー自身で負う事。
'% 　●多忙の為、バグ対応、質疑応答を行う時間がありませんので対応いたしません。
'% 　●本ライブラリの２次配付はお断りいたします。
'%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
'【説明】
'　本ライブラリは、Excel VBAに標準実装されているImageオブジェクトを使ってイメージ
'を表示したり、ファイルに保存することを前提に作られています。描画/読込/保存処理を
'簡単に行うために汎用性を捨てて、Windows標準のビットマップ形式のうちRGB24bit形式
'のみに特化してソフトウェアを製作いたしました。このため表示結果を保存するディスク
'スペースが他の画像形式に比べ、各段に必要となります。
'　動画制作の際には他のツール等を使ってbmpファイルをjpg, gifあるいはpng形式などに
'一括変換されることをお薦めいたします。予めご了承ください。
'--------------------------------------------------------------------------------
'【リファレンス】
'■Public Sub GetBitMapFile(cg As NTL_BitMap, FileName As String)
'【機能】指定ファイル(.bmp)のイメージをBitMapオブジェクトに読込む
'【事例】GetBitMapFile MyPic,"Doll.bmp" '【"Doll.bmp"をMyPicに読込む】
'
'■Public Sub BuildBitMap(Nx As Long, Ny As Long, cg As NTL_BitMap)
'【機能】指定サイズ(横画素数Nx, 縦画素数Ny)のBitMapオブジェクトを作成する
'【事例】BuildBitMap 320, 240, MyPic '【横320画素、縦240画素のBitMapオブジェクトを作る】
'
'■Public Sub FillBitMapImage(cg As NTL_BitMap, P As RGB24bitPixel)
'【機能】BitMapオブジェクトのFrameBufferを指定色Pで初期化する
'【事例】FillBitMapImage MyPic, MyPic.DarkGray '【BitMapオブジェクトのFrameBufferをDarkGrayで塗り潰す】
'
'■Public Sub CreateBitMapFile(cg As NTL_BitMap, FileName As String)
'【機能】BitMapオブジェクトのFrameBufferの画像をファイルにする
'【事例】CreateBitMapFile MyPic, "C:\Photo\Sample.bmp" '【BitMapオブジェクトの画像を"Sample.bmp"にファイル化】
'
'■Public Sub InitializeColorLookupTable(cg As NTL_BitMap, n As Long)
'【機能】BitMapオブジェクトにｎ色のカラールックアップテーブルを用意する
'【事例】InitializeColorLookupTable MyPic, 1024 '【BitMapオブジェクトに1024色のカラールックアップテーブルを用意する】
'
'■Public Sub CreateColor(cg As NTL_BitMap, i1 As Long, P1 As RGB24bitPixel, i2 As Long, P2 As RGB24bitPixel)
'【機能】BitMapオブジェクトのカラールックアップテーブル(CLUT)のi1番目からi2番目までにP1からP2までの色を直線補間で用意する
'【事例】CreateColor MyPic, 105, MyPic.Red, 210, MyPic.Yellow '【BitMapオブジェクトのCLUTの105番から210番までに赤から黄色まで用意する】
'
'■Public Sub CreateDipoleScale(cg As NTL_BitMap, n As Long)
'【機能】正負表示用のデフォルトのｎ色カラールックアップテーブルを準備する
'【事例】CreateDipoleScale MyPic, 2000 '【正負表示用のデフォルトのカラールックアップテーブルを2000色で準備する】
'
'■Public Sub CreateGrayScale(cg As NTL_BitMap, n As Long)
'【機能】ｎ階調グレイスケール用のカラールックアップテーブルを準備する
'【事例】CreateGrayScale MyPic, 256 '256階調グレイスケール用のカラールックアップテーブルを準備する
'
'■Public Sub CreateMonopoleScale(cg As NTL_BitMap, n As Long)
'【機能】絶対値表示用のデフォルトのｎ色カラールックアップテーブルを準備する
'【事例】CreateMonopoleScale MyPic, 1024 '【絶対値表示用の1024色カラールックアップテーブルを準備する】
'
'■Public Sub DrawPixel(cg As NTL_BitMap, x As Long, y As Long, P As RGB24bitPixel)
'【機能】BitMapオブジェクトの指定座標(x, y)に指定色Pの画素を描画する
'【事例】DrawPixel MyPic, 100, 20, MyPic.Orange '【BitMapオブジェクトの(100,20)に橙色の画素を描画する】
'
'■Public Function GetPixel(cg As NTL_BitMap, x As Long, y As Long) As RGB24bitPixel
'【機能】BitMapオブジェクトの指定座標(x, y)の画素の色情報を読込む
'【事例】P = GetPixel(MyPic, x, y) '【注】BitMapオブジェクトの範囲外の場合は黒色が返る
'
'■Public Sub DrawLine(cg As NTL_BitMap, X1 As Long, Y1 As Long, X2 As Long, Y2 As Long, P As RGB24bitPixel)
'【機能】BitMapオブジェクトの指定座標(X1,Y1)と(X2,Y2)を結ぶ指定色Pの直線を描画する
'【事例】DrawLine MyPic, 25, 70, 87, 10, MyPic.Brown '【(25,70)-(87,10)を結ぶ茶色の直線を描画する】
'
'■Public Sub DrawBox(cg As NTL_BitMap, X1 As Long, Y1 As Long, X2 As Long, Y2 As Long, P1 As RGB24bitPixel, P2 As RGB24bitPixel)
'【機能】(X1,Y1)と(X2,Y2)を対角線とするボックスを枠線色P1で描画し、内部をP2で塗り潰す
'【事例】DrawBox MyPic, 10, 10, 100, 20, MyPic.Blue, MyPic.Green
'
'■Public Sub ImageCopy(dstcg As NTL_BitMap, dstx As Long, dsty As Long, srccg As NTL_BitMap, srcx0 As Long, srcy0 As Long, srcxw As Long, srcyw As Long, T As RGB24bitPixel)
'【機能】ソース画像の(srcx0,srcy0)を始点とする幅(srcxw,srcyw)の領域をディスティネーション画像の(dstx, dsty)を始点とする領域にコピーする。
'        但し、指定色Tを持つ画素はコピーしないので透明として取り扱われる
'【事例】ImageCopy MyPic1, 10, 10, MyPic2, 100, 20, 16, 16, MyPic.Green

'【BitMap画素構造体】並び順に注意！BitMapのファイル形式(.bmp)に依存しています。
Public Type RGB24bitPixel ' 画素ﾃﾞｰﾀの構造(RGB24bitﾀｲﾌﾟ)
   Blue  As Byte ' 青(0～255)
   Green As Byte ' 緑(0～255)
   Red   As Byte ' 赤(0～255)
End Type

'【BitMapヘッダー情報構造体】BitMap形式のファイル(.bmp)のヘッダー情報部分
Public Type RGB24bitBitMapHeader ' RGB24bitﾀｲﾌﾟのBitMapﾌｧｲﾙのﾍｯﾀﾞｰ
  ID1 As Byte                    ' ﾌｧｲﾙ識別子 "B"
  ID2 As Byte                    '      :     "M"
  FileLength As Long             ' ﾌｧｲﾙの長さ = ﾍｯﾀﾞｰｻｲｽﾞ (54ﾊﾞｲﾄ) + ﾃﾞｰﾀｻｲｽﾞ( (x方向画素数×3に最も近い4の倍数)×(y方向画素数) )
  Null1 As Long                  ' 0
  HeaderSize As Long             ' ﾍｯﾀﾞｰ領域のｻｲｽﾞ        (54ﾊﾞｲﾄ)
  Offset As Long                 ' 画素ﾃﾞｰﾀまでのｵﾌｾｯﾄｻｲｽﾞ(40ﾊﾞｲﾄ)
  Nx As Long                     ' x方向画素数
  Ny As Long                     ' y方向画素数
  NumberOfPlanes As Integer      ' ﾌﾟﾚｰﾝの数              (1ﾌﾟﾚｰﾝ)
  BitsOfPixel As Integer         ' 1画素を構成するﾋﾞｯﾄ数  (24ﾋﾞｯﾄ)
  Null2 As Long                  ' 0
  SizeOfData As Long             ' 画素領域のﾊﾞｲﾄ･ｻｲｽﾞ    (x方向画素数×3に最も近い4の倍数)×(y方向画素数)
  Null3 As Long                  ' 0
  Null4 As Long                  ' 0
  Null5 As Long                  ' 0
  Null6 As Long                  ' 0
End Type

'【NTL ２次元描画用オブジェクトの構造体】
Public Type NTL_BitMap
'【BitMapファイル部分】
   Header As RGB24bitBitMapHeader ' RGB24bitﾀｲﾌﾟのBitMapﾍｯﾀﾞｰ
   PixelBuffer() As Byte          ' RGB画素データ格納領域
'【BitMap形式ファイルの横方向バイト数】x,y座標をアクセスしやすくするために２次元配列の領域
   BytesOfScaneLine As Long ' (x方向画素数×3に最も近い4の倍数) = (y方向に隣接する画素までのﾊﾞｲﾄ距離)
'【描画領域】x,y座標をアクセスしやすくするために２次元配列の領域
   FrameBuffer() As RGB24bitPixel '引数は自然数
'【描画用基本色】簡単に色指定するためのもの
'【基本色名称】                         R   G   B
   White     As RGB24bitPixel ' 白    255 255 255
   Black     As RGB24bitPixel ' 黒      0   0   0
   Red       As RGB24bitPixel ' 赤    255   0   0
   Orange    As RGB24bitPixel ' ｵﾚﾝｼﾞ 255 128   0
   Yellow    As RGB24bitPixel ' 黄    255 255   0
   Green     As RGB24bitPixel ' 緑      0 255   0
   Cyan      As RGB24bitPixel ' ｼｱﾝ     0 255 255
   Blue      As RGB24bitPixel ' 青      0   0 255
   Violet    As RGB24bitPixel ' 紫    255   0 128
   Magenta   As RGB24bitPixel ' ﾏｾﾞﾝﾀ 255   0 255
   Brown     As RGB24bitPixel ' 茶    128   0   0
   DarkGray  As RGB24bitPixel ' 濃灰　 50  50  50
   Gray      As RGB24bitPixel ' 灰　　128 128 128
   LightGray As RGB24bitPixel ' 淡灰  200 200 200
   DarkGreen As RGB24bitPixel ' 濃緑　  0 128   0
'【逆引きColorLookUpTable】数値を色に対応付けるためのデータテーブル
   nCLUT As Long           ' ColorLookUpTableの登録色数
   CLUT() As RGB24bitPixel ' ColorLookUpTableのﾃﾞｰﾀ領域
End Type

'%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
'%%% x方向画素数×3に最も近い4の倍数を求める。
'%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Private Function GetBytesOfScanLine(Nx As Long) As Long
   Dim i As Long, j As Long
   
   i = Nx * 3&               ' ScanLine方向の画素情報のバイト数
   j = i Mod 4&              ' ４バイト境界での余剰バイト数
   If j > 0 Then             ' ４バイト境界での余剰バイトがあるとき
      i = (i \ 4& + 1&) * 4& ' １ワード余計に必要
   End If
   GetBytesOfScanLine = i
End Function

'%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
'%%% 基本色を作る
'%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Private Sub SetBasicColor(cg As NTL_BitMap)
' Obj.色名.xxx                        R   G   B
   cg.White.Red = 255       ' 白    255 255 255
   cg.White.Green = 255     '
   cg.White.Blue = 255      '
   cg.Black.Red = 0         ' 黒      0   0   0
   cg.Black.Green = 0       '
   cg.Black.Blue = 0        '
   cg.Red.Red = 255         ' 赤    255   0   0
   cg.Red.Green = 0         '
   cg.Red.Blue = 0          '
   cg.Yellow.Red = 255      ' 黄    255 255   0
   cg.Yellow.Green = 255    '
   cg.Yellow.Blue = 0       '
   cg.Green.Red = 0         ' 緑      0 255   0
   cg.Green.Green = 255     '
   cg.Green.Blue = 0        '
   cg.Cyan.Red = 0          ' ｼｱﾝ     0 255 255
   cg.Cyan.Green = 255      '
   cg.Cyan.Blue = 255       '
   cg.Blue.Red = 0          ' 青      0   0 255
   cg.Blue.Green = 0        '
   cg.Blue.Blue = 255       '
   cg.Magenta.Red = 255     ' ﾏｾﾞﾝﾀ 255   0 255
   cg.Magenta.Green = 0     '
   cg.Magenta.Blue = 255    '
   cg.Orange.Red = 255      ' ｵﾚﾝｼﾞ 255 128   0
   cg.Orange.Green = 128    '
   cg.Orange.Blue = 0       '
   cg.Violet.Red = 128      ' 紫    128   0 255
   cg.Violet.Green = 0      '
   cg.Violet.Blue = 128     '
   cg.Brown.Red = 128       ' 茶    128   0   0
   cg.Brown.Green = 0       '
   cg.Brown.Blue = 0        '
   cg.Gray.Red = 128        ' 灰    128 128 128
   cg.Gray.Green = 128      '
   cg.Gray.Blue = 128       '
   cg.DarkGray.Red = 50     ' 濃灰　 50  50  50
   cg.DarkGray.Green = 50   '
   cg.DarkGray.Blue = 50    '
   cg.LightGray.Red = 200   ' 淡灰  200 200 200
   cg.LightGray.Green = 200 '
   cg.LightGray.Blue = 200  '
   cg.DarkGreen.Red = 0     ' 濃緑　  0 128   0
   cg.DarkGreen.Green = 128 '
   cg.DarkGreen.Blue = 0    '
End Sub

'%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
'%%% フレームバッファの画像データをピクセルバッファにコピーする
'%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
'【機能】NTL_BitMap形式のオブジェクト(cg)に指定ファイル(FileName)を読込む
Private Sub CopyFrameBufferDataToPixelBuffer(cg As NTL_BitMap)
   Dim x As Long, y As Long, i As Long, P As RGB24bitPixel
   
   For y = 0 To cg.Header.Ny - 1
      For x = 0 To cg.Header.Nx - 1
         i = x * 3& + cg.BytesOfScaneLine * y '【画素のバイト位置を求める】
         P = cg.FrameBuffer(x, y)             '【フレームバッファから画素データを読込む】
         cg.PixelBuffer(i) = P.Blue           '【画素バッファに青データをコピー】
         cg.PixelBuffer(i + 1) = P.Green      '【画素バッファに緑データをコピー】
         cg.PixelBuffer(i + 2) = P.Red        '【画素バッファに赤データをコピー】
      Next x
   Next y
End Sub

'%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
'%%% ピクセルバッファの画像データをフレームバッファにコピーする
'%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
'【機能】NTL_BitMap形式のオブジェクト(cg)に指定ファイル(FileName)を読込む
Private Sub CopyPixelBufferDataToFrameBuffer(cg As NTL_BitMap)
   Dim x As Long, y As Long, i As Long, P As RGB24bitPixel
   
   For y = 0 To cg.Header.Ny - 1
      For x = 0 To cg.Header.Nx - 1
         i = x * 3& + cg.BytesOfScaneLine * y '【画素のバイト位置を求める】
         P.Blue = cg.PixelBuffer(i)           '【画素バッファから青データを読込む】
         P.Green = cg.PixelBuffer(i + 1)      '【画素バッファから緑データを読込む】
         P.Red = cg.PixelBuffer(i + 2)        '【画素バッファから赤データを読込む】
         cg.FrameBuffer(x, y) = P             '【フレームバッファに画素データをコピー】
      Next x
   Next y
End Sub

'###########################################################
'### 指定画像ファイル(RGB24bit BMP形式)を読み込む
'###########################################################
'【機能】NTL_BitMap形式のオブジェクト(cg)に指定ファイル(FileName)を読込む
'【注意】ファイル形式は必ずRGB24bitの.bmp形式でなければならない
'　　　　動作速度を優先してエラーチェックしていないので注意の事
Public Sub GetBitMapFile(cg As NTL_BitMap, FileName As String)
   Open FileName For Binary Access Read As #1        '【ファイルを開く】
      Get #1, , cg.Header                            '【BitMap形式ファイルのヘッダー情報を読込む】
      ReDim cg.PixelBuffer(cg.Header.SizeOfData - 1) '【画素データの格納領域を確保する】
      Get #1, 55, cg.PixelBuffer                     '【画素データを読込む】
   Close #1
   cg.BytesOfScaneLine = GetBytesOfScanLine(cg.Header.Nx)   '【x方向画素数×3に最も近い4の倍数を求める】
   ReDim cg.FrameBuffer(cg.Header.Nx - 1, cg.Header.Ny - 1) '【フレームメモリの領域を確保する】
   CopyPixelBufferDataToFrameBuffer cg '【PixelBuffer⇒FrameBuffer】
   SetBasicColor cg                    '【基本色を設定する】
End Sub

'###########################################################
'### 指定サイズの画像ファイル(RGB24bit BMP形式)を作る
'###########################################################
'【機能】x方向画素数Nx、y方向画素数NyのNTL_BitMap形式のオブジェクト(cg)を作成する
Public Sub BuildBitMap(Nx As Long, Ny As Long, cg As NTL_BitMap)
   Dim i As Long, j As Long
   
   cg.BytesOfScaneLine = GetBytesOfScanLine(Nx)   '【ScanLine方向のバイト数を求める】
   SetBasicColor cg                                '【基本色を設定する】
   With cg.Header                                  '
      .ID1 = Asc("B")                              '【識別子の第１バイトは"B"】
      .ID2 = Asc("M")                              '【識別子の第２バイトは"M"】
      .FileLength = 54& + cg.BytesOfScaneLine * Ny '【.bmpファイル全体のバイト数】
      .Null1 = 0&                                  '
      .HeaderSize = 54                             '【.bmpファイルのヘッダー情報のバイト数】
      .Offset = 40                                 '【画素データまでのバイト数】
      .Nx = Nx                                     '【x方向の画素数】
      .Ny = Ny                                     '【y方向の画素数】
      .NumberOfPlanes = 1                          '【Plane数】
      .BitsOfPixel = 24                            '【１画素あたりのビット数】
      .Null2 = 0&                                  '
      .SizeOfData = cg.BytesOfScaneLine * Ny       '【画素バッファのバイト数】
      .Null3 = 0&                                  '
      .Null4 = 0&                                  '
      .Null5 = 0&                                  '
      .Null6 = 0&                                  '
   End With                                        '
   ReDim cg.PixelBuffer(cg.Header.SizeOfData - 1)  '【PixelBuffer領域確保】
   ReDim cg.FrameBuffer(Nx - 1, Ny - 1)            '【FrameBuffer領域確保】
   FillBitMapImage cg, cg.White                    '【FrameBufferを白色でクリアする】
End Sub

'###########################################################
'### フレームバッファのカレントイメージを指定色(Pixel)でクリアする
'###########################################################
Public Sub FillBitMapImage(cg As NTL_BitMap, P As RGB24bitPixel)
   Dim x As Long, y As Long
   
   For y = 0 To cg.Header.Ny - 1
     For x = 0 To cg.Header.Nx - 1
        cg.FrameBuffer(x, y) = P   '【指定色を代入】
     Next x
   Next y
End Sub

'###########################################################
'### フレームバッファのカレントイメージをファイル化する
'###########################################################
Public Sub CreateBitMapFile(cg As NTL_BitMap, FileName As String)
   CopyFrameBufferDataToPixelBuffer cg '【FrameBuffer⇒PixelBuffer】
   Open FileName For Binary Access Write As #1
      Put #1, 1, cg.Header
      Put #1, 55, cg.PixelBuffer
   Close #1
End Sub

'###########################################################
'### ColorLookUpTableを初期化してｎ色分の領域を確保する
'###########################################################
Public Sub InitializeColorLookupTable(cg As NTL_BitMap, n As Long)
   cg.nCLUT = n
   ReDim cg.CLUT(n - 1)
End Sub

'###########################################################
'### ColorLookUpTableの指定パレット番号i1からi2までを指定色P1からP2まで直線補間した色を作って設定する
'###########################################################
Public Sub CreateColor(cg As NTL_BitMap, i1 As Long, P1 As RGB24bitPixel, i2 As Long, P2 As RGB24bitPixel)
   Dim dR As Double, dG As Double, dB As Double, i As Long, d As Double

   d = 1# / CDbl(i2 - i1)
   dR = (CDbl(P2.Red) - CDbl(P1.Red)) * d
   dG = (CDbl(P2.Green) - CDbl(P1.Green)) * d
   dB = (CDbl(P2.Blue) - CDbl(P1.Blue)) * d
   For i = 0 To (i2 - i1) Step Sgn(i2 - i1)
      cg.CLUT(i1 + i).Red = P1.Red + dR * CDbl(i)
      cg.CLUT(i1 + i).Green = P1.Green + dG * CDbl(i)
      cg.CLUT(i1 + i).Blue = P1.Blue + dB * CDbl(i)
   Next i
End Sub

'###########################################################
'### ｶﾗｰﾊﾟﾚｯﾄを作る(ｎ階調のｸﾞﾚｲ･ｽｹｰﾙ)　白⇒黒
'###########################################################
Public Sub CreateGrayScale(cg As NTL_BitMap, n As Long)
   InitializeColorLookupTable cg, n
   CreateColor cg, 0, cg.Black, n - 1, cg.White
End Sub

'###########################################################
'### ｶﾗｰﾊﾟﾚｯﾄを作る(正負)　赤⇒橙⇒黄⇒緑|濃緑⇒ｼｱﾝ⇒青⇒紫
'###########################################################
Public Sub CreateDipoleScale(cg As NTL_BitMap, n As Long)
   Dim C0 As Long, C1 As Long, C2 As Long, C3 As Long, C4 As Long, C5 As Long
   
   C0 = n - 1
   C1 = n * 0.84
   C2 = n * 0.67
   C3 = n * 0.5
   C4 = n * 0.33
   C5 = n * 0.17
   InitializeColorLookupTable cg, n
   CreateColor cg, C0, cg.Red, C1, cg.Orange
   CreateColor cg, C1, cg.Orange, C2, cg.Yellow
   CreateColor cg, C2, cg.Yellow, C3, cg.Green
   CreateColor cg, C3, cg.DarkGreen, C4, cg.Cyan
   CreateColor cg, C4, cg.Cyan, C5, cg.Blue
   CreateColor cg, C5, cg.Blue, 0, cg.Violet
   cg.CLUT(C3) = cg.White
End Sub

'###########################################################
'### ｶﾗｰﾊﾟﾚｯﾄを作る(絶対値)
'###########################################################
Public Sub CreateMonopoleScale(cg As NTL_BitMap, n As Long)
   Dim C0 As Long, C1 As Long, C2 As Long, C3 As Long, C4 As Long, C5 As Long, C6 As Long
   
   C0 = n - 1
   C1 = n * 10 / 12
   C2 = n * 8 / 12
   C3 = n * 6 / 12
   C4 = n * 5 / 12
   C5 = n * 4 / 12
   C6 = n * 2 / 12
   InitializeColorLookupTable cg, n
   ' 表示色数　赤⇒橙⇒黄⇒緑⇒淡青⇒青⇒濃灰⇒灰⇒淡灰⇒白
   CreateColor cg, C0, cg.Red, C1, cg.Orange
   CreateColor cg, C1, cg.Orange, C2, cg.Yellow
   CreateColor cg, C2, cg.Yellow, C3, cg.Green
   CreateColor cg, C3, cg.Green, C4, cg.DarkGreen
   CreateColor cg, C4, cg.DarkGreen, C5, cg.Blue
   CreateColor cg, C5, cg.Blue, C6, cg.Violet
   CreateColor cg, C6, cg.DarkGray, 0, cg.White
End Sub

'###########################################################
'### カレントイメージに指定色で点を描画する
'###########################################################
Public Sub DrawPixel(cg As NTL_BitMap, x As Long, y As Long, P As RGB24bitPixel)
   If (x >= 0) And (x < cg.Header.Nx) And (y >= 0) And (y < cg.Header.Ny) Then
      cg.FrameBuffer(x, y) = P
   End If
End Sub

'###########################################################
'### カレントイメージの指定位置から色情報を取得する
'###########################################################
Public Function GetPixel(cg As NTL_BitMap, x As Long, y As Long) As RGB24bitPixel
   If (x >= 0) And (x < cg.Header.Nx) And (y >= 0) And (y < cg.Header.Ny) Then
      GetPixel = cg.FrameBuffer(x, y)
   Else
      GetPixel = cg.Black
   End If
End Function

'###########################################################
'### カレントイメージに指定色で直線を描画する
'###########################################################
Public Sub DrawLine(cg As NTL_BitMap, X1 As Long, Y1 As Long, X2 As Long, Y2 As Long, P As RGB24bitPixel)
   Dim x As Long, y As Long, k As Long, dX As Long, dY As Long, Sx As Long, Sy As Long

   dX = X2 - X1              ' 3+
   dY = Y2 - Y1              '  |
   Sx = Sgn(dX)              ' 2+
   If Sx = 0 Then Sx = 1&    '  |
   Sy = Sgn(dY)              ' 1+
   If Sy = 0 Then Sy = 1&    '  |
   dX = Abs(dX)              ' 0+---+---+---+---+
   dY = Abs(dY)              '  0   1   2   3   4
   If dX > dY Then           ' Major=X, Minor=Y
      y = Y1                 ' Minorの初期値
      k = 0
      For x = X1 To X2 Step Sx
         k = k + dY          ' Minor増加分だけ加算
         If k > dX Then      ' Major分を超えたらMajor分だけ減算
            k = k - dX       '
            y = y + Sy       ' Minor座標方向の移動バイト数を計算
         End If
         If (x >= 0) And (x < cg.Header.Nx) And (y >= 0) And (y < cg.Header.Ny) Then
            cg.FrameBuffer(x, y) = P
         End If
      Next x
   Else                      ' Major=Y, Mior=X
      x = X1                 ' Minorの初期値
      k = 0
      For y = Y1 To Y2 Step Sy
         k = k + dX          ' Minor増加分だけ加算
         If k > dY Then      ' Major分を超えたらMajor分だけ減算
            k = k - dY       '
            x = x + Sx       ' Minor座標方向の移動バイト数を計算
         End If
         If (x >= 0) And (x < cg.Header.Nx) And (y >= 0) And (y < cg.Header.Ny) Then
            cg.FrameBuffer(x, y) = P
         End If
      Next y
   End If
End Sub

'###########################################################
'### カレントイメージにボックスを描画する
'###########################################################
'【機能】枠線を指定色P1で描画し、枠内を指定色P2で塗り潰す
Public Sub DrawBox(cg As NTL_BitMap, X1 As Long, Y1 As Long, X2 As Long, Y2 As Long, P1 As RGB24bitPixel, P2 As RGB24bitPixel)
   Dim x As Long, y As Long, xx1 As Long, yy1 As Long, xx2 As Long, yy2 As Long
   
   xx1 = X1
   xx2 = X2
   yy1 = Y1
   yy2 = Y2
   If xx1 > xx2 Then
      i = xx1
      xx1 = xx2
      xx2 = i
   End If
   If yy1 > yy2 Then
      i = yy1
      yy1 = yy2
      yy2 = i
   End If
   For x = xx1 To xx2 '【上下枠線を指定色P1で描画】
      y = yy1
      If (x >= 0) And (x < cg.Header.Nx) And (y >= 0) And (y < cg.Header.Ny) Then
         cg.FrameBuffer(x, y) = P1
      End If
      y = yy2
      If (x >= 0) And (x < cg.Header.Nx) And (y >= 0) And (y < cg.Header.Ny) Then
         cg.FrameBuffer(x, y) = P1
      End If
   Next x
   For y = yy1 To yy2 '【左右枠線を指定色P1で描画】
      x = xx1
      If (x >= 0) And (x < cg.Header.Nx) And (y >= 0) And (y < cg.Header.Ny) Then
         cg.FrameBuffer(x, y) = P1
      End If
      x = xx2
      If (x >= 0) And (x < cg.Header.Nx) And (y >= 0) And (y < cg.Header.Ny) Then
         cg.FrameBuffer(x, y) = P1
      End If
   Next y
   For x = (xx1 + 1) To (xx2 - 1) '【枠内を指定色P2で塗り潰す】
      For y = (yy1 + 1) To (yy2 - 1)
         If (x >= 0) And (x < cg.Header.Nx) And (y >= 0) And (y < cg.Header.Ny) Then
            cg.FrameBuffer(x, y) = P2
         End If
      Next y
   Next x
End Sub

'###########################################################
'### イメージを領域コピーする
'###########################################################
'【機能】枠線を指定色P1で描画し、枠内を指定色P2で塗り潰す
Public Sub ImageCopy(dstcg As NTL_BitMap, dstx As Long, dsty As Long, srccg As NTL_BitMap, srcx0 As Long, srcy0 As Long, srcxw As Long, srcyw As Long, T As RGB24bitPixel)
   Dim Sx As Long, Sy As Long, dX As Long, dY As Long, P As RGB24bitPixel
   
   dY = dsty
   For Sy = srcy0 To srcy0 + srcyw - 1
      dX = dstx
      For Sx = srcx0 To srcx0 + srcxw - 1
         P = GetPixel(srccg, Sx, Sy)
         If (P.Red <> T.Red) Or (P.Green <> T.Green) Or (P.Blue <> T.Blue) Then
            DrawPixel dstcg, dX, dY, P
         End If
         dX = dX + 1
      Next Sx
      dY = dY + 1
   Next Sy
End Sub

【標準モジュール：NTL_Lib_FDTD_2D_basicのソース】

　２次元FDTDの主要演算部分です。
'###########################################################
'【2次元FDTD】単一物質(空気)
' 音圧格子(node)と粒子速度枝(branch)の接続関係
' Ny      v     v     v     v     v     v               v     v
'         |     |     |     |     |     |               |     |
' Ny-1 u--P--u--P--u--P--u--P--u--P--u--P-･･････････-u--P--u--P--u
'         |     |     |     |     |     |               |     |
' Ny-1    v     v     v     v     v     v               v     v
'         |     |     |     |     |     |               |     |
' Ny-2 u--P--u--P--u--P--u--P--u--P--u--P-･･････････-u--P--u--P--u
'         |     |     |     |     |     |               |     |
' Ny-2    v     v     v     v     v     v               v     v
'         |     |     |     |     |     |               |     |
'    :    :     :     :     :     :     :               :     :
'         |     |     |     |     |     |               |     |
'    2 u--P--u--P--u--P--u--P--u--P--u--P-･･････････-u--P--u--P--u
'         |     |     |     |     |     |               |     |
'    2    v     v     v     v     v     v               v     v
'         |     |     |     |     |     |               |     |
'    1 u--P--u--P--u--P--u--P--u--P--u--P-･･････････-u--P--u--P--u
'         |     |     |     |     |     |               |     |
'    1    v     v     v     v     v     v               v     v
'         |     |     |     |     |     |               |     |
'    0 u--P--u--P--u--P--u--P--u--P--u--P-･･････････-u--P--u--P--u
'         |     |     |     |     |     |               |     |
'y   0    v     v     v     v     v     v               v     v
'↑    0  0  1  1  2  2  3  3  4  4  5  5 ･･････････ Nx-2  Nx-1  Nx
'｜
'＋--->x
' 左右両端のuは完全吸収端とする
' 上下両端のvは更新せず、音源として駆動しない場合は v=0 とする
'###########################################################
Public Type FDTD_2D
   flag As Byte   ' 演算条件ﾌﾗｸﾞ
   dt As Double   ' 時間分解能[秒]
   dL As Double   ' 空間分解能[ｍ]
   Nx As Long     ' 音圧格子のx方向の個数
   Ny As Long     ' 音圧格子のy方向の個数
   m() As Byte    ' 音圧格子単位の物質番号
   P() As Double  ' 音圧成分
   u() As Double  ' 粒子速度x成分
   v() As Double  ' 粒子速度y成分
   A As Double    ' A=k*dt/dL (体積弾性率k=ρc^2)
   B As Double    ' B=dt/(ρ*dL)
End Type

'###########################################################
'### 単一物質2次元FDTDｼﾐｭﾚｰｼｮﾝ処理 【初期化処理】
'###########################################################
Public Sub Initialize_FDTD_2D(Model As FDTD_2D, MaterialName As String, dt As Double, dL As Double, Nx As Long, Ny As Long, flag As Byte)
   Dim i As Long, j As Long, M1 As Material ' 媒質の物質定数
   
   M1 = MaterialInformation(MaterialName, 20#)  ' 室温20℃での媒質密度ρ、音速、体積弾性率
   
   With Model ' FDTD2次元ﾓﾃﾞﾙ
      .dt = dt
      .dL = dL
      .Nx = Nx
      .Ny = Ny
      .A = M1.VolumeElasticity * dt / dL
      .B = dt / (M1.Density * dL)
      ReDim .P(.Nx - 1, .Ny - 1), .m(.Nx - 1, .Ny - 1), .u(.Nx, .Ny - 1), .v(.Nx - 1, .Ny)
      For j = 0 To .Ny - 1
         For i = 0 To .Nx - 1: .P(i, j) = 0#: .u(i, j) = 0#: .v(i, j) = 0#: Next i
      Next j
      For j = 0 To .Ny - 1: .u(.Nx, j) = 0#: Next j
      For i = 0 To .Nx - 1: .v(i, .Ny) = 0#: Next i
      .flag = flag
   End With
End Sub

'###########################################################
'### 単一物質2次元FDTDｼﾐｭﾚｰｼｮﾝ処理
'###########################################################
'【2次元FDTD】単一物質(空気)
Public Sub Cal_FDTD_2D(Model As FDTD_2D)
   Dim i As Long, j As Long
   
   With Model
      '【音圧格子Pの計算】
      For j = 0 To .Ny - 1
         For i = 0 To .Nx - 1
            .P(i, j) = .P(i, j) - .A * (.u(i + 1, j) - .u(i, j) + .v(i, j + 1) - .v(i, j))
         Next i
      Next j
      
      '【粒子速度格子u】
      For j = 0 To .Ny - 1
         For i = 1 To .Nx - 1
            .u(i, j) = .u(i, j) - .B * (.P(i, j) - .P(i - 1, j))
         Next i
      Next j
      '【粒子速度格子v】
      For j = 1 To .Ny - 1
         For i = 0 To .Nx - 1
            .v(i, j) = .v(i, j) - .B * (.P(i, j) - .P(i, j - 1))
         Next i
      Next j
   
      '【音圧格子が物体のときは固定端反射】
      For j = 0 To .Ny - 1
         For i = 0 To .Nx - 1
            If .m(i, j) <> 0 Then
               .v(i, j) = 0#
               .v(i, j + 1) = 0#
               .u(i, j) = 0#
               .u(i + 1, j) = 0#
            End If
         Next i
      Next j
   
   End With
End Sub

【標準モジュール：NTL_Lib_FDTD_definitionのソース】

　媒質データの定義や音源駆動関数のソース部分です。
Public Type SoundSource
   Freq As Double
   jTau As Long
   Omega As Double
   nWaveNum As Double
End Type

Public SoundInfo As SoundSource

Public Type Material
' 温度20℃の時、空気の体積密度ρは1.205kg/m3、音速ｃは343m/秒。
   MaterialName As String      ' 物質名
   Density  As Double          ' 媒質密度ρ[kg/ｍ3]
   Velocity As Double          ' 音速    ｃ[ｍ/秒]
   VolumeElasticity  As Double ' 体積弾性率＝ρｃ^2
   Absorption As Double        ' 吸収係数
End Type

'###########################################################
'### 物質定数
'###########################################################
Public Function MaterialInformation(MaterialName As String, Temperature As Double) As Material
   With MaterialInformation
      .MaterialName = MaterialName
   Select Case MaterialName
      Case "アンモニア", "NH3"
         .Velocity = 415# + 0.73 * Temperature
         .Density = 0.771
      Case "一酸化炭素", "CO"
         .Velocity = 337# + 0.604 * Temperature
         .Density = 1.2504
      Case "13A" ' 都市ガス
         .Velocity = 376.29 + 0.6323 * Temperature
         .Density = 1.2504
      Case "塩素", "Cl2"
         .Velocity = 205.3  ' 温度項の記述無し
         .Density = 3.214
      Case "空気", "Air"
         .Velocity = 331.45 + 0.607 * Temperature
         .Density = 1.205
      Case "酸素", "O2"
         .Velocity = 317.2 + 0.57 * Temperature
         .Density = 1.429
      Case "水蒸気"
         .Velocity = 404.8  ' 温度項の記述無し
         .Density = 0.598
      Case "水素", "H2"
         .Velocity = 1269.5 + 2# * Temperature
         .Density = 0.08988
      Case "窒素", "N2"
         .Velocity = 337# + 0.85 * Temperature
         .Density = 1.25055
      Case "二酸化炭素", "CO2"
         .Velocity = 258# + 0.87 * Temperature  ' (低周波)
     '   .Velocity = 268.6# + 0.87 * Temperature ' (低周波)
         .Density = 1.9769
      Case "エチルアルコール"
         .Velocity = 1207#   ' 温度項の記述無し
         .Density = 786#
      Case "クロロホルム"
         .Velocity = 995#   ' 温度項の記述無し
         .Density = 1490#
      Case "グリセリン"
         .Velocity = 1986#   ' 温度項の記述無し
         .Density = 1260#
      Case "水銀", "Hg"
         .Velocity = 1450#   ' 温度項の記述無し
         .Density = 13600#
      Case "水", "H2O"
         .Velocity = 1500#   ' 温度項の記述無し
         .Density = 1000#
      Case "海水"
         .Velocity = 1513#   ' 温度項の記述無し
         .Density = 1021#
      Case "ベンゼン"
         .Velocity = 1295#   ' 温度項の記述無し
         .Density = 870#
      Case "アルミニウム", "Al"
         .Velocity = 6420#   ' 温度項の記述無し
         .Density = 2690#
      Case "エボナイト"
         .Velocity = 2500#   ' 温度項の記述無し
         .Density = 1200#
      Case "黄銅"
         .Velocity = 4700#   ' 温度項の記述無し
         .Density = 8600#
      Case "ガラス"
         .Velocity = 5440#   ' 温度項の記述無し
         .Density = 2420#
      Case "金", "Au"
         .Velocity = 3240#   ' 温度項の記述無し
         .Density = 19320#
      Case "銀", "Ag"
         .Velocity = 3650#   ' 温度項の記述無し
         .Density = 10490#
      Case "氷"
         .Velocity = 3230#   ' 温度項の記述無し
         .Density = 917#
      Case "コンクリート"
         .Velocity = 4250#    ' 5250#  ' 温度項の記述無し
         .Density = 2200# '資料には記述無し
      Case "ステンレス"
         .Velocity = 5790#   ' 温度項の記述無し
         .Density = 7910#
      Case "大理石"
         .Velocity = 6100#   ' 温度項の記述無し
         .Density = 2650#
      Case "チタン", "Ti"
         .Velocity = 5990#   ' 温度項の記述無し
         .Density = 4580#
      Case "鉄", "Fe"
         .Velocity = 5950#   ' 温度項の記述無し
         .Density = 7860#
      Case "銅", "Cu"
         .Velocity = 5010#   ' 温度項の記述無し
         .Density = 8960#
      Case "鉛", "Pb"
         .Velocity = 1960#   ' 温度項の記述無し
         .Density = 11340#
      Case "ニッケル", "Ni"
         .Velocity = 6040#   ' 温度項の記述無し
         .Density = 8900#
   End Select
   .VolumeElasticity = .Density * .Velocity * .Velocity ' 体積弾性率k=ρc^2
   .Absorption = 0.47
   End With
End Function

'###########################################################
'### 音源駆動関数(粒子速度) 【初期化処理】
'###########################################################
Public Sub Initialize_SoundSource(nWaveLength As Long, dt As Double, nWaveNum As Double)
   Dim Pi As Double
' nWaveLength 1周期分のｻﾝﾌﾟﾙ数
' nWaveNum    駆動波数
   Pi = Atn(1#) * 4#
   With SoundInfo
      .nWaveNum = nWaveNum
      .Freq = 1# / (nWaveLength * dt)
      .jTau = nWaveLength * nWaveNum
      .Omega = 2# * Pi * .Freq * dt
   End With
End Sub

'###########################################################
'### 音源駆動関数(粒子速度)
'###########################################################
Public Function SoundSource(jTim As Long) As Double
   If jTim < SoundInfo.jTau Then
      SoundSource = Sin(SoundInfo.Omega * jTim) * (0.5 - 0.5 * Cos(SoundInfo.Omega * jTim / CDbl(SoundInfo.nWaveNum)))
   Else
      SoundSource = 0#
   End If
End Function

【Sheet2(２次元音響FDTD)のソース】


Private MyPic As NTL_BitMap            ' RGB24bit形式のBITMAP構造体(Colorﾃｰﾌﾞﾙも含む)
Private MyScale As NTL_BitMap          ' カラースケールのBitMap構造体を用意する
Private FileName As String             ' BITMAPﾌｧｲﾙの名称
Private rFileName As String            ' 結果保存用ﾌｧｲﾙ名
Private Const dt As Double = 0.00001   ' 時間ピッチ100ns
Private Const dL As Double = 0.05      '
Private Const Nx As Long = 200         ' FDTD演算結果表示用画像のx方向画素数
Private Const Ny As Long = 200         ' FDTD演算結果表示用画像のy方向画素数
Private Const nWaveWidth As Long = 200 ' 駆動音波の波長
Private Const nWave As Double = 0.5    ' 音源駆動波数
Private flag As Byte                   '
Private Model As FDTD_2D               '
Private Const nColor As Long = 180     ' ｶﾗｰﾊﾟﾚｯﾄのx方向画素数(色数)
Private Const nBarWidth As Long = 10   ' ｶﾗｰﾊﾟﾚｯﾄのy方向画素数(Barの表示幅)
Private Const Sx As Long = 85          ' 音源x座標
Private Const Sy As Long = 120         ' 音源y座標
Private angle As Integer
Private MaxValue As Double

Private Sub CommandButton1_Click()
   Dim i As Long, wk As Double, jTim As Long, ix As Long, iy As Long, L As Long, code As Long

'【ｶﾗｰﾊﾞｰ】
   BuildBitMap nColor, nBarWidth, MyPic ' ｶﾗｰﾊﾟﾚｯﾄ用のRGB24bit形式BITMAP(ｶﾗｰﾊﾞｰ)を作成する
   FillBitMapImage MyPic, MyPic.White   ' ｶﾗｰﾊﾞｰを指定色(白)で初期化する
   CreateMonopoleScale MyPic, nColor    ' 絶対値表示用のｶﾗｰﾊﾟﾚｯﾄを作成する
'  CreateDipoleScale MyPic, nColor      ' ±表示用のｶﾗｰﾊﾟﾚｯﾄを作成する
   For i = 0 To nColor - 1              '
      DrawLine MyPic, i, 0, i, nBarWidth - 1, MyPic.CLUT(i) '
   Next i                               '
   FileName = ThisWorkbook.Path + "\ColorBar" + Format(nColor) + "×" + Format(nBarWidth) + ".bmp"
   CreateBitMapFile MyPic, FileName       ' 指定名称で画像ﾌｧｲﾙ(180×10画素)を作成する
   Image2.Picture = LoadPicture(FileName) ' 指定名称の画像ﾌｧｲﾙを読み込んで表示する
   GetBitMapFile MyScale, FileName        ' スケールを読込む
   InitializeColorLookupTable MyScale, nColor
   CreateMonopoleScale MyScale, nColor    ' 絶対値表示用のｶﾗｰﾊﾟﾚｯﾄを作成する
   DrawLine MyScale, 0, 0, 179, 0, MyScale.DarkGray
   DrawLine MyScale, 179, 0, 179, 9, MyScale.DarkGray
   DrawLine MyScale, 179, 9, 0, 9, MyScale.DarkGray
   DrawLine MyScale, 0, 9, 0, 0, MyScale.DarkGray

'【FDTD演算結果表示部】
   SavePicture Image5.Picture, ThisWorkbook.Path + "\Temp000.bmp"
   GetBitMapFile MyPic, ThisWorkbook.Path + "\Temp000.bmp"
   FileName = ThisWorkbook.Path + "\BackGround" + Format(Nx) + "×" + Format(Ny) + ".bmp"
   CreateBitMapFile MyPic, FileName       ' 指定名称で画像ﾌｧｲﾙを作成する
   Image1.Picture = LoadPicture(FileName) ' 指定名称の画像ﾌｧｲﾙを読み込んで表示する

   Initialize_FDTD_2D Model, "空気", dt, dL, Nx, Ny, 0  '１[μ秒],１[mm]
   
   For iy = 0 To Ny - 1
      For ix = 0 To Nx - 1
         code = MyPic.FrameBuffer(ix, iy).Red * &H10000 + MyPic.FrameBuffer(ix, iy).Green * &H100& + MyPic.FrameBuffer(ix, iy).Blue
         If code = &HFFFFFF Then Model.m(ix, iy) = 0 Else Model.m(ix, iy) = 7
      Next ix
   Next iy
   Initialize_SoundSource nWaveWidth, dt, nWave

   MaxValue = 0
   flag = 1
   For jTim = 0 To 1500
      wk = SoundSource(jTim) '係数がないとオーバーフローするため
      Model.P(Sx, Sy) = wk   '【音源座標】ここでは座標(Sx,Sy)を『音圧』ソースと定義しています。
      Cal_FDTD_2D Model
      Label3.Caption = "Time Step = " + Format(jTim)
      
      If jTim Mod 5 = 0 Then '【表示間隔】５ステップ毎
         Image1.Picture = LoadPicture(FileName)
         DisplayImageAbsolute '【振幅表示に適したカラースケール】
'        DisplayImageValue
         ImageCopy MyPic, 15, 2, MyScale, 0, 0, 180, 10, MyScale.Brown
         CreateBitMapFile MyPic, ThisWorkbook.Path + "\Temp.bmp"
         Image1.Picture = LoadPicture(ThisWorkbook.Path + "\Temp.bmp")
         If jTim Mod 10 = 0 Then '【画像ファイル作成間隔】１０ステップ毎
            rFileName = ThisWorkbook.Path + "\FDTD2D_Result" + Format(jTim, "0000") + ".bmp"
            CreateBitMapFile MyPic, rFileName
         End If
      End If
      DoEvents
      If flag = 0 Then Exit Sub
   Next jTim
End Sub

Private Sub CommandButton2_Click()
   flag = 0 '【FDTDｼﾐｭﾚｰｼｮﾝ中断フラグ】
End Sub

'***********************************************************
'*** カレントイメージに音圧値を描画する
'***********************************************************
Private Sub DisplayImageAbsolute()
   Dim i As Long, j As Long, v As Double, vmax As Double, k As Long, L As Long

   vmax = -1E+100
   With Model
      For j = 0 To .Ny - 1
        For i = 0 To .Nx - 1
           v = Abs(.P(i, j))
           If v > vmax Then vmax = v
        Next i
      Next j
      If MaxValue < vmax Then MaxValue = vmax
      Sheet2.Cells(3, 8) = MaxValue '【過去最大値】
      Sheet2.Cells(4, 8) = vmax     '【現時刻最大値】
      If vmax < 0.15 Then vmax = 0.15 '固定表示
      If vmax = 0# Then
         vmax = 1#
      Else
         vmax = (nColor - 0.55) / vmax
      End If
      For j = 0 To .Ny - 1
         For i = 0 To .Nx - 1
            k = CLng(Abs(.P(i, j)) * vmax)
            If .m(i, j) = 0 Then DrawPixel MyPic, i, j, MyPic.CLUT(k)
         Next i
      Next j
   End With
End Sub

Private Sub DisplayImageValue()
   Dim i As Long, j As Long, v As Double, vmax As Double, vmin As Double, k As Long, L As Long

   v = 149 / 180 ' 180を149色に割り当てる
   With Model
      For j = 0 To .Ny - 1
         For i = 0 To .Nx - 1
            k = CLng(.P(i, j) * v + 150)
            If k >= nColor Then k = nColor - 1
            If k < 0 Then k = 0
            If .m(i, j) = 0 Then DrawPixel MyPic, i, j, MyPic.Color(k)
         Next i
      Next j
   End With
End Sub

Private Sub Image3_Click()
   Image5.Picture = Image3.Picture
End Sub

Private Sub Image4_Click()
   Image5.Picture = Image4.Picture
End Sub

Private Sub Image6_Click()
   Image5.Picture = Image6.Picture
End Sub

Private Sub Image7_Click()
   Image5.Picture = Image7.Picture
End Sub

【追加情報】【WebGLのページ】

JavaScriptでリアルタイム３次元グラフィックスが楽しめる時代の到来までもうすぐです。
現在、Google Chrome12、Mozilla FireFox5、及び開発版のSafari、Operaで使うことができます。
WebGLは、JavaScriptエンジンv8とDirectXを利用して高速化したOpenGLを融合して構成されたものです。
2011年5月に指摘されたセキュリティー問題はありますが、処理速度の高速性から利便性も高いので、将来的には普及するのではないでしょうか？
　●WebGLとは？　●事前準備：利用できるウェブブラウザは？　●事前準備：FireFox5でWebGLが動作しない場合には？
　●WebGLのサンプル　●WebGLのデモ　●クロスドメインテクスチャーを将来も利用できる方法