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
CL = 2*pi .* 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. 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–11 km): temperature drops at lapse rate L
% Stratosphere (11–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, 18, 300); % degrees
alpha_rad = alpha * pi/180;
% CL: thin aerofoil theory — CL increases at 2π per radian
CL = CL0 + 2*pi .* 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) / (2*pi) * 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 — R = (V/SFC) × (L/D) × ln(W₀/W₁) — is the fundamental formula connecting aircraft range to aerodynamic efficiency and fuel fraction. This script computes range, fuel burn, and then builds a contour map showing how range changes with every combination of L/D and SFC.
%% 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 ────────────────────────────────────────────────
R = (V / 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 ./ 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
eigs = eig(A);
fprintf('Open-loop eigenvalues:
');
for i = 1:4
wn = abs(eigs(i));
zeta = -real(eigs(i)) / wn;
fprintf(' λ%d = %+.4f%+.4fi wn=%.4f rad/s ζ=%.3f
', i, real(eigs(i)), imag(eigs(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, 1, 10, 100]); % heavily penalise pitch angle (theta)
R_lqr = 1; % elevator deflection cost
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 ──────────────────────────────────
theta_cmd = 5*pi/180; % 5 degree pitch up command (rad)
C = eye(4); % output all 4 states
sys_cl = ss(A_cl, B*K(4)*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 →