EXPERIMENT - 1

Draw Pole-Zero map of dynamic system model with plot customization option

Aim: To plot Pole-Zero map of the given dynamic system model with an option of plot
customization option.

clc; Theory:
clear;In mathematics, signal processing and control theory, a pole–zero plot is a graphical
close representation
all; of a rational transfer function in the complex plane which helps to convey certain
properties of the system such as Stability, Causal system / anti causal system, Region of
% Define Transfer
convergence Function
(ROC), Minimum phase / non minimum phase etc.
% -------------------------------
The zeros of the system are roots of the numerator polynomial and the poles of the system
num = [1 2 -4]; % s^2 + 2s - 4
are roots of the denominator polynomial. A pole-zero plots shows the location in the complex
den = [1 6 25]; % s^2 + 6s + 25
plane of the poles and zeros of the transfer function of a dynamic system, such as a controller,
compensator, sensor, equalizer, filter, or communications channel. By convention, the poles of
% -------------------------------
the system are indicated in the plot by an X while the zeros are indicated by a circle or O.
% Find Zeros and Poles S2+2S−4
% -------------------------------
Example: Consider a dynamic system with a transfer function, H(s) =
zeros_val = roots(num); S2+6S+25
Here, the
poles_val zeros are solutions of equation S2 + 2S − 4 = 0 and the poles are the solutions of
= roots(den);
equation S2 + 6S + 25 = 0
% Plot Pole-Zero Map
figure;

plot(real(zeros_val), imag(zeros_val), 'o', 'MarkerSize', 10,

'DisplayName', 'Zeros');
hold on;

plot(real(poles_val), imag(poles_val), 'x', 'MarkerSize', 10,

'DisplayName', 'Poles');

% Axis lines
xline(0);
yline(0);

grid on;

xlabel('Real Axis');
ylabel('Imaginary Axis');
title('Pole-Zero Map');
legend;

xlim([-6 3]);
ylim([-6 6]);
 EXPERIMENT - 2

Plot root locus with variables in transfer function and

Dynamic system though MATLAB / PYTHON
clc;
clear;
close all;

% -------------------------------
% Open-loop transfer function
% G(s) = 1 / [s (s^2 + 3s + 2)]
% -------------------------------

% Numerator and Denominator

N = [1]; % Numerator
D = [1 3 2 0]; % Denominator: s^3 + 3s^2 + 2s

% -------------------------------
% Zeros and Poles
% -------------------------------
Z = roots(N);
P = roots(D);

% -------------------------------
% Gain variation
% -------------------------------
k = 0;
x = [];

for i = 1:200

% Characteristic equation:
% s^3 + 3s^2 + 2s + k = 0
C = [1 3 2 k];

roots_C = roots(C);

