この記事は，Simutrans Advent Calendar 14日目の記事です．Simutransのソースコードを読んで最高速度を求めるという話をします．なおここで対象にするのはいわゆるstandard版と言われるものであり，OTRP版やExtendedについては触れません（もっともOTRP版についてはほとんど差異がありません）．対象としたコードはgithubからpullしてきたもので，リビジョン番号は10270です．
編成の加速やらを取り扱っているメソッドはsimconvoi.ccにあるconvoiクラスのcalc_accelerationである．ここで加速度に関係する値としてresidual_powerがあり，これを計算するのにres_power関数を用いている．
calc_accelerationメソッドで呼び出されている）res_power関数を見る
res_power関数は次のようになっている．残余出力\( W \)，現在の速度\( s \)，全出力\( p \)，総摩擦重量\( m \)および総重量\( M \)として，
\begin{equation}
W=p - \left\lfloor
\frac{\lfloor s\rfloor\times\left(ms/3125 + 1\right)}{2048}
+\frac{64M}{1000}
\right\rfloor
\end{equation}
ただし
速度\(s\)の単位は内部速度であり，これと時速との関係はsimutil.hで記述されている．
calc_accelerationメソッドを見る
……とはいってもちゃんと見るのは速度を計算する部分だけである．
これは761行のif(sum_gesamtweight != 0) {から始まる部分である．
residual_powerを求める
これはsint64 residual_power = res_power(...)と，計算した内容を代入しているだけなので省略．
delta_vを求める
これが速度の差分となる（正確にはこれの\(2^{-12}\)倍）．
delta_vを\(\Delta_v\)，前の計算したタイミングから経過した内部時間（だと思う）delta_tを\(\Delta_t\)として
\begin{equation}
\Delta_v=W\times\lfloor\Delta_t\rfloor\times1000/\lfloor M\rfloor
\end{equation}
これの後delta_vにprevious_delta_vを加える．
次の部分を見れば分かるが，このprevious_delta_vは前のタイミングで直接足しあわされなかった（言うならば小数点以下の）部分である．
previous_delta_vにdelta_vの下12bitを保存する
previous_delta_v = (uint16) (delta_v & 0x00000FFFll);として保存している．
akt_speedを求める
akt_speedは現在速度である．
akt_speedを\(v\)としてakt_speed_sollの\(2^{-4}\)または\(v+\lfloor\Delta_v/2^{12}\rfloor\)の大きい方に設定する．
上の内容を見ると，res_powerが\(=0\)となるときに加速度が0になり，速度が極値を取ると分かる．速度の極値は普通は車両が「理想的に到達可能な最高速度」であるから，res_powerが\(=0\)となるような速度\(v_\text{max}\)を求めると，これが理想的に到達可能な最高速度である筈である．
なお上で「理想的に到達可能な最高速度」と言ったのは，「車両に設定された最高速度」や「道路・線路・運河……などの道に設定された最高速度」などを無視した（あるいは十分に高いと見なす）ときの最高速度である．
res_power関数で\(W=0\)を代入して\(s\)について解くと「理想的に到達可能な最高速度」が分かる！ということである．というわけで解くと（ただしfloorについては無視して計算した），
\begin{equation}
s = \frac{-1+\sqrt{1+4 m \cdot \frac{2048}{3125}(p - 64M / 1000) }}{2m/3125}
\end{equation}
となる（ここで\(s\geq0\)を用いた）．摩擦が1である（すなわち直線上を走行している）のであれば\(M=m\)である．
また\(p\)について解くと「理想的に到達可能な最高速度」に到達するために必要な全出力が分かる．これは単に移項するだけで求めることが出来て \begin{equation} p = \frac{ s\times\left(ms/3125 + 1\right)}{2048} +\frac{64M}{1000} \end{equation} となる．
これを計算するためのexcelシートがあると便利ですよね？ というわけで作りました．
現実世界との対応関係も見ておこう． 時刻\(t\)での質点の位置\(x\)として，速度\(v\)が \begin{equation} v(t) := \dot{x}(t) \end{equation} と定義される．ここで文字の上にドットを付けたものは時間微分を表すものとする；すなわち \begin{equation} \dot{x}(t)= \frac{\mathrm{d}x(t)}{\mathrm{d}t} \end{equation} である．
質量\(m\)の質点について，速度\(v\)および質点に加えられた力\(F\)として \begin{equation} m\dot{v}(t) = F(t) \end{equation} と書くことが出来る．これを運動方程式という．
運動方程式の両辺を\(m\)で除算し，時刻\(t_0\)から\(t_1\)まで定積分することで \begin{equation} v(t_1)-v(t_0) = \int_{t_0}^{t_1} \mathrm{d}t\ \frac{F(t)}{m} \end{equation} と表せる． この式はすなわち時刻\(t_0\)から\(t_1\)での速度の増分（左辺）が右辺のように表せると言うことを意味する．
一方で運動方程式の両辺に\(v\)を掛けて時刻\(t_0\)から\(t_1\)まで定積分することで \[ \frac{1}{2}m\left( \left(v(t_1)\right)^2 - \left(v(t_0)\right)^2\right) = \int_{t_0}^{t_1} \mathrm{d}t \ F(t) v(t) \] という式を得る．ここで \begin{equation} T(t) := \frac{1}{2}m\left(v(t) \right)^2 ,\quad W(t,t_0) := \int_{t_0}^{t} \mathrm{d}\tau \ F(\tau) v(\tau) \end{equation} とすると上式は \begin{equation} T(t_1) - T(t_0) = W(t_1,t_0) \end{equation} と表される．このように定義された\(T(t)\)は時刻\(t\)における質点の運動エネルギーといい，\(W(t,t_0)\)は時刻\(t_0\)から\(t\)までに力が質点にした仕事という． また各時刻での単位時間あたりの仕事を仕事率といい，仕事率\(P(t)\)は \begin{equation} P(t):=\dot{W}(t,t_0) = F(t)v(t) \end{equation} と定義される．
さて，速度\(v\)がゼロではないときには \begin{equation} F(t) = \frac{P(t)}{v(t)} \end{equation} として時刻\(t\)での質点に加える力を速度と仕事率から導出することが出来て， 運動方程式から \begin{align} v(t_1)-v(t_0) &= \int_{t_0}^{t_1} \mathrm{d}t\ \frac{F(t)}{m} \notag\cr &= \frac{1}{m}\int_{t_0}^{t_1} \mathrm{d}t\ \frac{P(t)}{v(t)} \end{align} というように表すことが出来る． ここで仕事率が時刻によらず一定であり，かつ\(v_0 = v(t_0) \simeq v(t_1) \)と見なせる（速度の変化量が十分に小さい）と仮定すると \begin{equation} v(t_1) \simeq v_0 + \frac{P}{mv_0}(t_1-t_0) \label{速度の時間発展} \end{equation} と表すことが出来る．
(\ref{速度の時間発展})式とSimutransでの速度の時間発展を比較してみよう． Simutransでの速度の変化量は\(\Delta_v\)で表したものであり， \[ \Delta_v=W\times\lfloor\Delta_t\rfloor\times1000/\lfloor M\rfloor \] であった． このうち\(M\)は質点の質量\(m\)と\(\Delta_t\)は\(t_1-t_0\)と次元が同一である． 一方で現実の運動方程式に含まれていた\(1/v_0\)という部分がSimutransでの速度の変化量の式に含まれていないことに注目せよ．実際の所，SimutransでkWのように仕事率で表されている単位はむしろ力のように見なす方が自然で，このとき\(W\)は力の次元を持つ\(P/v_0\)と対応される．
この立場に立つと，「ギア比を大きな値とすると出力が増加する」という一見おかしな振る舞いも「ギア比を大きな値とすると（出力ではなく）力が増加する」と，まっとうな振る舞いのように認識することが出来る．