function [a,v,x] = posvelacc(t,t0list,vmaxlist,blist,v0list,x0list)
%posvel Calculate the position and velocity for const thrust models
%   This function was made from the TwoStepThrustAnalytical.m file
%   t is the array of times at which to calculate a, v, & x
%     For the numbering of parameters:
%     1 = pre-launch
%     2 = thrust 1
%     3 = thrust 2
%     4 = coast to apogee
%     5 = post-apogee
%   The t0list contains the transition times from phase to phase
%   The vmaxlist contains the vmax or vt for each phase
%   The blist contains the b or g for each phase
%   The v0list is the initial velocity at the start of each phase
%   The x0list is the initial altitude for each phase

% Split time array into pre-launch, Thrust 1, Thrust 2, Coast, Freefall
tpre = t(t<t0list(2));
t1 = t(t>=t0list(2) & t<t0list(3));
t2 = t(t>=t0list(3) & t<t0list(4));
tcoast = t(t>=t0list(4) & t<t0list(5));
tff = t(t>=t0list(5));

% Set the pre-launch conditions
apre = 0*tpre; % apre is assumed to be zero.
vpre = apre + v0; % v0 is most likely 0. If not, you need to modify the
                  % code for xpre
xpre = apre + x0; % Strictly speaking xpre can only be constant if
                  % vpre is zero

% Calculate the velocities first. They're needed for the accelerations.
v1 = velocity(v0list(2), vmaxlist(2), blist(2), t0list(2), t1);
v2 = velocity(v0list(3), vmaxlist(3), blist(3), t0list(3), t2);
vcoast = velcoast(v0list(4), vmaxlist(4), blist(4), t0list(4), tcoast);
vff = velocity(v0list(5), vmaxlist(5), blist(5), t0list(5), tff);
v = [vpre v1 v2 vcoast vff];

% Calculate the positions (altitudes) second
x1 = position(x0list(2), v0list(2), vmaxlist(2), blist(2), t0list(2), t1);
x2 = position(x0list(3), v0list(3), vmaxlist(3), blist(3), t0list(3), t2);
xcoast = poscoast(x0list(4), v0list(4), vmaxlist(4), blist(4), t0list(4), tcoast);
xff = position(x0list(5), v0list(5), vmaxlist(5), blist(5), t0list(5), tff);
x = [xpre x1 x2 xcoast xff];

% Calculate the accelerations
a1 = blist(2)*(1-v1.*v1/vmaxlist(2)^2);
a2 = blist(3)*(1-v2.*v2/vmaxlist(3)^2);
acoast = -blist(4)*(1+vcoast.*vcoast/vmaxlist(4)^2);
aff = blist(5)*(1-vff.*vff/vmaxlist(5)^2);
a = [apre a1 a2 acoast aff];

end

