Claude Codeでクロソイドの線形計算を書き、道路構造令の値で検算する

# 線形計算のルール
- 座標は x=東, y=北、角度の内部表現はラジアン、長さはメートル。
- 測量座標(X=北, Y=東)との変換は入出力の境界だけで行う。
- 左折・右折を符号 sgn(+1/-1)で一本化し、右折は左折の鏡像として扱う。
- 変更のたびに `python3 -m unittest` を実行し、全件通るまで完了としない。
- 級数と数値積分の差が 1e-6 m を超えたら、式ではなく係数の写し間違いを疑う。
- 道路構造令の値(第15条・第18条)は DESIGN_LIMITS に集約し、コードの途中に数字を直書きしない。

alignment.py を作って。入力は IP0→IP1→IP2 の座標、半径 R、パラメータ A。
標準ライブラリだけで、TS・SC・CS・ST の座標、TL、Lc、Ls と、
TS からの追加距離 s の中心線座標を返す関数を書いて。
test_alignment.py に、級数と Simpson 積分の一致、ST が IP から TL の位置にあること、
SC/CS が半径 R の円上にあること、右折が左折の鏡像になることのテストを書いて。
まずテストだけ書いて、実装はその後に。

import math
from dataclasses import dataclass
 
# 道路構造令 第15条(曲線半径)・第18条(緩和区間長)の表
# 設計速度 km/h -> (曲線半径 m 以上, 緩和区間長 m 以上)
DESIGN_LIMITS = {120: (710, 100), 100: (460, 85), 80: (280, 70), 60: (150, 50),
                 50: (100, 40), 40: (60, 35), 30: (30, 25), 20: (15, 20)}
 
 
def spiral_local(length, a):
    """TS を原点、接線方向を x 軸とする局所座標(級数 4 項)。"""
    l, a2 = length, a * a
    x = l - l**5/(40*a2**2) + l**9/(3456*a2**4) - l**13/(599040*a2**6)
    y = l**3/(6*a2) - l**7/(336*a2**3) + l**11/(42240*a2**5) - l**15/(9676800*a2**7)
    return x, y
 
 
def _place(origin, heading, lx, ly):
    c, s = math.cos(heading), math.sin(heading)
    return (origin[0] + lx*c - ly*s, origin[1] + lx*s + ly*c)
 
 
@dataclass
class Alignment:
    ts: tuple; sc: tuple; cs: tuple; st: tuple
    tl: float; ls: float; lc: float; delta: float
    d1: float; d2: float; sgn: int; radius: float; a: float
    center: tuple
 
    @property
    def total_length(self):
        return 2*self.ls + self.lc
 
    def point(self, s):
        """TS からの追加距離 s の中心線座標。"""
        ls, lc, r = self.ls, self.lc, self.radius
        if s <= ls:
            lx, ly = spiral_local(s, self.a)
            return _place(self.ts, self.d1, lx, self.sgn*ly)
        if s <= ls + lc:
            ang0 = math.atan2(self.sc[1]-self.center[1], self.sc[0]-self.center[0])
            ang = ang0 + self.sgn*(s - ls)/r
            return (self.center[0] + r*math.cos(ang), self.center[1] + r*math.sin(ang))
        lx, ly = spiral_local(self.total_length - s, self.a)
        # ST から逆向きに歩くと旋回方向が反転する
        return _place(self.st, self.d2 + math.pi, lx, -self.sgn*ly)
 
 
def solve(ip0, ip1, ip2, radius, a):
    ls = a*a/radius
    d1 = math.atan2(ip1[1]-ip0[1], ip1[0]-ip0[0])
    d2 = math.atan2(ip2[1]-ip1[1], ip2[0]-ip1[0])
    delta = (d2 - d1 + math.pi) % (2*math.pi) - math.pi
    sgn = 1 if delta > 0 else -1
    i = abs(delta)
    tau = ls/(2*radius)
    xs, ys = spiral_local(ls, a)
    dr = ys - radius*(1 - math.cos(tau))
    xm = xs - radius*math.sin(tau)
    tl = (radius + dr)*math.tan(i/2) + xm
    lc = radius*(i - 2*tau)
    if lc < 0:
        raise ValueError("交角に対して R と A が大きすぎる(円曲線長が負)")
    ts = (ip1[0] - tl*math.cos(d1), ip1[1] - tl*math.sin(d1))
    st = (ip1[0] + tl*math.cos(d2), ip1[1] + tl*math.sin(d2))
    sc = _place(ts, d1, xs, sgn*ys)
    cs = _place(st, d2 + math.pi, xs, -sgn*ys)
    h = d1 + sgn*tau
    center = (sc[0] - sgn*radius*math.sin(h), sc[1] + sgn*radius*math.cos(h))
    return Alignment(ts, sc, cs, st, tl, ls, lc, delta, d1, d2, sgn, radius, a, center)
 
 
def check_design_speed(radius, ls, speed):
    rmin, lmin = DESIGN_LIMITS[speed]
    return radius >= rmin and ls >= lmin

AssertionError: 832.8297414353683 != 1042.043242878848 within 6 places (209.2135014434798 difference)

test_arc_endpoints_continuous ... ok
test_cs_matches_circle_walk ... ok
test_design_speed_60 ... ok
test_right_turn_is_mirror ... ok
test_sc_cs_on_circle ... ok
test_series_matches_integration ... ok
test_tangent_points_lie_on_tangent_lines ... ok
 
Ran 7 tests in 0.003s
OK

EXPECTED = {  # 線形計算書の主要点(単位 m)。値は自分の計算書のものに置き換える
    "TS": (840.929, 0.000),
    "SC": (940.652, 5.545),
}
 
def test_against_report(self):
    for name, (x, y) in EXPECTED.items():
        px, py = getattr(self.c, name.lower())
        self.assertAlmostEqual(px, x, places=3)
        self.assertAlmostEqual(py, y, places=3)