MATLAB is the language aerospace engineers actually use. Every example on this page is complete, working, copy-and-run code covering ISA atmosphere, lift and drag, orbit propagation, FFT analysis, and LQR flight control.
MATLAB is not complex. It has a small number of rules that govern all its behaviour. Master these five and every script you ever write will make sense.
A scalar is a 1×1 matrix. A vector is a 1×n matrix. An array of 1,000 values is a 1×1000 matrix. This consistency means every operation (arithmetic, indexing, functions) works identically regardless of size. Writing a calculation for one altitude value or 10,000 altitude values is the same code.
a = 5 creates a 1×1 matrix.
v = [10 20 30] creates a 1×3 row vector.
h = linspace(0,20000,500) creates 500 equally-spaced altitudes in one command.
A * B is matrix multiplication. Dimensions must be compatible (n×m times m×p). Used for linear algebra.
A .* B is element-wise multiplication. Each element multiplied by the corresponding element. Used for array calculations: CL at every angle of attack, density at every altitude, thrust at every time step.
Same rule applies to .^ (element-wise power) and ./ (element-wise divide). When working with arrays, you almost always want the dot version. Forgetting the dot is the single most common MATLAB error.
x = 3, assigns 3 to x AND prints x = 3 to the command window.
x = 3;, assigns 3 to x silently. No output.
In scripts computing thousands of values, always use semicolons or your command window floods with numbers. Temporarily remove a semicolon to inspect an intermediate value when debugging. Use fprintf to print specific formatted results.
First element is v(1), not v(0). Last element is v(end). Slice with colon: v(2:5) gives elements 2 through 5. v(end-2:end) gives the last three.
Everything after % on a line is a comment, ignored by MATLAB. Use %% at the start of a line to create a section break (runs as a cell block in the editor).
The most powerful functions you will use daily: linspace, zeros, ones, plot, fprintf, max, min, find, ode45, fft, eig, and \ (backslash linear solver).
%% MATLAB SYNTAX QUICK REFERENCE, copy this into MATLAB and run it
% ── Creating arrays ──────────────────────────────────────────────────
alpha = linspace(-5, 20, 100); % 100 evenly spaced AoA values (degrees)
h = 0:100:10000; % altitudes from 0 to 10,000 m in 100 m steps
M = zeros(3,3); % 3×3 matrix of zeros
I = eye(4); % 4×4 identity matrix
% ── Element-wise operations (the dot operators) ──────────────────────
alpha_r = alpha * pi/180; % scalar * vector: fine without dot
a3d = 2*pi/(1+2/(8*0.85)); % finite-wing slope, AR 8, e 0.85
CL = a3d .* alpha_r; % .* element-wise: CL for each AoA
CD = 0.012 + CL.^2 / (pi*8*0.85); % .^ element-wise power
LD = CL ./ CD; % ./ element-wise divide
% ── Indexing ─────────────────────────────────────────────────────────
first = alpha(1); % first element: MATLAB starts at 1, NOT 0
last = alpha(end); % last element
mid = alpha(20:40); % elements 20 through 40
tail = alpha(end-4:end); % last 5 elements
% ── Finding values ───────────────────────────────────────────────────
[LD_max, idx] = max(LD); % max value AND its index in one call
fprintf('Max L/D = %.2f at alpha = %.1f deg
', LD_max, alpha(idx));
% ── Logical indexing ─────────────────────────────────────────────────
good_alpha = alpha(LD > 12); % all AoA values where L/D > 12
% ── Solving linear systems ───────────────────────────────────────────
A = [2 1; 5 3]; b = [8; 19];
x = A \ b; % backslash = MATLAB's linear solver (Ax = b)
fprintf('x1=%.1f, x2=%.1f
', x(1), x(2));
Every script below solves a real aerospace engineering problem. All code is complete and runnable, copy it into MATLAB, press Run, and it works.
The International Standard Atmosphere defines temperature, pressure, density, and speed of sound at every altitude. One subtlety this script glosses over: ISA is defined on geopotential altitude H = REz/(RE+z), not geometric altitude z. Treating them as equal costs about 0.3% in pressure at 11 km and more above, which is why the values here sit slightly high against published ISA tables. Every performance calculation (stall speed, climb rate, Mach number, engine thrust) depends on these values. MATLAB has a built-in atmosisa() function but this script builds it from first principles so you understand every line.
%% ISA STANDARD ATMOSPHERE: from scratch, 0 to 25 km
% ── Physical constants ──────────────────────────────────────────────
T0 = 288.15; % Sea-level temperature (K)
P0 = 101325; % Sea-level pressure (Pa)
L = 0.0065; % Temperature lapse rate in troposphere (K/m)
g = 9.80665; % Gravity (m/s²)
R = 287.058; % Specific gas constant for air (J/kg·K)
gam = 1.4; % Ratio of specific heats
% ── Altitude array: 0 to 25,000 m ──────────────────────────────────
h = linspace(0, 25000, 1000);
% ── Temperature ─────────────────────────────────────────────────────
% Troposphere (0 to 11 km): temperature drops at lapse rate L
% Stratosphere (11 to 25 km): isothermal at 216.65 K
T = T0 - L .* h;
T(h > 11000) = 216.65;
% ── Pressure ─────────────────────────────────────────────────────────
% Troposphere: hydrostatic equation gives a power-law profile
P = P0 .* (T / T0).^(g / (L * R));
% Stratosphere: exponential decay (isothermal, constant scale height)
i11 = find(h > 11000, 1);
P(h > 11000) = P(i11) .* exp(-g * (h(h > 11000) - 11000) / (R * 216.65));
% ── Density and speed of sound ───────────────────────────────────────
rho = P ./ (R .* T); % Ideal gas: rho = P / (R·T)
a = sqrt(gam * R .* T); % Speed of sound: a = sqrt(γRT)
% ── Print values at key altitudes ────────────────────────────────────
targets = [0, 1524, 5000, 10668, 11000, 20000]; % 0, 5000ft, 16404ft, 35000ft
fprintf('%-10s %-8s %-10s %-12s %-9s
','Alt(m)','T(K)','P(Pa)','rho(kg/m3)','a(m/s)');
fprintf('%s
', repmat('-',1,54));
for alt = targets
[~,i] = min(abs(h - alt));
fprintf('%-10g %-8.2f %-10.1f %-12.4f %-9.1f
', alt,T(i),P(i),rho(i),a(i));
end
% ── 4-panel plot ─────────────────────────────────────────────────────
figure('Name','ISA Atmosphere');
subplot(2,2,1); plot(T,h/1000,'b-','LineWidth',2); xlabel('T (K)'); ylabel('Alt (km)'); grid on;
subplot(2,2,2); plot(P/1000,h/1000,'r-','LineWidth',2); xlabel('P (kPa)'); ylabel('Alt (km)'); grid on;
subplot(2,2,3); plot(rho,h/1000,'g-','LineWidth',2); xlabel('rho (kg/m³)'); ylabel('Alt (km)'); grid on;
subplot(2,2,4); plot(a,h/1000,'m-','LineWidth',2); xlabel('a (m/s)'); ylabel('Alt (km)'); grid on;
Computes CL, CD, L/D ratio, and actual aerodynamic forces across a full angle-of-attack sweep. Uses the parabolic drag polar. The equation every aircraft performance module is built on. The script finds cruise AoA, maximum L/D, and minimum drag speed.
%% LIFT, DRAG & PERFORMANCE. A320-like aircraft at cruise
% ── Aircraft and wing parameters ─────────────────────────────────────
S = 122.4; % Wing reference area (m²)
AR = 9.5; % Aspect ratio
e = 0.82; % Oswald span efficiency factor
CD0 = 0.024; % Zero-lift drag coefficient (parasite drag)
CL0 = 0.28; % CL at zero AoA (from wing camber)
% ── Cruise flight conditions ─────────────────────────────────────────
rho = 0.3796; % Air density at 35,000 ft from ISA (kg/m³)
V_cruise = 232; % True airspeed (m/s) ≈ M0.78 at cruise
W = 68000 * 9.81; % Aircraft weight mid-cruise (N)
q = 0.5 * rho * V_cruise^2; % Dynamic pressure (Pa)
% ── AoA sweep ────────────────────────────────────────────────────────
alpha = linspace(-4, 14, 300); % degrees - capped at stall
% The linear model has no stall. Sweeping past ~14 deg would report CL and
% L/D values that do not exist on a real wing, so the sweep stops there.
% A model is only worth what its validity range is worth.
alpha_rad = alpha * pi/180;
% CL: thin-aerofoil theory gives the TWO-DIMENSIONAL slope a0 = 2*pi/rad.
% A finite wing is always shallower, because trailing vorticity induces
% downwash that cuts the effective incidence (Prandtl lifting line):
% a = a0 / (1 + a0/(pi*AR*e))
% At AR 9.5 this is 0.0873/deg against the 2D value of 0.1097/deg, so
% using 2*pi on a wing over-predicts lift by about 20%.
a0 = 2*pi; % 2D lift slope (per radian)
a = a0 / (1 + a0/(pi*AR*e)); % 3D lift slope (per radian)
CL = CL0 + a .* alpha_rad;
% CD: parabolic drag polar: CD = CD0 + CL²/(π·AR·e)
CD = CD0 + CL.^2 ./ (pi * AR * e);
% Actual forces (N)
Lift = q * S .* CL;
Drag = q * S .* CD;
LD = CL ./ CD;
% ── Find cruise AoA (where Lift = Weight) ────────────────────────────
CL_cruise = W / (q * S);
alpha_cruise = (CL_cruise - CL0) / a * 180/pi;
CD_cruise = CD0 + CL_cruise^2 / (pi*AR*e);
LD_cruise = CL_cruise / CD_cruise;
Drag_cruise = q * S * CD_cruise;
% ── Maximum L/D ──────────────────────────────────────────────────────
[LD_max, imax] = max(LD);
fprintf('--- Cruise Results ---
');
fprintf('CL_cruise = %.4f
', CL_cruise);
fprintf('CD_cruise = %.5f
', CD_cruise);
fprintf('L/D at cruise = %.2f
', LD_cruise);
fprintf('Cruise AoA = %.2f deg
', alpha_cruise);
fprintf('Drag force = %.0f N
', Drag_cruise);
fprintf('Max L/D = %.2f at alpha = %.2f deg
', LD_max, alpha(imax));
The Breguet equation connects aircraft range to aerodynamic efficiency and fuel fraction. Watch the units: if SFC is a mass-specific consumption in kg/(N·s), the equation is R = V·(L/D)·ln(W₀/W₁) / (g₀·SFC), because the weight flow is g₀·SFC·T. If SFC is quoted as a weight-specific consumption in 1/s, g₀ drops out. Getting this wrong inflates the answer by a factor of 9.81, which is exactly the kind of error a sanity check catches: an A320 burning 11.5 t of fuel does not fly 38,000 km. Assumptions: constant L/D, constant SFC, and cruise-climb at constant CL. Real cruise is stepped, so this is an upper bound before climb, descent, reserves and routing are deducted.
%% BREGUET RANGE EQUATION + PARAMETRIC STUDY
% R = (V / SFC) * (L/D) * ln(W0/W1)
% ── Aircraft parameters ──────────────────────────────────────────────
V = 232; % Cruise true airspeed (m/s) at M0.78 / 35,000 ft
SFC = 1.5e-5; % Thrust-specific fuel consumption (kg/N/s), turbofan
LD = 14.5; % Lift-to-drag ratio at cruise
W0 = 73500 * 9.81; % Max take-off weight (N)
W1 = 62000 * 9.81; % Landing weight (N)
% ── Range calculation ────────────────────────────────────────────────
% Weight-flow, not mass-flow. With SFC in kg/(N·s) the fuel WEIGHT
% flow is g0*SFC*T, so g0 must appear in the denominator. Dropping it
% inflates the range by a factor of 9.81.
g0 = 9.80665; % standard gravity (m/s^2)
R = (V / (g0 * SFC)) * LD * log(W0 / W1); % log = natural log in MATLAB
R_km = R / 1000;
R_nm = R / 1852;
W_fuel = (W0 - W1) / 9.81; % Fuel burned (kg)
fuel_frac = W_fuel / (W0/9.81) * 100; % Fuel fraction (%)
fprintf('Range = %.0f km (%.0f nm)
', R_km, R_nm);
fprintf('Fuel burned = %.0f kg (%.1f%% of MTOW)
', W_fuel, fuel_frac);
% ── Parametric study: Range vs L/D and SFC ───────────────────────────
LD_vec = linspace(8, 20, 60);
SFC_vec = linspace(0.8e-5, 2.5e-5, 60);
% meshgrid creates 2D grids from two 1D vectors
[LD_g, SFC_g] = meshgrid(LD_vec, SFC_vec);
% Vectorised Breguet over every combination
R_g = (V ./ (g0*SFC_g)) .* LD_g .* log(W0/W1) / 1000; % km
figure('Name','Breguet Parametric Study');
contourf(LD_g, SFC_g*1e5, R_g, 20); % 20 contour levels
colorbar; xlabel('L/D ratio'); ylabel('SFC (× 10⁻⁵ kg/N/s)');
title('Breguet Range (km)'); hold on;
plot(LD, SFC*1e5, 'w*', 'MarkerSize',14); % mark our aircraft
text(LD+0.3, SFC*1e5, ' Our aircraft', 'Color','white');
ode45 is MATLAB's built-in Runge-Kutta integrator, arguably the most useful aerospace function that exists. You provide the equations of motion and initial conditions; it integrates them numerically through time. Used for orbit propagation, rocket trajectories, aircraft 6-DOF simulation: anything involving differential equations.
%% TWO-BODY ORBIT PROPAGATION: ISS orbit, 3 revolutions
% ode45 integrates the equations of motion: [r_dot; v_dot] = [v; a]
mu = 3.986004418e14; % Earth gravitational parameter μ = GM (m³/s²)
RE = 6.371e6; % Earth mean radius (m)
% ── ISS initial conditions ────────────────────────────────────────────
alt = 408e3; % Orbital altitude (m)
r0 = RE + alt; % Orbital radius (m)
v0 = sqrt(mu / r0); % Circular orbit speed from vis-viva: v = sqrt(μ/r)
T = 2*pi * sqrt(r0^3/mu); % Orbital period (s)
% State vector: [x, y, vx, vy]: working in 2D orbital plane
y0 = [r0; 0; 0; v0]; % Start on the x-axis, velocity in y-direction
% ── Equations of motion (Newton's law of gravitation) ─────────────────
% State derivative: dx/dt = [vx; vy; ax; ay]
% Acceleration: a = -μ * r_vec / |r|³
eom = @(t,y) [y(3);
y(4);
-mu*y(1) / (y(1)^2 + y(2)^2)^(3/2);
-mu*y(2) / (y(1)^2 + y(2)^2)^(3/2)];
% ── Integrate for 3 complete orbits ──────────────────────────────────
% Tight tolerances are essential for orbital mechanics accuracy
opts = odeset('RelTol',1e-12, 'AbsTol',1e-9);
[t, Y] = ode45(eom, [0, 3*T], y0, opts);
% ── Results ──────────────────────────────────────────────────────────
r_error = norm(Y(end,1:2) - y0(1:2)'); % closure error after 3 orbits
fprintf('Orbital velocity = %.2f m/s
', v0);
fprintf('Orbital period = %.3f min
', T/60);
fprintf('Orbit closure err = %.4f m (3 revolutions)
', r_error);
% ── Plot orbit ───────────────────────────────────────────────────────
figure('Name','ISS Orbit');
plot(Y(:,1)/1e6, Y(:,2)/1e6, 'b-', 'LineWidth',1.5); hold on;
theta = linspace(0,2*pi,300);
fill(RE/1e6*cos(theta), RE/1e6*sin(theta), [0.2 0.3 0.8]); % Earth
axis equal; grid on; xlabel('x (Mm)'); ylabel('y (Mm)');
title(sprintf('ISS Orbit: %.0f km altitude, T = %.1f min', alt/1000, T/60));
The Fast Fourier Transform converts a time-domain signal into its frequency components. Essential for: finding structural resonance frequencies, processing accelerometer data from rockets or UAVs, diagnosing engine vibration, and validating sensor data quality.
%% FFT: VIBRATION FREQUENCY ANALYSIS
% Simulates accelerometer data from a rocket and finds resonance frequencies
% ── Simulate sensor data ──────────────────────────────────────────────
fs = 1000; % Sample rate (Hz), typical IMU rate
t = 0 : 1/fs : 3; % 3 seconds of data
N = length(t);
% Signal = structural bending mode + engine vibration + sensor noise
f1 = 12; A1 = 1.2; % 12 Hz structural mode, amplitude 1.2 g
f2 = 62; A2 = 0.4; % 62 Hz engine vibration, amplitude 0.4 g
f3 = 187; A3 = 0.15; % 187 Hz higher harmonic
noise_amp = 0.1; % White noise amplitude
signal = A1*sin(2*pi*f1*t) + A2*sin(2*pi*f2*t) + A3*sin(2*pi*f3*t) + noise_amp*randn(1,N);
% ── Compute FFT ───────────────────────────────────────────────────────
Y = fft(signal); % complex frequency spectrum
P2 = abs(Y/N); % two-sided amplitude spectrum
P1 = P2(1 : floor(N/2)+1); % single-sided (positive frequencies only)
P1(2:end-1) = 2*P1(2:end-1); % correct amplitude (halved by single-siding)
f = fs * (0 : floor(N/2)) / N; % frequency axis (Hz)
% ── Find peaks ───────────────────────────────────────────────────────
[pks, locs] = findpeaks(P1, 'MinPeakHeight',0.05, 'SortStr','descend');
fprintf('Dominant frequencies:
');
for k = 1:min(5,length(locs))
fprintf(' %6.1f Hz | amplitude = %.3f g
', f(locs(k)), pks(k));
end
% ── Plot time and frequency domain ───────────────────────────────────
figure('Name','Vibration FFT');
subplot(2,1,1);
plot(t(1:500), signal(1:500));
xlabel('Time (s)'); ylabel('Acceleration (g)'); grid on; title('Time Domain');
subplot(2,1,2);
plot(f, P1, 'b-', 'LineWidth',1.5); xlim([0 250]);
xlabel('Frequency (Hz)'); ylabel('Amplitude (g)'); grid on; title('Frequency Domain (FFT)');
State-space representation of aircraft longitudinal dynamics, the mathematics behind fly-by-wire autopilots. This script builds the A and B matrices from stability derivatives, analyses open-loop modes via eigenvalues, designs an LQR optimal controller, and simulates the step response to a pitch command.
%% LONGITUDINAL DYNAMICS: STATE SPACE + LQR AUTOPILOT DESIGN
%
% States: x = [u, w, q, theta]
% u = forward speed perturbation (m/s)
% w = heave velocity perturbation (m/s)
% q = pitch rate (rad/s)
% theta = pitch angle perturbation (rad)
%
% Control input: delta_e = elevator deflection (rad)
% dx/dt = A*x + B*u
V0 = 200; % Trim speed (m/s)
g = 9.81;
% ── A matrix (longitudinal stability derivatives) ─────────────────────
A = [-0.0188, 0.0482, 0, -g; % u-dot
-0.5700, -1.1000, V0, 0; % w-dot
0.0014, -0.0488, -0.527, 0; % q-dot
0, 0, 1, 0]; % theta-dot = q
% ── B matrix (elevator control effectiveness) ─────────────────────────
B = [0; -7.6; -15.7; 0];
% ── Open-loop eigenvalue analysis ────────────────────────────────────
% Eigenvalues of A reveal natural flight modes
% Negative real parts = stable mode
lam = eig(A); % do not name this 'eigs' - that shadows a built-in
fprintf('Open-loop eigenvalues:
');
for i = 1:4
wn = abs(lam(i));
zeta = -real(lam(i)) / wn;
fprintf(' λ%d = %+.4f%+.4fi wn=%.4f rad/s ζ=%.3f
', i, real(lam(i)), imag(lam(i)), wn, zeta);
end
% ── LQR controller design ─────────────────────────────────────────────
% Q weights state deviations, R weights control effort
% Tune Q to make the controller prioritise different states
Q = diag([0.01, 0.1, 1, 50]); % state weights
R_lqr = 2000; % elevator cost - see note below
% Tuning matters physically, not just numerically. With R = 1 this problem
% returns closed-loop poles near 47 rad/s and demands roughly 600 deg of
% elevator for a 5 deg pitch command. That is not a controller, it is a
% number. Raise R until the bandwidth and the control deflection are things
% a real actuator could deliver, then check the damping against MIL-STD-1797.
K = lqr(A, B, Q, R_lqr); % optimal feedback gain [1×4 vector]
fprintf('
LQR gain K = [%.4f %.4f %.4f %.4f]
', K);
% Closed-loop: A_cl = A - B*K
A_cl = A - B*K;
fprintf('Closed-loop eigenvalues:
');
disp(eig(A_cl)) % all negative real parts = stable closed loop
% ── Step response: 5° pitch command ──────────────────────────────────
% State feedback alone does NOT track a command: u = -Kx drives the states
% to zero, not to a setpoint. To make theta follow a 5 deg command you need
% a feedforward gain Nbar chosen so the closed-loop DC gain to theta is 1.
% Scaling B by K(4) is not that, and gives no guarantee about steady state.
theta_cmd = 5*pi/180; % 5 degree pitch up command (rad)
C_theta = [0 0 0 1]; % we care about theta only
Nbar = 1/dcgain(ss(A_cl, B, C_theta, 0));
C = eye(4); % output all 4 states
sys_cl = ss(A_cl, B*Nbar*theta_cmd, C, zeros(4,1));
figure('Name','LQR Step Response');
step(sys_cl(4,1)); % plot pitch angle (state 4) response
ylabel('Pitch angle (rad)'); title('LQR: 5° Pitch Step Response');
Import Cp distributions from ANSYS Fluent into MATLAB, integrate for CL and CD, and plot validation against NACA data.
Log sensor data with Arduino, then import the CSV files into MATLAB and apply FFT, filtering, and statistical analysis.
The stability derivatives, phugoid mode, and LQR control design in the code examples come directly from flight mechanics theory.
Propagate orbits with ode45, implement Tsiolkovsky, and simulate Hohmann transfers: all covered in the orbit propagation example.
Orbit propagation studies, gas turbine cycle optimisation, fatigue crack growth: all implemented in MATLAB with full code.
Understand the engineering context behind every MATLAB script, what the equations mean in the real world.
SheCodes Lab teaches Python and C++ from scratch, side by side, free, no experience needed. Includes an engineering module covering NumPy, pandas, ISA models, cost index, and flight data analysis. The same tools used to build the calculators on this site.
shecodeslab.com →