aircraft stability in matlab

main.m
% Define aircraft parameters
mass = 1000; % kg
CD0 = 0.05; % drag coefficient
CL_alpha = 0.1; % lift coefficient
gravity = 9.81; % m/s^2
wing_area = 20; % m^2
wing_span = 10; % m
tau = 0.1; % time constant for control input

% Define initial conditions
initial_altitude = 1000; % m
initial_velocity = 50; % m/s
initial_heading = 0; % degrees
initial_pitch = 0; % degrees

% Define simulation time and time step
t_end = 100; % seconds
dt = 0.01; % seconds
t = 0:dt:t_end;

% Define control input function
delta_e = @(t) 2*deg2rad(sin(0.1*t));

% Simulate aircraft model
altitude = zeros(size(t));
velocity = zeros(size(t));
heading = zeros(size(t));
pitch = zeros(size(t));
altitude(1) = initial_altitude;
velocity(1) = initial_velocity;
heading(1) = deg2rad(initial_heading);
pitch(1) = deg2rad(initial_pitch);
for i = 2:length(t)
    % Compute aerodynamic forces
    lift = 0.5*CL_alpha*wing_area*velocity(i-1)^2;
    drag = 0.5*CD0*wing_area*velocity(i-1)^2;
    
    % Compute gravitational force
    weight = mass*gravity;
    
    % Compute control input
    delta_e_i = delta_e(t(i-1));
    delta_theta_i = delta_e_i*exp(-dt/tau);
    
    % Compute aircraft motion
    altitude(i) = altitude(i-1) + velocity(i-1)*sin(pitch(i-1))*dt;
    velocity(i) = velocity(i-1) + (-drag-weight*sin(pitch(i-1))/mass)*dt;
    heading(i) = heading(i-1) + velocity(i-1)*cos(pitch(i-1))/altitude(i-1)*dt;
    pitch(i) = pitch(i-1) + delta_theta_i*dt;
end

% Plot aircraft motion
figure;
subplot(2,2,1);
plot(t,altitude);
xlabel('Time (s)');
ylabel('Altitude (m)');
subplot(2,2,2);
plot(t,velocity);
xlabel('Time (s)');
ylabel('Velocity (m/s)');
subplot(2,2,3);
plot(t,rad2deg(heading));
xlabel('Time (s)');
ylabel('Heading (deg)');
subplot(2,2,4);
plot(t,rad2deg(pitch));
xlabel('Time (s)');
ylabel('Pitch (deg)');
1804 chars
70 lines