1998/6/15作成
1999/1/3に追加
1999/9/16に追加
NekoFlightのモデル計算部は、C++で組まれています。このプログラムの元になったUNIX版のCプログラムはNekoFlightの歴史のページで公開していますので、詳細を知りたい方はそちらを解析してください。人に見せるとは思ってなかったのでとんでもなくひどいソースです(^^;
この手のシミュレータ系のプログラムで基本となるのは、運動を表す常微分方程式を解くことです。といってもそんなに難しいことではありません。
まず物体の位置をp、速度をv、加速度をaとすれば、
dp/dt = v, dv/dt = aという常微分方程式が基本となります。
すなわち、位置の変化率が速度で、速度の変化率が加速度です。
この常微分方程式をどうやって解くかですが、ゲームなどでは計算精度がそれほど重要でないのでオイラー法で十分です。オイラー法なんていう凄そうな名前がついてますが話は簡単で、
{
double dt = 0.1;
double x = 0; y = 0; vx = 1; vy = 0; ax = 0; ay = 0.98;
for (;;) {
x += vx * dt;
y += vy * dt;
vx += ax * dt;
vy += ay * dt;
}
}
これだけです(^^;。単に加速度に時間幅を掛けて速度に足して、同様に位置に速度を足すだけ。本当は精度を気にするといろいろ問題あるのですが、ゲームならこんなのでＯＫです。
この例だと、ｙ方向への加速度一定ですので放物運動となります。
なお、オイラー法というのは常微分方程式の１次解法です。すなわち、線形近似です。常微分方程式の解放には他にルンゲクッタなどありますので精度を気にする方はこちらを使いましょう。また、オイラー法は単独１階の常微分方程式の解放ですが、上で示したように高階でも連立に分解すればＯＫですし、連立の場合も、それぞれをオイラー法で解けばＯＫです。
このあたりの細かいことは、数値解法に関する教科書を買ってきましょう。
ちょとしたシミュレータなら、ここで書いた程度の知識でＯＫですが。
この手のフライトシミュレータで問題になるのが、機体の向きをどの様に表すかです。
グラフィックスライブラリなどでは、オブジェクトの向きとオブジェクトの上方の二つの単位ベクトルで表現することが多いです。NekoFlightでは、ミサイルの向きにこの表現法を使ってますが、フライトシミュレータでの向き表現には、実際はオイラー角を用いることが一般的です。
−−−−−−−−−−−−−−−−−−−
'99/1/3に追加
向 博幸 氏（nbe00760@nifty.ne.jp）から、物体回転処理に関して良いことを聞いたので補足。オイラー角だと、真上や真下で極が発生するとか変換式が複雑などの問題があります。ですので、回転処理に関しては向さんのＨＰにあるような行列を用いる方法がお勧めです。ということで、Ver0.9の時点でのNekoFlightでは機体の回転周りオイラー角で処理してますが、そのうち機体回転角を行列で表すようにする予定です。
一応参考になるかもしれないので、前に書いたオイラー角の話も残しておきます。NekoFlightでオイラー角を使っている経緯に関しては、NekoFlightの歴史のページを参照してください。
−−−−−−−−−−−−−−−−−−−
'99/9/16に追加
行列処理の具体的な話を向氏のＨＰの方で知ることができます。参照してください。
−−−−−−−−−−−−−−−−−−−
オイラー角では、３つのひねり角で機体の向きを表します。機体座標からワールド座標への変換を行うには、まずピッチ方向にひねり、ロール方向にひねり、地面垂直方向にひねります。
ですので、機体のロール角がa、ピッチ角がb、機首方位がcの場合、機体座標→ワールド座標の変換行列要素は
{
y00 = cos(b) * cos(c) - sin(a) * sin(c) * sin(b);
y01 = -cos(b) * sin(c) - sin(a) * cos(c) * sin(b);
y02 = -sin(b) * cos(a);
y10 = cos(a) * sin(c);
y11 = cos(a) * cos(c);
y12 = -sin(a);
y20 = sin(b) * cos(c) + sin(a) * sin(c) * cos(b);
y21 = -sin(b) * sin(c) + sin(a) * cos(c) * cos(b);
y22 = cos(b) * cos(a);
}
となり、実際に座標変換するには、
void change3d3(double x, double y, double z, double* x1, double*
y1, double* z1)
{
*x1 = x * y00 + y * y01 + z * y02;
*y1 = x * y10 + y * y11 + z * y12;
*z1 = x * y20 + y * y21 + z * y22;
}
となります。
オイラー角を使えば、上のaがロール角、bがピッチ角、cが機首方位となり、ＨＵＤの表示に使うにはこれでちょうど良いですが、機体の実際の回転計算にはこのままでは使えません。機体の座標系におけるロール、ピッチ、ヨー方向の回転変化をオイラー角変化に変換する必要があります。
具体的には、ロール角速度、ピッチ角速度、ヨー角速度をvaのx,y,zメンバ、ロール角変化率（操縦桿の左右）をam.x、ピッチ角変化率（操縦桿の前後）をam.y、ヨー角変化率（フットペダルの位置）をam.z、各軸の慣性モーメント量をix, iy, iz、機体のオイラー角をa.x, a.y, a.zとすると、次の様にして機体の向きが計算できます。
{
va.x += am.x / ix * dt;
va.y += am.y / iy * dt;
va.z += am.z / iz * dt;
double ax = a.x;
double ay = a.y;
a.x += (va.x * cos(ay) + va.z * sin(ay)) * dt;
a.y += (va.y + (va.x * sin(ay) - va.z * cos(ay)) * sin(ax) / cos(ax)) * dt;
a.z += (-va.x * sin(ay) + va.z * cos(ay)) * dt / cos(ax);
}
実際のNekoFlightでは、amは操縦桿の位置そのままではなく、角翼の計算結果から求まった角軸の回転速度変化率となります。amに操縦桿やペダルの位置を用いて、vaを適当に減衰させてもある程度それらしくなります（実はNekoFlightのVTOL時などはこれで処理してます）。
ここでいう減衰処理とは、たとえば、空気抵抗で速度がゆっくり減少するというような処理をどうやるのかということです。
シミュレータにおいては、この手の減衰処理にダンパーを仮定すると綺麗に表現できます。
ダンパーとは、速度に比例して反対方向に力がかかるようなもので、バイクのショックアブソーバーのダンパーと概念的には同じものです。
具体的には、速度をその速度自身に比例させて減少させると実際の減衰運動っぽくなります。たとえば、平面上をオブジェクトが動き回っており、慣性が働き、同時に抵抗も働いているという平面状のビー玉のような動きを簡単に表すには、
{
double dt = 0.1;
double vx = 1.23; vy = 4.2; x = 0; y = 0;
for (;;) {
x += vx * dt;
y += vy * dt;
// ここでキー操作などに伴いvxやvyを変化
vx -= vx * 0.5; // 速度減衰率
vy -= vy * 0.5;
}
}
などとすると、慣性を持って動いている物体らしい動きになります。
この処理はなかなか便利で、NekoFlightではミサイルシーカーの動きや機銃の向き、操縦桿の位置の処理などにこれを用いています。
この減衰処理と似たものとして、振動処理があります。ダンパーに加えて、位置に比例して力を受けるばねもモデルに含ませたものです。これを用いると、たとえば単振動や他のものを追尾するような動きが表現できます。
{
double dt = 0.1;
double vx = 1.23; vy = 4.2; x = 0; y = 0;
for (;;) {
x += vx * dt;
y += vy * dt;
vx -= x * 0.1; // ばね定数
vy -= y * 0.1;
vx -= vx * 0.5; // 減衰率
vy -= vy * 0.5;
}
}
これは減衰ありばねのシミュレーションになりますので、結局減衰振動のシミュレーションになります。パラメータによってはオーバーシュートなどが表現できます。
まじめにシミュレーションするには、力学を考慮し、各力点に加わる力の計算を行う必要があります。
このあたりは高校物理そのまま（結局、a=F/m）で、物体に働く力が求まったら上に書いた常微分方程式の数値解法で、数値積分し物体の速度と位置を求めます。
回転を無視する場合これで終わりですが、物体の回転も考える場合、剛体計算も考慮する必要があります（これに関しては大学の教養で習います）。
といってもこれも簡単で、重心（or支点）と作用点間の位置ベクトルと、力のベクトルの外積を求め、これを慣性モーメントで割ることで、各軸の回転速度変化量が求まります。あとは物体の位置計算同様、二回数値積分することで、物体の向きが求まります（上の機体の回転参照）。
このあたりは、力学の教科書を参照しましょう。
フライトシミュレータでは、翼の力学計算を行う必要があります。まじめに行うといろいろ面倒ですが、ゲームに使う程度ならば、揚力と抗力を迎角と翼速度から単純に求めることができます。
まず、翼に働く力は次の図の様に３つの力があります。ただし、エンジンがある場合、これに推力も加わります。
ここで、抗力は翼の空気抵抗から発生する力で向きは翼の進行方向の反対方向、揚力は抗力に垂直、翼上向へ働く力です。重力はいうまでもなく地面の方向に働きます。また、迎角というのは、翼の進行方向と翼の基本面との角度です。この図では表してませんが、実際は速度には奥行き方向への速度成分もあり、翼水平方向における進行方向と翼垂直面との角度というのも考えられます。これは横滑り角と呼ばれます。
ここで、各力の求め方は、
抗力＝抗力係数＊翼速度＾2＊翼面積＊空気密度／2
揚力＝揚力係数＊翼速度＾2＊翼面積＊空気密度／2
重力＝翼の質量＊重力加速度
です。なおここでの翼速度というのは、翼の水平方向成分を無視した進行方向への速度です（つまり横滑り成分に関しては無視します）。結局、抗力、揚力ともに、速度の二乗に比例します。
この式で問題になるのが抗力係数と揚力係数です。これはいろいろな要因で変化するのですが、細かいことを無視すれば、主に迎角で決まります。迎角とこれらの係数の関係は次のようになります。
揚力係数は大体、迎角に比例し、抗力係数は迎角の２乗に比例します。
ただしこれは失速角度内の話で、迎角が失速角（大体±0.4rad程度）を越えると、揚力係数が０になり、抗力係数が非常に大きくなります。
これだけ知っていれば単純翼の計算が可能ですが、詳細はNekoFlightの歴史で紹介した航空機力学入門などを参考にしましょう。
機体全体の計算を行うには、手抜き計算で良いのであれば、機体全体を一つの翼とみなして、上の翼計算に推力を追加するだけでＯＫです。ただしこの方法だと、機体の回転周りはまったく計算できませんので、操縦桿を右に倒したら単にロール速度を増すなどのいんちき計算でごまかします(^^;これで処理する場合、上で書いた減衰処理などが役に立ちます。
もう少しまじめに計算したい場合は、主翼の左右と水平、垂直尾翼をそれぞれ別の翼とし計算し、あとはそれらの剛体計算で処理します。これだと、操縦特性もある程度それらしくなります。
各軸の操縦に関しては、操縦桿やペダルにあわせて、各翼をひねります。あとは翼単位で上で述べた計算を行い、力学計算の部分で述べた剛体計算で処理します。
ここで注意することは、機体自体も回転しているため、各翼の速度は、機体重心の速度と差があるということです。これを無視しても平気ですが、幾つか操作副作用（ピッチ操作がヨー操作を伴うなど）の計算で変なことになります。
ミサイル周りやランドスケープ処理、あたり判定などの話もしたいのですが、とりあえず休憩(^^;
そのうち追加します。