クロソイド曲線で遊ぶ

# 背景
最近、「アートで魅せる数学の世界」（岡本 健太郎）という本を読んでいます。色々と楽しい図形を描いて遊んでいます。その中でクロソイド曲線が出てきました。概要は知っていたのですが、自分で書いたことがないな、と思ったので、少し勉強してみることにしました。

# クロソイド曲線とは
自動車を運転していて、一定の角速度でハンドルを回したときに自動車の描く軌跡と言われます。数学的には、曲率が一定の割合で変化する曲線のことです。

# クロソイド曲線を導出する
数学的には、曲率が一定の割合で変化する曲線ですので、自動車を持ち出す必要はありません。が、ハンドル角と曲率の関係も気になるので、自動車を持ち出すことにします。
「新版交通工学」（福田正編）によれば、

> ハンドルを$\phi'$回転すると、車の前輪は$\phi=k\phi'$($k$は定数)だけ向きを変える。

そうです。そして、前輪と後輪の間隔を$b$とすれば、後輪の中心が描く軌跡の曲線半径$R$は
$$
R = b/\tan\phi
$$
で与えられるそうです。
ここで、走行距離$L$と、接線の角度$\tau$を導入します。接線の角度が$\Delta \tau$だけ変化すると、走行距離の微小変化は
$$
\Delta L = R\Delta \tau
$$
となります。現在の接線角度を$\tau$とすれば、$x, y$方向の微小変化量は
```math
\begin{align}
\Delta x &= \Delta L \cos\tau \\
\Delta y &= \Delta L \sin\tau
\end{align}
```
となります。速さは一定であるとすれば$L=vt$となります。
これで必要な関係式は出そろいました。

それでは、まず$\tau$を求めてみます。積分してやると、
```math
\begin{align}
\tau(t) &= \int \frac{1}{R}dL\\
&= \int \frac{\tan\phi(t)}{b}dL \\
&= v\int \frac{\tan\phi(t)}{b}dt
\end{align}
```
となります。これを用いて、$(x,y)$を求めることが出来ます。
```math
\begin{aligned}
x(t) &= v\int \cos\left(v\int \frac{\tan\phi(t')}{b}dt'\right) dt \\
y(t) &= v\int \sin\left(v\int \frac{\tan\phi(t')}{b}dt'\right) dt
\end{aligned}
```
あとは、ハンドル角の時間依存性を決定してやれば軌跡を描くことが出来ます。
クロソイド曲線では、ハンドルを等角速度$\alpha$で回転します。すると、
$$
\phi(t) = k\alpha t
$$
と書けます。さらに、ハンドル角は小さいものとして、
$$
\tan \phi(t) \simeq \phi(t) = k\alpha t
$$
と近似します。すると、
$$
\tau(t) = \frac{vk\alpha}{b}\frac{t^2}{2}
$$
となります。これから、
```math
\begin{aligned}
x(t) &= v\int \cos\left( \frac{vk\alpha}{b}\frac{t^2}{2} \right) dt \\
y(t) &= v\int \sin\left( \frac{vk\alpha}{b}\frac{t^2}{2} \right) dt \\
\end{aligned}
```
となります。これがクロソイド曲線です。

# クロソイド曲線を描画する
それでは、クロソイド曲線を描画してみます。パラメータは全て$1$とします。
描画するためには、微分形式の方が扱いやすいので、
```math
\begin{aligned}
\Delta x &= \cos \tau(t) \Delta t \\
\Delta y &= \sin \tau(t) \Delta t \\
\tau(t) &= \frac{t^2}{2}
\end{aligned}
```
を用います。オイラー法を使って描画してみます。
```c++
#define _USE_MATH_DEFINES
#include <iostream>
#include <iomanip>
#include <cmath>

int main()
{
	int N = 1000000;
	double dt = 0.00001;
	double x = 0, y = 0, tau = 0;
	
	for (int i = 0; i < N; i++)
	{
		double t = i * dt;
		if (i % 100 == 0) std::cout << std::setprecision(15) << x << " " << y << std::endl;
		x += cos(tau) * dt;
		y += sin(tau) * dt;
		tau = t * t / 2.;
	}

	return 0;
}
```
![20165ade9510c71c8fe6b178eb68b425.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/520857/3e99e31a-4775-175b-34de-e352094da8f3.png)
出来ました。ハンドル角が小さいという近似の範囲を超えて描画していることには注意が必要です。

# ハンドル角が小さいという近似を撤廃してみる
一般の場合には次の式
```math
\begin{aligned}
\Delta x &= \cos \tau(t) \Delta t \\
\Delta y &= \sin \tau(t) \Delta t \\
\tau(t) &= \int \tan\phi(t)dt
\end{aligned}
```
を用います。クロソイドの時と同様に、$\phi(t)=t$とします。つまり、ハンドル角が小さいという近似のみ撤廃します。すると、
$$
\tau(t) = -\ln \cos t
$$
となります。明らかに$t=\pi/2$で発散してしまいますので、$t<\pi/2$がモデルの適用範囲です。
```c++
#define _USE_MATH_DEFINES
#include <iostream>
#include <iomanip>
#include <cmath>

int main()
{
	int N = 1000000;
	double dt = 0.00001;
	double x = 0, y = 0, tau = 0;
	
	for (int i = 0; i < N; i++)
	{
		double t = i * dt;
		if (t < M_PI / 2.)
		{
			if (i % 100 == 0) std::cout << std::setprecision(15) << x << " " << y << std::endl;
			x += cos(tau) * dt;
			y += sin(tau) * dt;
			tau = -log(cos(t));
		}
	}

	return 0;
}
```
![2e39a20362f27a5818b71123e5f8018d.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/520857/09d3ac63-6975-f3d3-f448-e200c7809218.png)

ハンドル角が小さいときはクロソイド曲線に近い軌跡となっています。$t<\pi/2$に限定しているので、クロソイドのように渦を巻いたりしません。