minimumClimbDae(soldae,iphase);

function daeout = minimumClimbDae(soldae,iphase);
global CoF;

t = soldae{1};
s = soldae{2};
u = soldae{3};
p = soldae{4};
h = s(:,1);
E = s(:,2);
fpa = s(:,3);

g = 9.80665;
m = 37000./2.2;
S = 60;
hbar = h./1000;


%rho calculation
z = (-3.48643241e-2).*hbar+(3.50991865e-3).*hbar.^2-(8.33000535e-5).*hbar.^3+(1.15219733e-6).*hbar.^4;
r = 1.0228066.*exp(-z);
y = -0.12122693.*hbar+r-1.0228055;
rho = 1.225.*exp(y);
%rho calculation

%speed of sound%
theta = 292.1-8.87743.*hbar+0.193315.*hbar.^2+(3.72e-3).*hbar.^3;
a     = 20.0468.*sqrt(theta);
%speed of sound%

%Velocity and mach
Eminush = (E-h).^2;
Eminush = sqrt(Eminush);
v = sqrt(2.*g.*Eminush);

M = v./a;
%Velocity and mach

%Who are lift and drag
q = 0.5.*rho.*v.*v.*S;
L = m.*g.*u;
numeratorCD0 = CoF(1,1).*M.^0+CoF(1,2).*M.^1+CoF(1,3).*M.^2+CoF(1,4).*M.^3+CoF(1,5).*M.^4;
denominatorCD0 = CoF(2,1).*M.^0+CoF(2,2).*M.^1+CoF(2,3).*M.^2+CoF(2,4).*M.^3+CoF(2,5).*M.^4;
Cd0 = numeratorCD0./denominatorCD0;
numeratorK = CoF(3,1).*M.^0+CoF(3,2).*M.^1+CoF(3,3).*M.^2+CoF(3,4).*M.^3+CoF(3,5).*M.^4;
denominatorK = CoF(4,1).*M.^0+CoF(4,2).*M.^1+CoF(4,3).*M.^2+CoF(4,4).*M.^3+CoF(4,5).*M.^4+CoF(4,6).*M.^5;
K   = numeratorK./denominatorK;
D = q.*(Cd0+K.*((m.^2).*(g.^2)./(q.^2)).*(u.^2));
%Who are lift and drag

%Who is thrust
e0 = CoF(5,1).*M.^0+CoF(6,1).*M.^1+CoF(7,1).*M.^2+CoF(8,1).*M.^3+CoF(9,1).*M.^4+CoF(10,1).*M.^5;
e1 = CoF(5,2).*M.^0+CoF(6,2).*M.^1+CoF(7,2).*M.^2+CoF(8,2).*M.^3+CoF(9,2).*M.^4+CoF(10,2).*M.^5;
e2 = CoF(5,3).*M.^0+CoF(6,3).*M.^1+CoF(7,3).*M.^2+CoF(8,3).*M.^3+CoF(9,3).*M.^4+CoF(10,3).*M.^5;
e3 = CoF(5,4).*M.^0+CoF(6,4).*M.^1+CoF(7,4).*M.^2+CoF(8,4).*M.^3+CoF(9,4).*M.^4+CoF(10,4).*M.^5;
e4 = CoF(5,5).*M.^0+CoF(6,5).*M.^1+CoF(7,5).*M.^2+CoF(8,5).*M.^3+CoF(9,5).*M.^4+CoF(10,5).*M.^5;
e5 = CoF(5,6).*M.^0+CoF(6,6).*M.^1+CoF(7,6).*M.^2+CoF(8,6).*M.^3+CoF(9,6).*M.^4+CoF(10,6).*M.^5;

T = (e0.*hbar.^0+e1.*hbar.^1+e2.*hbar.^2+e3.*hbar.^3+e4.*hbar.^4+e5.*hbar.^5).*9.80665./2.2;
%Who is thrust

%************************************************

%User Input - Give me f in dot(x) = f
%*********************

Edot = v./(m.*g).*(T-D);
hdot = v.*sin(fpa);
fpadot = g./v.*(L./(m.*g)-cos(fpa));


daeout = [hdot Edot fpadot];