AutoLISPでのクロソイド曲線の計算

## クロソイド曲線とは
角加速度、速度が一定の曲線です。
道路設計の平面設計において、急なハンドル操作とならないように
小さい曲線半径の手前に設置する緩和曲線として知られます。
詳細な解説は『道路構造令の解説と運用』や他のWebページをご覧ください。


## 道路設計におけるクロソイド曲線の決定
道路設計の流れとしては、
**曲線半径の決定 → (緩和曲線長の設定) → パラメータの決定 → クロソイド作図**
が一般的かと思います。
今回は、片側を直線とした一般的なクロソイドを想定しています。

## CAD上での描画
CADではクロソイド曲線を表現するために任意の精度(分割数、端点)で
クロソイド曲線上の座標をプロットし、それを直線で結んだポリラインとして作成します。
分割数を増やせば精度は上がりますが、その分データ量が増えます。
100ぐらいで分割されていれば道路設計の範囲であれば事足りると思います。

## AutoLISPでの処理
プログラム全体の流れです。
1.　曲線半径、パラメータ、分割数、精度の入力
2.　曲線長の計算　※1から求まる。
3.　各座標の計算（再帰関数、サブルーチン）
4.　クロソイド曲線の座標リストの作成
5.　作図処理　※今回は省略

#### 各座標の計算
まず、クロソイド曲線上における座標の計算です。
流れとしては曲線半径、パラメータなどの入力ですが
先に再帰関数を用いた各座標値の計算のサブルーチンを作成します。

こちらの「[あの良寛 クロソイド曲線 VBScript](https://www.minc.ne.jp/~ryokan/vbscript/vb14.html)」様を参考にして作成しています。

`helix`：角度(τ)
`unit`：精度、有効数字
を事前に与えておくことが必要です。

後述のコードを参照ですが、
`helix`はクロソイドの条件と分割数から自ずと算出されます。

~~~lisp
(defun Pxloop (k)
  (if (= k 1)
    1.0
    (/ (* (Pxloop (1- k)) -1 (round (expt helix 2) unit)) (* (- (* 2 k) 2) (- (* 2 k) 3)))
  )
)

(defun Qxloop (k)
  (/ (Pxloop k) (- (* 4 k) 3))
)

(defun Pyloop (k)
  (if (= k 1)
    1.0
    (/ (* (Pyloop (1- k)) -1 (round (expt helix 2) unit)) (* (- (* 2 k) 1) (- (* 2 k) 2)))
  )
)

(defun Qyloop (k)
  (/ (Pyloop k) (* (- (* 4 k) 1)))
)
~~~

`round`関数を使用していますがAutoLISPでは標準の関数にないので、自作が必要です。
~~~lisp
(defun round (num unit / int)
  (setq int (expt 10 unit))
  (setq num (* num int))
  (setq num (+ num 0.5))
  (setq num (float (fix num)))
  (setq num (/ num int))
)
~~~

#### クロソイド関数の作成
`klothoide`関数を作成します。

入力としては
`rad`：曲線半径
`param`：パラメータ
`split`：分割数
`unit`：精度
が必要です。

重複しますが流れとしては以下の通りです。
1.　`total-len`（曲線長の計算）
2.　`total-len`を`split`で等分し、`unit-len`とする。
3.　等分位置での`helix`（角度）の計算
4.　`helix`から`Qxloop`、`Qyloop`による座標の計算
5.　繰り返してクロソイド座標をリストとして返す。

クロソイドの第N項と第N+1項の誤差が`margin`の範囲内に収まるように
再帰関数による計算をかけています。

~~~lisp
(defun klothoide
  (rad param split unit
   /
   margin unit-len len helix i pt ptlst
   k yk ptx pty xs ys xqk yqk
  )
  (setq total-len (round (/ (expt param 2) rad) unit))
  (setq unit-len (round (/ (float total-len) split) unit))
  (setq ptlst (list (list 0.0 0.0)))
  (setq len 0.0)
  (setq margin (round (expt 0.1 unit) unit))
  (repeat split
    (setq len (+ len unit-len))
    (setq rad (round (/ (expt param 2) len) unit))
    (setq helix (round (/ len (* 2 rad)) unit));τ=L/2R
    (setq  i 1
           xs 1.0
           xk (round (* param (sqrt (* 2.0 helix))) unit)
           ys (round (/ 1.0 3.0) unit)
           yk (round (* helix param (sqrt (* 2.0 helix))) unit)
           yqk (Qyloop 1)
           xqk (Qxloop 1)
    )
    (while (and (> (abs (* xk xqk)) margin)
                (> (abs (* yk yqk)) margin)
           )
      (setq i (1+ i))
      (setq xqk (Qxloop i)
            yqk (Qyloop i)
            xs (+ xs xqk)
            ys (+ ys yqk)
            pty (* yk ys)
            ptx (* xk xs)
      )
    );while
    (setq pt (list ptx pty))
    (setq ptlst (cons pt ptlst))
  );repeat
  (setq ptlst (reverse ptlst))
);klothoide
~~~

# 作図処理
作図処理は色々ありますが、[こちら](https://www.theswamp.org/index.php?topic=32148.0)が参考になるかと思います。
ただし、クロソイド曲線の座標計算が原点を元に計算しているため
現在のコードでは作図後にさらにCAD上で動かさなければいけません。

プログラムに組む際は、座標の移動や回転処理を追加する必要があります。
また直線から始まるクロソイドのみで曲線間同士をつなぐクロソイド曲線には対応していません。

# 雑記
初めてQiitaに投稿してみました。
AutoLISPの投稿が少ないので自分へのメモ書きとしての記載です。

実際に使うなら、起点-IP点-終点の三点指定して、
曲線半径とパラメータを入力して作図処理となりそう。

といいつつもCivil3DやAPS-MarkⅣ、V-roadの線形ソフトが既にあるから
今後も要らないんだろうなあと思いつつ後学のために。

## 参考文献

[クロソイド：あの良寛 クロソイド曲線 VBScript](https://www.minc.ne.jp/~ryokan/vbscript/vb14.html)
[クロソイド：DRA-CADアソシエイト クロソイドの設計と作図](https://www.kozo.co.jp/dracad/member/doc.php?p2=83)
[クロソイド：千三つさんが教える土木工学 緩和曲線](https://sites.google.com/view/senmitsu/%E6%B8%AC%E9%87%8F%E7%B3%BB/10-3-%E7%B7%A9%E5%92%8C%E6%9B%B2%E7%B7%9A)
[再帰関数：Manual Chair Japan](http://manualchair.html.xdomain.jp/index.html)
[作図処理：the swamp Entmake Functions (Lee Mac)](https://www.theswamp.org/index.php?topic=32148.0)