function j = hcstep(s, v, a, y, e, m, dt)
%% One step of harmonic predictive control
% s,v,a,j are stock, velocity, acceleration, jerk.
% y < e are target yield and extinction; assume e < s < y.
% m > 0 is maximum allowed acceleration
% dt is time granularity
if v < 0; F = e; else F = y; end                % Determine focus
if abs(v*dt) > abs(F - s); j = 0; return; end	% Proximity overflow
p = v/(F - s);                  % Proximity p satisfies 0 <= p <= 1/dt
h = [m - a,  -(3*a + 2*v*p)];   % Headroom h (both cases)
j = p*h(1 + (F == y));                        % Jerk = Proximity times Headroom
end

 
