launchDae(sol,iphase);

function dae = launchDae(sol,iphase);

global CONSTANTS;
t = sol{1};
x = sol{2};
u = sol{3};
p = sol{4};
r = x(:,1:3);
v = x(:,4:6);
m = x(:,7);
zz = zeros(size(t));
rad = sqrt(sum(r.*r,2));
omega_matrix = CONSTANTS.omega_matrix;
% omega_matrix = [0 -CONSTANTS.omega 0; CONSTANTS.omega 0 0; 0 0 0];
vrel = v-r*omega_matrix.';
speedrel = sqrt(sum(vrel.*vrel,2));
altitude = rad-CONSTANTS.Re;
rho = CONSTANTS.rho0*exp(-altitude/CONSTANTS.H);
a1 = CONSTANTS.rho0*CONSTANTS.sa*CONSTANTS.cd;
a2 = a1*exp(-altitude/CONSTANTS.H);
bc = a2.*(rho./(2*m));
bc  = (rho./(2*m)).*CONSTANTS.sa*CONSTANTS.cd;

bcspeedmat = [zz zz zz];
bcspeed = bc.*speedrel;
bcspeedmat = [bcspeed bcspeed bcspeed];
Drag = [zz zz zz];
Drag = -bcspeedmat.*vrel;
muoverradcubedmat = [zz zz zz];
muoverradcubed = CONSTANTS.mu./rad.^3;
muoverradcubedmat = [muoverradcubed muoverradcubed muoverradcubed];
grav = [zz zz zz];
grav = -muoverradcubedmat.*r;

if iphase==1,
    T_srb = 6*CONSTANTS.thrust_srb*ones(size(t));
    T_first = CONSTANTS.thrust_first*ones(size(t));
    T_tot = T_srb+T_first;
    m1dot = -T_srb./(CONSTANTS.g0*CONSTANTS.ISP_srb);
    m2dot = -T_first./(CONSTANTS.g0*CONSTANTS.ISP_first);
    mdot = m1dot+m2dot;
elseif iphase==2,
    T_srb = 3*CONSTANTS.thrust_srb*ones(size(t));
    T_first = CONSTANTS.thrust_first*ones(size(t));
    T_tot = T_srb+T_first;
    m1dot = -T_srb./(CONSTANTS.g0*CONSTANTS.ISP_srb);
    m2dot = -T_first./(CONSTANTS.g0*CONSTANTS.ISP_first);
    mdot = m1dot+m2dot;    
elseif iphase==3
    T_first = CONSTANTS.thrust_first*ones(size(t));
    T_tot = T_first;
    mdot = -T_first./(CONSTANTS.g0*CONSTANTS.ISP_first);
elseif iphase==4,
    T_second = CONSTANTS.thrust_second*ones(size(t));
    T_tot = T_second;
    mdot = -T_second./(CONSTANTS.g0*CONSTANTS.ISP_second);
end;

path = sum(u.*u,2);
Tovermmat = [zz zz zz];
Toverm = T_tot./m;
Tovermmat = [Toverm Toverm Toverm];
thrust = Tovermmat.*u;

rdot = v;
vdot = thrust+Drag+grav;

dae(:,8) = path;
dae(:,1:7) = [rdot vdot mdot];