M = 194295; % 列车质量 kg
highSpeed = 50/3.6; % 最初加速到的最高速度
L = 10000;
S = linspace(s0,s1,L);
Ssize = size(S);
V = zeros(Ssize);
E = zeros(Ssize);
F = zeros(Ssize);
% 终点制动曲线
[ endBrakingCurveV,endBrakingCurveS ] = brakingCurveFun( s0, s1, 0, gradient,
curvature );
E(1) = 0; % 消耗能量初始化
V(1) = 0; % 初始速度为0
% 牵引阶段
i = 2;
while (i<length(V) && V(i-1) < highSpeed)
[ Fmax ] = maxTractionFun( V(i - 1) );
[W] = totalResistanceFun(V(i - 1), S(i-1), gradient, curvature);
capacityMaxA = (Fmax - W) / M; % 能够达到的最大加速度
if capacityMaxA > 1
a = 1; % 实际加速度, 因为题目限制最大加速度不能超过1
totalF = M * a; % 合力
F(i - 1) = totalF + W; % 实际牵引力
else
a = capacityMaxA;
F(i - 1) = Fmax;
end
V(i) = sqrt((V(i-1))^2 + 2 * a * (S(i-1) - S(i)));
if checkMeetEndBrakingCurve(S(i),V(i),endBrakingCurveS,endBrakingCurveV)
V(i-1:end) = interp1(endBrakingCurveS,endBrakingCurveV,S(i-1:end),'pchip');
brakingTerminal = [brakingTerminal;S(i-1),S(end)]; % 记录制动区间
E(i:end) = E(i-1);
i = length(V);
elseif
checkMeetBrakingCurve(S(i),V(i),brakingCurveS,brakingCurveV,curveTerminal)
index = find(curveTerminal(:,1)>S(i),1);
tempEnd = curveTerminal(index,2);