Use
[t, vecste] = ode45(@(t,vecste)eom(t,vecste,spec,aero,thrust), tspan, stvecinit);
instead of
[t, vecste,spec,aero,thrust] = ode45('eom', tspan, stvecinit);
run('database.m')
% For Initial Values
uinit = 50;
vinit = 0;
winit = 1;
pinit = 0;
qinit = 0;
rinit = 0;
phiinit = 0;
thetainit = 0;
psiinit = 0;
xinit = 0;
yinit = 0;
H = 3000;
zinit = -H;
tspan = [0 3];
stvecinit = [uinit vinit winit pinit qinit rinit phiinit thetainit psiinit xinit yinit zinit];
[t, vecste] = ode45(@(t,vecste)eom(t,vecste,spec,aero,thrust), tspan, stvecinit);
%% Results Calculations
u = vecste(:,1);
v = vecste(:,2);
w = vecste(:,3);
p = vecste(:,4);
q = vecste(:,5);
r = vecste(:,6);
phi = vecste(:,7)*180/pi;
theta = vecste(:,8)*180/pi;
psi = vecste(:,9)*180/pi;
xx = vecste(:,10);
y = vecste(:,11);
z = vecste(:,12);
V = sqrt(u.^2+v.^2+w.^2);
alpha = atand(w./u);
beta = asin(v./V);
figure (1)
plot(t,alpha, 'LineWidth', 1.5)
grid on; grid minor
title('Variation of AoA with Time')
xlabel('Time (in sec)')
ylabel('AoA (in degrees)')
figure (2)
plot(t,beta, 'LineWidth', 1.5)
grid on; grid minor
title('Variation of pitch attitude with time')
xlabel('Time (in sec)')
ylabel('\theta (in degrees)')
figure (3)
plot(t,V, 'LineWidth', 1.5)
grid on; grid minor
title('Variation of Velocity with Time')
xlabel('Time (in sec)')
ylabel('Velocity (in m/sec)')
figure (4)
plot(xx,-z, 'LineWidth', 1.5)
grid on; grid minor
title('X vs Z')
xlabel('X (in m)')
ylabel('Z (in m)')
figure (5)
plot3(xx,y,-z, 'LineWidth', 1.5)
grid on; grid minor
title('3-D plot of the Trajectory')
xlabel('X (in m)')
ylabel('Y (in m)')
zlabel('Z (in m)')
function x = eom(t,vecste,spec,aero,thrust)
m = spec.m.table{1};
g = spec.g.table{1};
Ixx = spec.Ixx.table{1};
Iyy = spec.Iyy.table{1};
Izz = spec.Izz.table{1};
Ixz = spec.Ixz.table{1};
u = vecste(1);
v = vecste(2);
w = vecste(3);
p = vecste(4);
q = vecste(5);
r = vecste(6);
phi = vecste(7);
theta = vecste(8);
psi = vecste(9);
xx = vecste(10);
y = vecste(11);
z = vecste(12);
h = -z;
V = sqrt(u^2+v^2+w^2);
alpha = atan(w/u)
beta = asin(v/V);
V = sqrt(u^2+v^2+w^2);
[density,temperature,pressure] = atmosphere(h);
[Thrust] = thrustmodel(V, density,thrust);
[Fax,Fay,Faz,Mx,My,Mz] = aeromodel(p,q,r,V,alpha,beta,aero,spec,density);
[FGx,FGy,FGz] = weightmodel(phi,theta,psi,spec);
%Forces
X = Fax + FGx+Thrust;
Y = Fay + FGy;
Z = Faz + FGz;
% Moments already defined with same name using aeromodel function(Mx,My,Mz)
% 1. Translational Dynamics
udot = X/m-g*sin(theta)-q*w+r*v;
vdot = Y/m+g*sin(theta)*sin(phi)-r*u+p*w;
wdot = Z/m+g*cos(theta)*cos(phi)-p*v+q*u;
% 2. Rotional Dynamics
pdot = Mx/Ixx + ((Iyy-Izz) * (q*r/Ixx)); %+ (( (x(6,1) + p*q) /Ixx )*Ixz);
qdot = My/Iyy + ((Izz-Ixx) * (r*p/Iyy)); %- ((p^2 + r^2)/Iyy)*Ixz;
rdot = Mz/Izz + ((Ixx-Iyy) * (p*q/Izz)); %- ((q*r-x(4,1))/Izz)*Ixz;
% Euler dots
attmatrix = [ 1 tan(theta)*sin(phi) cos(phi)*tan(theta);
0 cos(phi) -sin(phi) ;
0 sec(theta)*sin(phi) sec(theta)*cos(phi)
];
omega = [p;q;r];
attvec = attmatrix*omega;
phidot = attvec(1);
thetadot = attvec(2);
psidot = attvec(3);
% x y z
vector = [u;v;w];
transvector = ned2body(phi,theta,psi);
accvector = transvector*vector;
xdot = accvector(1);
ydot = accvector(2);
zdot = accvector(3);
x(1,1) = udot;
x(2,1) = vdot;
x(3,1) = wdot;
x(4,1) = pdot;
x(5,1) = qdot;
x(6,1) = rdot;
x(7,1) = phidot;
x(8,1) = thetadot;
x(9,1) = psidot;
x(10,1)= xdot;
x(11,1)= ydot;
x(12,1)= zdot;
end
