function [t0list,vmaxlist,blist,v0list,x0list] = makelists(tTlist, ...
    Tlist, mlist, CD, Ap, rho, t0, v0, x0, CDp, App)
%makelists Create the lists of transitions and initial conditions
%   The functions used on each part of the flight change with each
%   transition. The regions are, Pre-launch, initial thrust, final thrust,
%   coast to apogee, parachute descent from apogee. CDp and App are for the
%   parachute. For free fall (ballistic descent) make them the same as for
%   the thrust and coasting phase.
% For the numbering of parameters:
% 1 = pre-launch
% 2 = thrust 1
% 3 = thrust 2
% 4 = coast to apogee
% 5 = post-apogee

gloc = 9.80665; % local acceleration of gravity in m/s^2

t0list(1) = t0 - 100; % Start of pre-launch, not actually used in calculations
t0list(2) = t0; % End of pre-launch, start of thrust 1
t0list(3) = t0list(2) + tTlist(1); % End of thrust 1, start of thrust 2
t0list(4) = t0list(3) + tTlist(2); % End of thrust 2, start of coast

vmaxlist(1) = 0; % pre-launch
vmaxlist(2) = sqrt(2*(Tlist(1)- mlist(2)*gloc)/(rho*CD*Ap));
vmaxlist(3) = sqrt(2*(Tlist(2)- mlist(3)*gloc)/(rho*CD*Ap));
vmaxlist(4) = sqrt((2*mlist(4)*gloc)/(rho*CD*Ap)); % Terminal Velocity
vmaxlist(5) = sqrt((2*mlist(5)*gloc)/(rho*CDp*App)); % Terminal Velocity Parachute

blist(1) = -gloc;
blist(2) = (Tlist(1)/mlist(2) - gloc);
blist(3) = (Tlist(2)/mlist(3) - gloc);
blist(4) = gloc;
blist(5) = -gloc;

v0list(1) = v0;
v0list(2) = v0;
v0list(3) = velocity(v0list(2), vmaxlist(2), blist(2), 0, tTlist(1));
v0list(4) = velocity(v0list(3), vmaxlist(3), blist(3), 0, tTlist(2));
v0list(5) = 0; % vapogee = 0 by definition

tapogee = vmaxlist(4)/gloc*atan(v0list(4)/vmaxlist(4)); % length of stage 4
t0list(5) = t0list(4) + tapogee; % Start of descent from apogee

x0list(1) = x0;
x0list(2) = x0;
x0list(3) = position(x0list(2), v0list(2), vmaxlist(2), blist(2), 0, tTlist(1));
x0list(4) = position(x0list(3), v0list(3), vmaxlist(3), blist(3), 0, tTlist(2));
x0list(5) = vmaxlist(4)^2/gloc*log(sqrt(1+v0list(4)^2/vmaxlist(4)^2))+x0list(4);

arg1 = gloc*x0list(5)/vmaxlist(5);
if arg1 > 100 % if exp(arg1) => Inf
    arg2 = log(2) + arg1; % Use simplified large number version
else
    arg2 = acosh(exp(arg1)); % Use full expression
end
timpact = vmaxlist(5)/gloc*arg2;
t0list(6) = t0list(5) + timpact; % full time to impact
v0list(6) = velocity(v0list(5), vmaxlist(5), blist(5), 0, timpact); %

end