x = [x; roots_C.']; % store roots row-wise

k = k + 0.1;
end
% -------------------------------
% Root Locus Plot
% -------------------------------
figure;

plot(real(x), imag(x), '.', 'DisplayName', 'Root Locus');

hold on;

plot(real(Z), imag(Z), 'o', 'MarkerSize', 8, 'DisplayName', 'Zeros');

plot(real(P), imag(P), '*', 'MarkerSize', 10, 'DisplayName', 'Poles');

% Axis lines
xline(0);
yline(0);

title('Root Locus');
xlabel('Real Axis');
ylabel('Imaginary Axis');
grid on;
legend;
 Experiment -3

Aim: To draw Bode plot from a transfer function in MATLAB/PYTHON and

explain the gain and phase margins
clc;
clear;
close all;

% -------------------------------------------------
% Define transfer function H(s)
% H(s) = (37.5s + 75) / (s^3 + 16s^2 + 100s)
% -------------------------------------------------
num = [37.5 75];
den = [1 16 100 0];

sys = tf(num, den);

% -------------------------------------------------
% Frequency range
% -------------------------------------------------
w = logspace(-1, 3, 1000); % 0.1 to 1000 rad/s

% -------------------------------------------------
% Bode data
% -------------------------------------------------
[mag, phase, wout] = bode(sys, w);

% Convert to vectors
mag = squeeze(mag);
phase = squeeze(phase);

% Convert magnitude to dB
mag_db = 20*log10(mag);

% -------------------------------------------------
% Plot Bode Diagram
% -------------------------------------------------
figure;

subplot(2,1,1)
semilogx(wout, mag_db)
ylabel('Magnitude (dB)')
title('Bode Plot')
grid on
subplot(2,1,2)
semilogx(wout, phase)
xlabel('Frequency (rad/s)')
ylabel('Phase (deg)')
grid on

% -------------------------------------------------
% Gain Margin & Phase Margin Calculation
% -------------------------------------------------

% Gain crossover frequency (|G| = 0 dB)

[~, idx_gc] = min(abs(mag_db));
wg = wout(idx_gc);

% Phase crossover frequency (phase = -180 deg)

[~, idx_pc] = min(abs(phase + 180));
wp = wout(idx_pc);

% Gain Margin (dB)

GM = -mag_db(idx_pc);

% Phase Margin (deg)

PM = 180 + phase(idx_gc);

% -------------------------------------------------
% Display Results
% -------------------------------------------------
fprintf('Gain Margin (GM): %.2f dB\n', GM);
fprintf('Phase Margin (PM): %.2f degrees\n', PM);
fprintf('Gain crossover frequency (wg): %.2f rad/s\n', wg);
fprintf('Phase crossover frequency (wp): %.2f rad/s\n', wp);
 EXPERIMENT - 4

Simulate a spring- mass- damper system with and without a forcing

function though SIMULINK

Aim: To generate and analyze the response of a spring mass motion to a step and ramp input.

Theory:
clc;
clear;
close all;

% ------------------
% System Parameters
% ------------------
M = 2;
B = 2;
K = 2;

% ------------------
The Input
% Select Mathematical
Heremodel of the syste m is described by:
% ------------------ x = 1 (−Bx − Kx +f(t))
M
input_type = 'step'; % change to 'ramp' or 'none'
Let Mass M=2kg; Stiffness Constant K=2 N/m; B = 2 Ns/m. Then,
% ------------------ 1
% Time Span x = 2 (−2x − 2x + f(t))
% ------------------
tspan = [0 30];

% Initial Conditions [x, v]

y0 = [0 0];

% ------------------
% Solve ODE
% ------------------
[t, y] = ode45(@(t,y) model(t, y, M, B, K, input_type), tspan, y0);

% ------------------
% Plot
% ------------------
plot(t, y(:,1), 'LineWidth', 2)
title(['Spring Mass Damper - ', input_type, ' input'])
xlabel('Time (s)')
ylabel('Displacement')
grid on
% ------------------
% Function Definition
% ------------------
function dydt = model(t, y, M, B, K, input_type)

x = y(1); % displacement
v = y(2); % velocity

% ---- Input Force ----

if strcmp(input_type, 'step')
f = 1;
elseif strcmp(input_type, 'ramp')
f = t;
else
f = 0;
end

% ---- Equations ----

dxdt = v;
dvdt = (f - B*v - K*x) / M;

dydt = [dxdt; dvdt];

end
 EXPERIMENT - 5

Simulation of simple servo-mechanism feedback system in S domain

clc;
clear;
close all;

% --------------------------------------------
% System Parameters
% --------------------------------------------
A = 10;
K = 1;
B = 1;
J = 0.01;
F = 0.1;
R = 1;

% --------------------------------------------
% Transfer Function
% G(s) = KAB / (JR s^2 + FR s + KAB)
% --------------------------------------------
num = [K*A*B];
den = [J*R, F*R, K*A*B];

% Create transfer function

sys = tf(num, den);

% --------------------------------------------
% Step Response
% --------------------------------------------
t = 0:0.01:10; % time vector
[y, t] = step(sys, t);

% --------------------------------------------
% Plot
% --------------------------------------------
plot(t, y, 'LineWidth', 2)
title('Servo Mechanism Step Response')
xlabel('Time (seconds)')
ylabel('Output C(t)')
grid on
 EXPERIMENT - 6

Aim
To simulate the trajectory of a bomb dropped from an aircraft on a moving tank under pure
pursuit motion and to plot the bomb trajectory.

clc;
clear;
close all;

% --------------------------------
% Initial Conditions
% --------------------------------
x_aircraft = 0; % Aircraft initial x position (m)
y_aircraft = 50; % Aircraft height (m)

v_aircraft = 20; % Aircraft horizontal speed (m/s)

v_tank = 10; % Tank speed (m/s)

g = 9.81; % Gravity (m/s^2)

% --------------------------------
% Time of Flight Calculation
% y = y0 - (1/2)gt^2
% --------------------------------
t_flight = sqrt((2 * y_aircraft) / g);

% Time array
t = linspace(0, t_flight, 200);

% --------------------------------
% Bomb Motion (Projectile)
% --------------------------------
x_bomb = x_aircraft + v_aircraft * t;
y_bomb = y_aircraft - 0.5 * g * t.^2;

% --------------------------------
% Tank Motion
% --------------------------------
x_tank_initial = 20;
x_tank = x_tank_initial + v_tank * t;
y_tank = zeros(size(t));
% --------------------------------
% Plotting
% --------------------------------
figure;

plot(x_bomb, y_bomb, 'LineWidth', 2)

hold on

plot(x_tank, y_tank, 'LineWidth', 2)

scatter(x_bomb(end), 0, 50, 'r', 'filled')

xlabel('Horizontal Distance (m)')

ylabel('Vertical Height (m)')
title('Bomb Drop on Moving Tank (Pure Pursuit Motion)')
legend('Bomb Trajectory', 'Tank Motion', 'Impact Point')

grid on
EXPERIMENT - 7

Simulate aircraft Take-off and Landing with trajectory tracing

Aim: To simulate an aircraft Take-off and Landing with trajectory tracing

Take-off code
clc;
clear all;
close all;

% -------------------------------
% Initial Conditions
% -------------------------------
x = 0;
z = 0;

v0 = 0; % Initial velocity
vl = 64; % Lift-off velocity (m/s)
n = 1.68; % Load factor
a = 2.285; % Acceleration (m/s^2)
z1 = 50; % Target height (m)
g = 9.81;

% -------------------------------
% Ground Roll Calculations
% -------------------------------
sg = vl^2/(2*a); % Ground roll distance
tg = vl/a; % Time to lift-off

% -------------------------------
% Circular Climb
% -------------------------------
R = vl^2/(g*(n-1)); % Radius of curvature
w = g*(n-1)/vl; % Angular velocity

tr = (1/w)*acos(1-(z1/R)); % Time to reach height z1

xr = sg + R*sin(w*tr); % x at end of circular climb
theta = w*tr; % Climb angle

% -------------------------------
% Time Simulation
% -------------------------------
t = 0:1:60;

figure;
for i = 1:length(t)

v = v0 + a*t(i);

if (v <= vl) % Ground Roll Phase

 x(i) = v0*t(i) + 0.5*a*t(i)^2;
z(i) = 0;

else
if (z(i-1) <= z1) % Circular Climb
x(i) = sg + R*sin(w*(t(i)-tg));
z(i) = R*(1 - cos(w*(t(i)-tg)));
else % Straight Climb
x(i) = xr + (t(i)-(tr+tg))*vl*cos(theta);
z(i) = z1 + (t(i)-(tr+tg))*vl*sin(theta);
end
end

plot(x(i), z(i), 'r*')

axis([0 2000 0 2000])
grid on
drawnow
pause(0.05)
end

figure
plot(x, z, 'b', 'LineWidth', 2)
axis([0 2000 0 2000])
grid on
title('Take-Off Trajectory')
xlabel('Range (m)')
ylabel('Altitude (m)')

disp(['Ground Roll Distance sg = ', num2str(sg)])

Landing Code
clc;
clear all;
close all;

% -------------------------------
% Initial Conditions
% -------------------------------
z0 = 500; % Initial altitude
x0 = 0;
Va = 60; % Approach speed
g = 9.81;
n = 1.6;
a = -2.28; % Deceleration
z1 = 150; % Flare height
% -------------------------------
% Circular Flare Calculations
% -------------------------------
R = Va^2/(g*(n-1));
w = g*(n-1)/Va;

theta = acos(1 - z1/R);

ta = (z0 - z1)/(Va*sin(theta));
xa = x0 + ta*Va*cos(theta);

tr = theta/w;
xr = R*sin(theta);

tg = -Va/a;

% -------------------------------
% Time Simulation
% -------------------------------
t = 0:1:(ta + tr + tg);

figure;
for i = 1:length(t)

if (t(i) <= ta) % Glide Slope

x(i) = x0 + t(i)*Va*cos(theta);
z(i) = z0 - t(i)*Va*sin(theta);

elseif (t(i) <= ta + tr) % Circular Flare

x(i) = xa + R*(sin(theta) - sin(theta - w*(t(i)-ta)));
z(i) = R*(1 - cos(theta - w*(t(i)-ta)));

else % Ground Roll

x(i) = xa + xr + Va*(t(i)-(ta+tr)) + 0.5*a*(t(i)-(ta+tr))^2;
z(i) = 0;
end

plot(x(i), z(i), 'r*')

axis([0 2000 0 2000])
grid on
drawnow
pause(0.05)
end

figure
plot(x, z, 'b', 'LineWidth', 2)
axis([0 2000 0 2000])
grid on
title('Landing Trajectory')
xlabel('Range (m)')
ylabel('Altitude (m)')
EXPERIMENT-8

Simulate stall of aircraft and show the effect of variation in static margin on stalling
characteristics.

Aim: To simulate stall of aircraft and show the effect of variation in static margin on stalling
characteristics.

MATLAB CODE 1
(Chord Line Animation + Static Margin Calculation)

clc;
clear;
close all;

filename = 'EXP10_data.xlsx';
data = readmatrix(filename); % Updated function

n = size(data,1);

x1 = [-2.5 7.5];

figure;

subplot(2,2,1)

for j = 1:n

angle = data(j,1) * pi/180;

x = [x1(1)*cos(angle) x1(2)*cos(angle)];
y = [-x1(1)*sin(angle) -x1(2)*sin(angle)];

plot(x,y,'b','LineWidth',2)
title('Chord Line vs AoA')
axis([-10 10 -5 5])
grid on
drawnow
pause(0.3)
end

% Static Margin Calculation

for i = 1:n-1
 Claplha(i) = (data(i+1,2)-data(i,2)) / ...
(data(i+1,1)-data(i,1));

Cmalpha(i) = (data(i+1,3)-data(i,3)) / ...

(data(i+1,1)-data(i,1));

SM(i) = -Cmalpha(i)/Claplha(i);
end

SM(n) = SM(n-1);

subplot(2,2,2)
plot(data(:,1),SM,'r','LineWidth',2)
xlabel('AoA')
ylabel('Static Margin')
grid on

subplot(2,2,3)
plot(data(:,1),data(:,2),'g','LineWidth',2)
xlabel('AoA')
ylabel('Lift Coefficient')
grid on

subplot(2,2,4)
plot(data(:,1),data(:,3),'k','LineWidth',2)
xlabel('AoA')
ylabel('Pitching Moment')
grid on

MATLAB CODE 2
(Static Margin and Stall Characteristics)

clc;
clear all;
close all;

% ---------------------------------------
% Aircraft Parameters
% ---------------------------------------
mass = 5000; % kg
wingArea = 25; % m^2
wingSpan = 10; % m
wingChord = wingArea / wingSpan;
incidenceAngle = 0; % degrees
liftCurveSlope = 0.1; % per degree
cgPosition = 0.3; % fraction of chord

% ---------------------------------------
% Atmospheric Conditions
% ---------------------------------------
rho = 1.225; % kg/m^3
g = 9.81; % m/s^2
weight = mass * g;

% ---------------------------------------
% Simulation Parameters
% ---------------------------------------
alpha = -10:0.5:20; % AoA range
staticMarginRange = -0.2:0.05:0.2;

% ---------------------------------------
% Stall Detection
% ---------------------------------------
for i = 1:length(staticMarginRange)

staticMargin = staticMarginRange(i) * wingChord;

neutralPoint = cgPosition * wingChord + staticMargin;

for j = 1:length(alpha)

aerodynamicCenter = neutralPoint - ...

(liftCurveSlope / (2*pi)) * wingChord;

cm = -(cgPosition*wingChord - aerodynamicCenter) ...

* liftCurveSlope * (alpha(j) - incidenceAngle);

cl = liftCurveSlope * (alpha(j) - incidenceAngle);

lift = 0.5 * rho * wingArea * cl * (50^2); % assumed V = 50 m/s

if lift < weight

fprintf('Stall at alpha = %.2f deg, Static Margin = %.2f m\n', ...
alpha(j), staticMargin);
break
end
end
end

% ---------------------------------------
% Static Margin vs Stall Speed
% ---------------------------------------
wingArea = 20; % m^2
meanAerodynamicChord = 2.5; % m
weight = 5000 * g; % N

staticMargin = -0.1:0.05:0.1;

stallSpeed = sqrt((2 * weight) ./ ...

(rho * wingArea .* ...
(1 + staticMargin)));
figure
plot(staticMargin, stallSpeed, 'b-o', 'LineWidth', 2)
xlabel('Static Margin')
ylabel('Stall Speed (m/s)')
title('Effect of Static Margin on Stall Speed')
grid on