独り言日記（2006/06）

独り言日記

GPUで交差判定を試してみました（2006/06/30）

定番（？）ってことで、GPUで三角形とレイの交差判定を実験。自動的に並列化されて速くなるのでは、と勝手に期待していたのですが、結果はPentium4 2.8GHzマシンとGeForce6600GTボードでチェックしてみたのですが、大差なし。う〜ん、残念。

D3DFMT_A32B32G32R32FのRGBAそれぞれfloatタイプのテクスチャを作成して、3枚に三角形の頂点座標をバラにして格納、TARGET（出力先）として同じくfloat*4のテクスチャに計算後のレイから交点までの距離t、三角形での内分点情報U/Vを格納、としたのですが・・・。たとえば256x256ピクセルのテクスチャだとすると、65536ポリゴンとレイの交差判定がGPU任せで1回の描画で求まることになります。三角形とレイの交差判定処理は、そのままピクセルシェーダーにのせることができました（ただ、なんか挙動不審な予感）。

PC側のメモリからGPU管理下のメモリへのブロック転送はボトルネックになるかと思ったのですが（たしかになりはしますが）、まだマシなほう。むしろ、ピクセルシェーダーでのレンダリング自身がボトルネック化してます。後、除算・ifによる分岐は、やはり足を引っ張りますね。

