Steering axis and Ackermann

%% Steering axis and Ackermann: the line a wheel steers about, and how far each wheel steers
%
% Companion to the documentation page "Steering axis and Ackermann". Runs unchanged in
% GNU Octave (7 or later) and MATLAB. Edit the car in example_suspension.m, or type your
% own hardpoints from the Kinematics tool, and run:
%
%   >> steering_axis
%
% What this computes:
%
% 1. The steering axis, through the two ball joints, at the static ride (pd_metrics):
%    its caster and inclination, where it meets the road against the contact patch
%    (scrub radius and mechanical trail), and against the wheel centre at hub height.
% 2. What steering does to the wheel. The lower wishbone is held still and the rack moved
%    2 mm either way, as a wheel turns with the car standing on it: camber per degree of
%    steer, and how far the patch drops, which lifts the car (jacking).
% 3. Ackermann. Each front wheel's steer at a rack position, the wheel centre held at the
%    static ride. The inside wheel of a turn should steer more: with both wheels rolling
%    about one point on the rear axle's line,
%        cot(outer) - cot(inner) = track / wheelbase,
%    and the app's percentage is the left side over the right. The app finds the static
%    ride at each rack position by a search that stops within 1e-4 mm of it; the script
%    solves where it stopped (example_suspension.m gives it), so the two agree to their
%    last digits. Solved at the ride itself, the 10 mm row moves in its sixth digit.
%
% Part of the PaceDynamics documentation, https://pacedyn.dev/docs
% SPDX-License-Identifier: MIT

[car, app] = example_suspension();
g = 9.81;
Fz = car.massTotal * g * [car.frontMassFrac, 1 - car.frontMassFrac] / 2;   % N, per wheel
R0 = car.tyreFreeRadius - Fz / car.tyreRate;                  % mm, the radius each axle stands on

axles = {pd_corner(car.front), pd_corner(car.rear)};
read = @(a, p, R) pd_metrics(p.at(a.ubj), p.at(a.lbj), p.at(a.C), p.n, p.at(a.tieOuter), R);
for ax = 1:2
  a = axles{ax};
  s(ax) = read(a, pd_upright(a, 'travel', 0), R0(ax));
end

% Steering with the lower wishbone held where it stands (its angle 0), rack 2 mm each way.
f = axles{1};
lo = read(f, pd_upright(f, 'beta', 0, -2), R0(1));
hi = read(f, pd_upright(f, 'beta', 0,  2), R0(1));
camberPerSteer = (hi.camber - lo.camber) / (hi.toe - lo.toe);
jacking = -(hi.patch(3) - lo.patch(3)) / (hi.toe - lo.toe);   % mm per degree, + = the car rises

% Ackermann. The rack moves the left tie rod inboard as it moves the right one outboard;
% each corner's own frame has y outboard, so the left wheel sees -rack and the right +rack.
% A wheel steering left turns its spin axis toe-out on the left and toe-in on the right.
track = 2 * abs(car.front.wheelCentre(2));
ideal = track / car.wheelbase;
cotd_ = @(x) 1 ./ tand(x);
for i = 1:numel(app.steering.ackermann.rack)
  r = app.steering.ackermann.rack(i);
  pl = pd_upright(f, 'travel', app.steering.ackermann.travelInner(i), -r);
  pr = pd_upright(f, 'travel', app.steering.ackermann.travelOuter(i),  r);
  ml = read(f, pl, R0(1));  mr = read(f, pr, R0(1));
  left  = -(ml.toe - s(1).toe);                   % deg into a left turn
  right =   mr.toe - s(1).toe;
  inner(i) = left;  outer(i) = right;             % turning left: the left wheel is inside
  pct(i) = 100 * (cotd_(outer(i)) - cotd_(inner(i))) / ideal;
end

%% Against the app
ax = {'front', 'rear'};
for k = 1:2
  A = app.steering.(ax{k});
  pd_compare(['Caster, ' ax{k} ' [deg]'],             s(k).caster,        A.caster);
  pd_compare(['KPI, ' ax{k} ' [deg]'],                s(k).kpi,           A.kpi);
  pd_compare(['Scrub radius, ' ax{k} ' [mm]'],        s(k).scrub,         A.scrubRadius);
  pd_compare(['Mechanical trail, ' ax{k} ' [mm]'],    s(k).trail,         A.mechanicalTrail);
  pd_compare(['Kingpin offset, ' ax{k} ' [mm]'],      s(k).kingpinOffset, A.kingpinOffset);
  pd_compare(['Caster offset, ' ax{k} ' [mm]'],       s(k).casterOffset,  A.casterOffset);
end
pd_compare('Steering arm, front [mm]',                s(1).steeringArm,   app.steering.front.steeringArm);
pd_compare('Camber per steer [deg/deg]',              camberPerSteer,     app.steering.front.camberPerSteer);
pd_compare('Steer jacking [mm/deg]',                  jacking,            app.steering.front.steerJacking);
for i = 1:numel(app.steering.ackermann.rack)
  r = app.steering.ackermann.rack(i);
  pd_compare(sprintf('Inside wheel at %g mm rack [deg]', r),  inner(i), app.steering.ackermann.inner(i));
  pd_compare(sprintf('Outside wheel at %g mm rack [deg]', r), outer(i), app.steering.ackermann.outer(i));
  pd_compare(sprintf('Ackermann at %g mm rack [%%]', r),      pct(i),   app.steering.ackermann.percent(i));
end
pd_compare();