ただ、GPUに放り投げればそれなりの速度で交差判定計算をしてくれてはいますので、CPU/GPUにてスレッドで並列処理できるとしたら、これは速度アップにつながるかな？それよりも、根本のアルゴリズムを推敲しないといかんでしょうね(^_^;;。まだまだ私自身修行が必要なようです。

ソフトウェア開発環境展（2006/06/28）

東京ビッグサイトで開催中です。で、行ってまいりました。

http://www.sodec.jp/SODEC/

これは、ソフトウェア・セキュリティ・組み込みなど、各社が製品を展示・商談するためのイベントです。何社か、詳しく製品についての話を聞いてきました。

去年も行ってきたのですが、今年はソフトウェアでは特に目立ったのはなかった気が。心なしか、規模も縮小した？ソフトウェアに関しては飽和してるな、というのが感じられました。後、目立って見えた点としてEclipse上の開発も増えてましたね。ただ、ここに参加していない会社で面白い技術を持ったところはさまざまありますので（Googleとか）、あくまでも基幹系として見たとすると、次の新しいネタを考えないといけないかなぁと思ったりします。CRMって技術ではなくて、あくまでも囲い込みの手段でしかないかな、というのが正直な感想です。それをかっこよく３文字英語にしました、というのはそう長くは続かないのでは・・・と思ってます（個人的な意見です）。というよりも、ごく当たり前のことをしているのがCRMなんですよね。むしろ、環境を提供するもの、ってのはネットが普及しきった現在では、そろそろ「便利になった」が当たり前になって、それ以上のものは難しいんでないかなぁと。ネットワーク（インターネット）の普及レベルのビッグバンがないと、技術的なことや便利さに関しては、ネタ探しに苦労しそうです。

セキュリティは詳しくは分からないのですが、Winny関連を出しているところが増えましたね。後、ブースも増えたかな。指紋認証とかその他ハードがらみのがあって、セキュリティに関してはこれからのびそうですね。

個人的には、組み込み系が面白かったです。FPGAを各社押し出しているのが目立ってました。Xilinx社もPCI-Express対応を大きく前に出してましたねぇ、後、10万以下のボード（個人向けか）を紹介してる会社も見受けられました。私が過去実験してたのもあって、ボトルネックについてとあるブースで話をしていたのですが、このへん、意見が食い違っていたようです。正直、２Ｄ画像処理ですと、GPUでできてしまう時代になってます（ただ、FPGAでの画像処理を出展してるブースも多々ありました）。なんで、もう一歩進むとすると、汎用的な並列演算機としてのFPGAが紹介されていれば、なんて思ったり（なんとなくですが、客層がFPGAを触りたいユーザではなくて、それで実現したシステムを使いたい、がメインだからか、数値比較みたいなのはなかったです）。ただ、HDLがC言語ライクになってたり、XilinxのようにPCI-Express対応をうたってたり、と、進化は確実にしていってますね。私は純粋にハードウェアとしてソフトウェア的にプログラムができる、っていう点に非常に好かれるものがあるのですが、第三者の視点で見たらどうなんでしょうねぇ。利用価値という点では、まだまだ未開拓な分野であると思ってます。まっさらな「白紙」がFPGAであるとすると、それに乗っける機能は自由に作れるわけなんで、その先はアイデア次第ですかね。時間があればFPGAはもっと掘ってみたいところです（今は、遊びでGPUのほうをいじってますが、なんとなくGPUで「なぜこんな制限があるんだろう」ってのが分かってきた気がします。もしかしたら、FPGAみたいな使い方もできるかな）。

Light Space Perspective Shadow Maps（2006/06/26）

シャドウマップの方法として、「Light Space Perspective Shadow Maps」というのがあります（略して「LiSPSM」。「LSPSM」と書いているサイトもあります）。

御本家。論文のあるサイト（サンプルソースもあります）

http://www.cg.tuwien.ac.at/research/vr/lispsm/

LiSPSMについてのゲームでの解説

http://watch.impress.co.jp/game%2Fdocs/20050929/3dinis.htm

日本のサイトでは、ぐっと深く突っ込んだサイトがなかなかなかったので、いろいろ試行錯誤した結果を。まず、御本家の論文およびサンプルソースを読むと、何をしているのか、というのは理解はできると思います。特に、サンプルソースはUniformなシャドウマップ（一番簡単なシンプルシャドウマップ）、LiSPSMの両方の実装が書いてますので、見比べながらだと分かりやすいかと。

シャドウマップの基本はどれも同じで、

光源から見たときのオブジェクトを光源を中心にレンダリング
  （光源ビューでの深度を得る）。シャドウマップを生成する。
視点から見たときの影を求めたい交点を光源座標に変換し、
  シャドウマップに格納された深度よりも大きな値が入っている場合は影と判断。

となります。２パスの考え方です。

Uniformなシャドウマップでは（以下USM）、単純に上記を行っています。ただ、欠点としてはシーンが大きくなった場合に、視点に近い部分の影がガタガタになる、というのがあります。

このときのシャドウマップ（512x512 pixel）は以下。

これは、視点から見たときの近遠も関係なく「（視点を考慮しない）光源から見たシャドウマップ」を生成しているからでもあります。

これを改善したものとしてPSMがあります（すみません、これはアルゴリズムを掘ってません(^_^;;）。で、飛んでLiSPSMでは、これにPerspective補正＋αがかかっているものになります。

このときのシャドウマップ（512x512 pixel）は以下。

シャドウマップ自身にパースペクティブがかかっている・光源座標で新たに変換行列を構築している、以外にも、視点から見たときのビューの四角錐、およびシーン全体のAABB（バウンディングボックス）、そして光線ベクトルでクリッピングしているという点が大きな最適化のように感じます。それらはサンプルソースに書いてますね。非常に参考になります。御本家の論文サイトのサンプルは平行光源の場合の実装になります。

点光源の場合はキューブマップ的に、スポットライトの場合は範囲を調整して、とするといけそうですが、平行光源のシャドウマップよりも負荷はかかりそうですね。

LiSPSMは、シャドウマップに毛が生えたもの程度にしか考えてなかったのですが、なかなか奥が深い・・・。結構理解するまでに時間がかかりました（実際はシャドウマップ以外の部分（クリッピングとか）が肝かなぁと思ったりします）。と、今回はさわりだけの解説でした。

ASUSのPhysXボード（2006/06/25）

ASUSびいきな私としては、ちょっと店頭で拝んできたかったのですが、置いてるところを見つけられませんでした。

http://plusd.itmedia.co.jp/pcuser/articles/0606/25/news001_3.html

まだ4万円もするのね・・・。たぶん物理演算のハードウェア化の普及やバス幅（PCI==>PCI-Express）も過渡期だと思うので、見送ってる人もいるのかなぁとか思ったりします（私もそうです）。

この記事で「世の中には純粋に物理演算を楽しみたい人、がいるらしい」って・・・。いやそうなんですが、マニア向けなのかな(^_^;;。ワンチップMSXのFPGAもマニア向けだったのか結局発売までこぎつけられませんでしたね。個人的にはボードとPhysXSDKだけでかなり遊び甲斐があると思うのですが・・・。しかし、果たして普及するかどうか、ある意味楽しみではあります。

GeForce6600GT（2006/06/24）

GeForce6600GTのボードを秋葉で買ってきました。いまさらかよ、って感じですが(^_^;;。自宅のメインマシンはノートPCなので使えない、ということで仕事場のに入れてます。しかし、ひさびさに秋葉に行ってグラフィックカードを見ていたのですが、AGP対応のがえらい少なくなりましたね。もう、時代はPCI-Expressか・・・。マザーボードのほうも、AGPをサポートしていないのも増えてきてました。私の使っているマシンは、ノーマルなAGP+PCIなので古いのを探すしかありませんでした（わざわざ高いGeForce7系を買うつもりもなかったですし、最低限の用件を満たすので満足）。

で、6600GTです。ようやくVertex/PixelShader3.0にネイティブ対応、そして浮動小数点テクスチャが使える・・・。シャドウマップの実験をしてみたいのと、ちょっとGPUレイトレのネタが思いつきまして、これが実際にハードウェアに乗るかどうかテストしたいと思ってます（でも、例によってメモリ転送でもたつくかなぁ）。これがうまいこといくと、汎用的なレイトレからGI（パストレとか）にも使えるかもしれない、と思ってるのですが・・・さて。しかし、HLSL触ってたのはもう4年前か・・・えらいブランクがあきました。

Googleで「HLSL」で検索すると、なんかうちの過去のサイトがMSの次に出てきてました。なつかしっ(^_^;;！！

あれからずいぶん進化してますね＜VS3.0/PS3.0

小さいことからコツコツと（2006/06/18）

とはよくいったものですね、ほんとにそう思います。ということで、仕事の小さなひと山（技術的な部分）を１つクリア。早速次の山に取り掛からねば。今年はそんな調子で、考えては実装して、の繰り返しです。何回アルゴリズムを練ったことやら。何かに書き留めておかないと忘れてはダメなので、メモを残すようにはしてますが、その管理自体が散漫なので、記録したこと自体忘れてしまいそうです（汗）。ノートはすでに何冊かなぁ、20冊くらいはたまったかな。大事なネタ帳だったりします。しかし、持ち運びがそろそろ大変になったので早めにデジタル化したいところですが・・・

しかし、ワールドカップがう〜んな結果ですね。がんばれ日本、ですが、、2002年と比べると・・・。じっくり観戦したいんだけど、なかなか仕事しながらだとどっちも集中できない(^_^;;。

物理から見る弾道計算（2006/06/18）

今話題なので・・・（苦笑）。あくまでも物理・数学の勉強ということで。このサイトが非常に詳しいですね。

http://homepage3.nifty.com/kubota01/

後、数学・物理がよく分かってないプログラマ（含む、私）のためか、ソースコードも公開されてます。・・・わかってらっしゃる！！一般のプログラマは、数式をプログラムで示してくれると一番よく分かったりするもんなんです。

この分野はよく分からなかったりしますが、数学・物理による卓上計算が一番実践的に使われている分野の１つであると感じますね。さて、例のミサイルの行方はどうなるのだろう・・・18日あたりが焦点になってますが・・・（願わくは、打ち上げ中止か、前のように日本を飛び越えて海に着水、となってほしいところです）

ところで、GoogleEarthをダウンロードしました。

http://earth.google.com/download-earth.html

こいつは面白い！！はじめに地球儀が出て、回転したりどんどん拡大していってその地域の衛星または航空写真が拝めます。写真を撮ったメーカー（？）が、場所によって変わりますので世界的に情報を集めたんでしょうかね。場所によっては画像が荒かったりします。東京と大阪を拡大してみると、東京は鮮明だけど大阪は粗い(^_^;;。

地理が苦手だったので（未だに日本の都道府県を覚えてない（汗））、これ見ながらだと覚えれそうな気がします(^_^;;。ちっこい島を拡大して地名が出てきたりだとか、RPGみたいな発見がありますね。へぇ、沖縄ってこんなに遠い位置にあったのか、とか（＜無知丸出し）

IntelMacでのUniversalBinary対応 in コマンドライン（2006/06/17）

アップルのサイトでは、必要な情報が探しづらいですねぇ(^_^;;。デベロッパ向けのUB関連のページをようやく再発掘。

http://developer.apple.com/jp/documentation/MacOSX/Conceptual/universal_binary/index.html

XCodeでは、プロジェクトのオプションを変更することでPPC/Intel両方に対応するUB対応のアプリケーションが組めるのですが（上記URLにやり方書いてます）、コマンドラインでビルドしたい場合はどうだろう？

調べてみました。上記サイトにも載ってはいるのですが、「-arch ppc -arch i386」をgcc(g++)オプションに加えるといいです。ただ、-archオプションは「-E/-M」関連のオプションとかぶるため、既存の./configureで作成されるMakefileではコンパイルがうまくいかない、ということもどうやら存在する模様。

以下、自分確認用メモ。UB対応のビルドを行うためのMakefile例です。ここでは静的ライブラリを作成しています。

CC = g++
CINCLUDE = -Iinclude -Iinclude/ogl
CFLAGS   = -DPARAM1 -DPARAM2 -O2 -fomit-frame-pointer -ffast-math -arch ppc -arch i386

OBJS = src01.o src02.o
TARGET = libxxx.a
all: $(OBJS) $(TARGET)

src01.o: src01.cpp
   $(CC) $(CFLAGS) $(CINCLUDE) -c -o src01.o \
 src01.cpp

src02.o: src02.cpp
   $(CC) $(CFLAGS) $(CINCLUDE) -c -o src02.o \
 src02.cpp

$(TARGET): $(OBJS)
   ar cru $(TARGET) $(OBJS)
   ranlib $(TARGET)

「CINCLUDE」にて、参照するインクルードディレクトリを指定してます。この場合は、「./include」と「./include/ogl」を参照。「CFLAGS」はコンパイラオプションで、「-D」に続く「PARAM1/PARAM2」はプリプロセッサ定義（#define）、「-O2」で最適化できるだけ最大に、「-f」は最適化のオプション、「-arch ppc -arch i386」でPPCとIntelつまりUB対応、となります。

後は普通のMakefileですが、OSX含むUNIX系の場合は、

xxx.o: xxx.cpp

の次の行での実行命令は、TABで空白をあける必要があります（BCCでは、スペースでした）。

最後の、「ar」と「ranlib」命令はライブラリを生成する命令です。

でうまくビルドが成功すると、単純にUBでないバイナリと比べると２倍くらいのサイズになりますね。PPCとIntelのバイナリがちゃんぽんされた状態です。

続・ODEでの空気抵抗（2006/06/13）

以前、ODEでの空気抵抗対応がなんかうまいこといかない、とか書いてたのですが、速度と位置を数値検証していたらどうにかいけるようになりました。

ODEでは「dBodySetLinearVel」にてリニアな（直線的な）速度ベクトルを得ることができるのですが、これの長さを速度vとして

R = (S * Cd * P * (v * v) / 2 - 6 * π * r * μ * v)

の式に入れてあげると空気抵抗値が求まります（ただし、この式で計算できるのは半径rの球体の場合のみ）。

「dBodySetLinearVel」で得たベクトルを(vx, vy, vz)としたときに

// ベクトルを正規化
d = sqrt(vx * vx + vy * vy + vz * vz);
vx = vx / d;
vy = vy / d;
vz = vz / d;

// 空気抵抗のベクトル
rx = -vx * R;
ry = -vy * R;
rz = -vz * R;

で求まった(rx, ry, rz)を「dBodyAddForce(bodyID, rx, ry, rz);」として入れると、リファレンス的にRunge-Kutta法で計算した空気抵抗ありの速度・位置とほぼ同じになりました。以前はODEでの力の計算を間違ってたかな。結局は力「F」を加えれば、意図したとおりに動作してますね。

後、PhysX SDKでも同じようにできそうですね。抵抗は空気抵抗だけでなくて「転がり抵抗」というのもあるらしいので、これも加えるとよりもっともらしくなるかも。抵抗としてもう１つ「摩擦」があるのですが、ODEの場合は静止摩擦だけのようです（動摩擦は対応していない模様）。PhysXでは静止・動摩擦両方を指定できますね。

関係ないけど、このサイトのWikiが重いですね(^_^;;。データ量も多いのもあって、そろそろ限界か・・・。そのままデータをPHP+MySQLのWikiとかに移行できるといいのですが・・・

PhysX SDKでの空気抵抗（2006/06/09）

AGEIA社のサポートページより、物理演算エンジンの「PhysX SDK」がダウンロードできます。

http://support.ageia.com/

ただし、先にユーザ登録する必要があります。このサイト自身のユーザ登録をすると、ダウンロードページにアクセスできるようになります。非商用で個人で使用する分には自由に利用することができます。

機能として魅力的だったので、ダウンロードしてテストしてみました。ODEよりも（ソフトウェアでの）計算が速いですね。これのドキュメントがよく出来てます。チュートリアルがDirectXのドキュメントみたいに充実してますね。さらっとしか見れてませんが、構成としてはODEと似ていて、シミュレーションで衝突などの「イベント」が起こった段階で駆動する仕組みっぽいです。

空気抵抗はPhysXでもサポートされてませんねぇ。サイトのKnowledge Baseにそのことが書いてました。空気抵抗計算では速度ベクトルの方向にてぶち当たる断面積の計算が必要となりますので、この部分をまじめに計算してしまうとかなりボトルネックになりそうではあります。基本の物理演算のさらに上位の機能、と位置づけているのかもしれませんね。

この空気抵抗ですが、対象オブジェクトの周囲のいくつかの方向からサンプリングして、3D==>2D圧縮して断面積を求める。そのサンプリング点を補間して周囲360度どの方向からの速度ベクトルが来ても、だいたいの断面積が求まる「関数」を作ったらたぶんリアルタイムでも複雑な形状での空気抵抗対応ができるんじゃないかな、と思ってます。空気抵抗こそ、SDK（とハードウェア）に組み込んでほしかったですねぇ。

単位系（2006/06/08）

最近は物理エンジンが面白くてよくいじってるのですが、完全に物理空間で制御された箱庭が出来れば面白そうだなぁ、と思ったりします。

で、3DCGソフトと物理の世界の差はなんだろう、と見ると「単位の違い」が結構目に付くように思います。今の3DCGの空間では、位置や時間の概念はあっても、質量や密度・速度はあまり関係していません（速度は存在するツールもあるかも）。

たとえばShadeだと（プラグインSDKでの）距離や位置の基本単位はmmです。物理の世界ではmが多いのかな。で、基本的な要素を比較したら以下のようになります。

---	Shade	物理空間
距離	mm	m
密度	---	kg/m^3
質量	---	kg
速度	---	m/s
時間	sec	sec
加速度	---	m/s^2
力	---	N

ほかにも電気が加わったりしたとすると、アンペアなどいろいろパラメータは考えられそうですね。ただ、単位といえども大きく「CGS単位系」「MKS単位系」「SI単位系」と３つが存在するように千差万別ではあります。SI単位系は国際単位と呼ばれ、いわゆる単位界のUTF-8みたいなもん（かな？）です。国際標準。

「CGS単位系」では、cm/g/secを基本としてます。「MKS単位系」はm/kg/secを基本としてます。「SI単位系」はm/kg/secを基本としてその他電気系のが統一されてたりするようです。

http://ja.wikipedia.org/wiki/SI%E5%8D%98%E4%BD%8D%E7%B3%BB

というわけで物理空間だとSI単位系で統一していくのが妥当かなと思ってます。しかし、いろいろサイトを回って物理の勉強をしてみると、単位がバラバラで結構迷いますね(^_^;;。分野によって偏りもあるようで。

ODE 0.6 RC2のソース修正分（2006/06/07）

「ODE 0.6(RC2)」の2点のソースに手を加えたものを置いてます。

ode_0_6_rc2_quickstep_util_src_20060607.zip(218)

ODEのソースを入れ替えてビルドし直すとOK。日本語のコメント入れてますので、どなたか、それを消してご本家を突っついていただければ幸いです（人任せ、すみません(^_^;;）。

「quickstep.cpp/util.cpp」のメモリ確保・解放部分にてヒープを使うように修正です。なお、これは0.6 RC2であることに注意してください。2006/06/07段階では、ご本家（ode.org）にてRC3が公開されています（でも、このスタックオーバーフローの修正はされていません）。

step処理部分で「dWorldQuickStep」を使うことで、今まで落ちていた部分がきっちりと動くようになってます。10000以上の衝突でも問題なしでした（修正前は1000くらいで落ちてました）。トライアングルメッシュ部（サンプルのtest_trimesh.cpp）でもdWorldStepでなくてdWorldQuickStepに変更、で上記修正（リビルド）済みのode.dllを使うことで落ちることはなくなってます。dWorldStep処理でも、やはりスタックオーバーフローが潜んでますので、同様にヒープに作業メモリを作成することで回避できるものと思われます（ただ、dWorldQuickStepを使うことに一本化するといじらなくても問題なしですが）。

これで、physXみたいに数千以上のボックスを積み上げて、そこに鉄球を落として崩壊させる、ってのができそうですね。

ODEのスタックオーバーフローバグ（2006/06/07）

ODEでのバグを追ってます。衝突数が1000を超えたあたりで落ちてしまう原因ですが、初歩的な部分でやっちゃってますね。

ODEのソース「quickstep.cpp」内のdxQuickStepper()が、dWorldQuickStep()関数から呼ばれる部分となります。ここで、

dRealAllocaArray(J, m * 12);

みたいな箇所が頻繁にあります。#defineを展開すると

dReal *J = (dReal*) ALLOCA ((m * 12)*sizeof(dReal));

となり、ALLOCAは16バイト境界でallocaするマクロです（ODE内で定義されてます）。

dReal *J = (dReal*) alloca (((m * 12)*sizeof(dReal) + 15) & 0xfffffff0);

な感じですね。sizeof(dReal)は8を返します。

で、「alloca」はスタックに対してメモリを確保し、freeしなくても自動解放するメモリ確保関数です。このときのmに14553とかが入ってきた場合、「(m * 12)*sizeof(dReal)」の計算で1397088バイト消費してしまいます。つまり、1MB近くをスタックにて消費です（VC++でのデフォルトのスタックサイズは1MBなんで、これでアウト）。数千の衝突でこのようにスタックオーバーフローですので、これをmallocでヒープに作業メモリを確保することで、たぶん落ちバグは回避できるはず。トライアングルメッシュ衝突時に落ちることがあるのもこれっぽい。スタックとヒープのメモリがうまいこと操れてないなぁと感じました。ただ、オープンソースの強みで ソースをいじって修正可能ですのでいじってみます。

だれか英語分かる方、ODEの作者様に報告してあげてほしいです(^_^;;。allocaを使っている箇所が随所にありますが、これだとすぐにあふれてしまいますよ〜〜。

ODEのdWorldQuickStep（2006/06/06）

自分用メモ。物理エンジン「ODE」での衝突計算を反復的に行うものとして「dWorldQuickStep()」があります。これは、「N x 衝突時のジョイントのアタッチ数」のメモリを消費します。Nは反復回数で「dWorldSetQuickStepNumIterations()」にて指定できます。デフォルトの反復回数は20です。このdWorldQuickStep関数は速度や位置を求める処理を近似しているといえます。

この「dWorldQuickStep()」は「dWorldStep()」関数のメモリ消費量と速度を改善した関数です。dWorldStepの場合は「衝突時のジョイントのアタッチ数 x 衝突時のジョイントのアタッチ数」のメモリを消費してしまうため現実的でないです。

ただし、「dWorldQuickStep()」はまだ大規模な処理には向かない感じです。衝突時のジョイントのアタッチ数（衝突の数）が1300くらいになるとODE0.6ではクラッシュしてしまいます。ODE0.5だともう少しアタッチ数が増えても問題なし。でも、1800くらいが限度っぽい。

数千の衝突は普通に存在しますので（というか、「接地」でも衝突判定はされるのでこれを入れてるとすぐにオーバーフロー）、ちょっとこのへんは実用では厳しいですね。ODEのソースも追っかけてみましたが、スタックオーバーフローしてる模様。メモリ管理の最適化が甘いのかなぁ。exe自身のメモリ消費量はさほど食ってないので、ローカルメモリで食ってしまってるような気が。ということで、ODEのソースをちとデバッグしてみよう（読むのが面倒な場合は独自でエンジン作ったほうが楽かも。要推敲です）。

ODEサンプルの「test_crash.cpp」はブロックを積み上げてそれを崩すサンプルですが、これ、ODE0.6ではそのままでエラーになってます（落ちます）。同様の問題と思われます（ODE0.5ではスタックギリギリで間に合ってるので問題なし）。

空気抵抗ありの世界での自由落下テスト（2006/06/06）

それでは、現実に近い話で。以下の条件で球体の自由落下を試みてみます。

球の半径       : 0.05 m (5cm)

重力加速度     : 9.8 m/s^2
空気抵抗係数   : 0.3
空気密度       : 1.292
空気粘性率     : 0.000018

水の質量密度   : 1000.0 kg/m^3
木材の質量密度 :  380.0 kg/m^3

空気抵抗係数や空気密度、空気の粘性率は温度などの環境によって変わりますので、上記定数はあくまで参考程度としてください。木材は、木の種類や水分の含有率によっても変わりますのでこれも参考程度。

上記条件にて、水と同じ質量の球体を高さ100m上空から初速度0で自由落下させたとします。計算結果は、約4.73秒後に38.517m/sの速度で地面に到達しました。

木材の球体で同様に自由落下させると、約5.1秒後に31.656m/sの速度で地面に到達しました。

ということで、同じ大きさの形状であっても 重い（質量が多い、密度が高い）物体ほど先に地面に到達することになります。

以下、今まで行った分のCの実験ソース（.NET2003のプロジェクト）です。DOS窓で動きます。空気抵抗なしのRunge-Kutta法にて速度を求める処理も参考としていれてます（これと、v=v0-gt / x = vo*t - 1/2*g*t^2、とで見比べたら分かりやすいかも）。

AirResistanceTest_src_20060606.zip(228) (29KB)

速度値をExcelなどでグラフ化したら、以前の日記のように終端速度に向かって曲線を描くのを確認できるかと思います。

で、話変わって物理演算エンジンのODEでは空気抵抗が考慮されていません。これに空気抵抗を実装するには、外力として「dBodyAddForce」で力を与えてあげる手がありますが・・・これってニュートン力のベクトル（空気抵抗値そのもの）を渡せばいいのだろうか？どうも外力で押したとしても質量を考慮した抵抗がかかるようで、思い通りに空気抵抗を加えられない・・・。「dBodySetLinearVel」で直接速度変更もできるのですが、これだとODEの計算と自前の速度計算が同時に入ってしまって上下に揺れる変な動きができてしまうし・・・。

Runge-Kutta法（2006/06/03）

さて、微分方程式を近似する手法「Runge-Kutta」についての説明です。

http://maya.phys.kyushu-u.ac.jp/~knomura/education/numerical-physics/text6/node2.html

が参考になります。

空気抵抗ありの最終的に物体（球体）にかかる力の式を以下であらわすとします。

m * (dv/dt) = (m * g) - (S * Cd * P * (v * v) / 2 - 6 * π * r * μ * v)
 
m  = 質量 (kg)
g  = 重力加速度 (m/s^2)
r  = 球の半径（m）
S  = 断面積（rであらわされる円の面積） (m^2)
Cd = 空気抵抗係数
P  = 空気密度 (kg/m^3)
π  = 円周率。3.1415926535...
μ  = 粘性率

とします。それぞれ定数または値はあらかじめ計算されてます。

これを、(dv/dt)の式にするために両辺をmで割ります。

(dv/dt) = g - (S * Cd * P * (v * v) / 2 - 6 * π * r * μ * v) / m

このときの式は微分方程式となり、右サイドにはvという変数が与えられた「(dv/dt) = f(v, t)」の関数と見ることができます。vは速度(m/s)、tは経過時間（秒）をあらわします。

この「(dv/dt) = f(v, t)」をRunge-Kutta法で近似できます。(dv/dt)は、単位時間における速度の差分、つまりは「加速度」になります。関数「f(v, t)」には加速度を求める計算を入れます。最終的に求めるのは速度vです。ここで初速度をvoとし以下の関数を作成しましょう。以下の「dt」は、tを細かく分割した十分に小さい値を指定します。tを30分割する（1秒間に30回の計算を行う）場合は、「dt = 1.0 / 30.0;」を指定します。

double calcVel(double v, double t) {
  double k1, k2, k3, k4;
  k1 = dt * f(v, t);
  k2 = dt * f(v + k1 * 0.5, t + dt * 0.5);
  k3 = dt * f(v + k2 * 0.5, t + dt * 0.5);
  k4 = dt * f(v + k3, t + dt);
  
  return (v + (k1 + 2.0 * k2 + 2.0 * k3 + k4) / 6.0);
}

実際はf(v, t)のtの引数は、加速度計算の式に含まれてませんので無視してもOKですね。これを使って「v = calcVel(v0, t);」と呼ぶことで、初速度voが与えられたときの次の加速度vが求まります。以降は、

v = calcVel(v, t);

と繰り返していくと、空気抵抗がある場合の自由落下の速度の変化を求めていくことができます。

速度はこれでOK。移動距離は速度vと時間dtが存在するので「x = v * dt」で算出できますね。このRunge-Kutta法は、速度の式だけでなくて、微分方程式にて精度を求めた近似を行いたい場合に結構有用に使えそうです。

精度はどんなもんか、というのを「v = vo + gt」のニュートンの式と、「F = ma」すなわち「m * (dv/dt) = m * g」つまり「(dv/dt) = g」をRunge-Kutta法で解いたものとを数値比較してみました（これをみても、空気抵抗がない場合は質量および落下の高さが関係ないのは分かりますね）。

100m上空から球体を初速度0で自由落下させた場合です。

経過時間(sec)	速度(v=vo + gt)	速度（Runge-Kutta法を使用）
0.0	0.0	0.0
1.0	9.8	9.8
2.0	19.6	19.6
3.0	29.4	29.4
4.0	39.2	39.2
5.0	49.0	49.0
6.0	58.8	58.8

0.000001くらいの誤差は丸めてます（まるめちゃ意味ないのですが(^_^;;）。で、結局同じ値です。ぶらぼ〜〜！ただし、分割の時間dtを十分小さい値をとらないと精度は変わるようです。上記は、dt = 1.0/30.0として、繰り返し計算したときの1.0秒ごとの結果を抜き出したものです。あと、加速度の式「(dv/dt) = f(v, t) = g」なもんでvが存在していないのも精度が高くみえる原因ではありますね。

では、次回はいよいよ「質量だけが異なる同じ大きさの球体を 同じ高さから自由落下させた場合に、どちらが先に地面に着地するか」の証明をしていこうかと思います。

空気抵抗ナシとアリのときの考察（2006/06/02）

もうちょっと深く掘ってみました。たとえば雨粒が落下する場合を考えます。

雨粒１つを「r = 半径 0.005m（5mm）、d = 質量密度1000kg/m^3」でできている球体とします。質量密度が1000kg/m^3ということは、これは「水」の密度をあらわしています。

この場合、半径と密度より 質量は球体の体積x密度で計算できます。

m = ((4.0 / 3.0) * π * r^3) * d

これを上空100mの高さから自由落下させたときの地面からの距離と速度を数値化します。

空気抵抗のない場合

ニュートンの式より、以下で速度と距離は求まります。

V = Vo + g * t;
X = Vo * t + 0.5 * g * t * t;

自由落下させるので初速度Voは0とし、鉛直下方向をプラスとしてます。

ただ、後で空気抵抗も考慮する都合上、

m * (dv/dt) = m * g
(dx/dt) = V

の微分方程式より算出することにします。このとき「Runge-Kutta法」を使用しました（それらの解説についてはまた後で）。これ、近似で展開できる方法のようです。(dv/dt)は加速度ですが、これより速度を求めることができます。同様に(dx/dt)は速度ですが、これより移動距離を求めることができます。

ニュートンの式で分かるとおり、時間tに比例して速度は増加します。線形になってますね。この速度をグラフ化しました。

縦方向が速度値、横方向が1/30秒ごとの時間です（+1で(1/30秒経過)）。

空気抵抗のある場合

前日に書いた「 m・(dv/dt) = (m x g) - (Ra + Rb)」より、(dv/dt)の式に変換して、Runge-Kutta法で速度と移動距離を求めていきます。

この速度のグラフは以下のとおり。

最終的に18.2〜18.3 m/sあたりで値が安定しています。これは何を表しているのか、というと時間がたてば落下する雨粒は等速になる、ってことなんです。ようは、空気抵抗があるから、雨が降っていたときに無尽蔵に速度アップする雨粒に脳天を貫かれなくてすむわけです。この等速値を「終端速度」と呼びます。

しかし、「Runge-Kutta法」って便利ですねぇ。これで、ようやく質量と空気抵抗の関係を検証することができるようになりました。

空気抵抗と戯れる（2006/06/01）

早くも今年も真ん中に突入しました。ただいま、物理の世界にどっぷりつかってます。

さて、問題。同じ大きさのボールがあるとして、１つが軟球、もう１つが鉄球とします。つまり、質量だけが異なる物体です。これを同じ高さから同時に初速度0で落下させるとどちらが先に着地するでしょう？

よくある質問なんですが、答えは「空気抵抗がない場合（真空内など）は、同時」です。では空気抵抗があればどうでしょう？

空気抵抗は「慣性抵抗」と「粘性抵抗」の２種類があります。方向としては、進行方向とは常に逆方向です。

慣性抵抗

S  : 進行方向に接する断面積（m^2）
Cd : 空気抵抗係数。雨粒だと0.2くらい、車だと0.35-0.40くらい
p  : 媒体である空気の密度（kg/m^3）。空気中なら約 1.2
V  : 速度（m/s）

このとき、慣性抵抗の力は
Ra = S x Cd x P x (V^2) / 2

粘性抵抗

これは、形状の形により異なるみたいです。球の場合は以下のとおり。

μ : 粘性率。媒体が空気の場合、1.8 x 10^(-5)
r : 球の半径（m）

このとき、粘性抵抗の力は
Rb = 6 x π x r x μ x V

で、この「慣性抵抗 + 粘性抵抗」が進行方向にはむかう空気抵抗となります。ですが、この空気抵抗は質量が絡んでいませんよね。ってことは、質量のみが異なる球体は、空気抵抗の影響は関係なさそうです。ん？ということは、空気抵抗が存在した世界でも、質量のみが異なる同じ大きさの物体を同じ高さから落としたとすると、同時に着地するの？

ここで、重力加速度gで下方向に進む質量mの球体は、

m・(dv/dt) = (m x g) - (Ra + Rb)

の式が成り立ちます。ここで質量が出てきました。つまり、ここでようやく「空気抵抗がある場合は、質量が重い物体のほうが速く落ちる」ということが証明できそうです。左の加速度がキーポイントかな？ごめんなさい、まだ、このあたりは飲み込めてないです。

ということで、空気抵抗がある世界で、同じ大きさの質量だけ異なる球を落下させた場合は重いほうが先に地面に着地する、らしい。

う〜む、正直このへん理解してなかったのですが(^_^;;、行き着く先は「F=ma」なんですな。これだけのために結構調べまくりました。間違いがあったらすみません(^_^;;。ちなみに、速度が遅い場合は粘性抵抗がよく効き、速度が速くなると慣性抵抗が良く効くようになるとのこと。ミサイルとか車とかなら（速いんで）、慣性抵抗をいかに抑えるかが大事なんですね。速度を求めるためには、進行方向の表面積と空気抵抗係数をいかに抑えるか、がいろいろ研究されているみたいで。なんとなく、空気抵抗の仕組みが分かってきた、かな。ネットで調べていたら、このへんは高校レベル、なんてのが多くって結構へこみました（苦笑）。私のときは、高校ではやってなかったかなぁ、う〜む。