Intro to Programming for Engineers

Example 1.1 Compute the weight of a 1,100 kg aircraft and check the types.

mass_kg = 1100
g = 9.81
weight = mass_kg * g
print(weight)
print(type(mass_kg), type(weight))

10791.0
<class 'int'> <class 'float'>

Example 1.2 How much has Mauna Loa CO₂ risen from the 1959 to the 2025 annual mean (NOAA)?

ppm_1959 = 315.98
ppm_2025 = 427.35
rise = ppm_2025 - ppm_1959
print(f"CO2 rose {rise:.2f} ppm, or {rise / ppm_1959:.1%}")

CO2 rose 111.37 ppm, or 35.2%

Example 1.3 Explore operators and float surprises.

print(7 / 2, 7 // 2, 7 % 2, 2 ** 10)
print(0.1 + 0.2 == 0.3, abs((0.1 + 0.2) - 0.3) < 1e-9)
print(int("42") + 1, str(42) + "1")

3.5 3 1 1024
False True
43 421

Example 2.1 Classify NASA’s 2024 global temperature anomaly.

anomaly = 1.28   # NASA GISTEMP, 2024, deg C vs 1951-1980
if anomaly >= 1.0:
    print("at least 1.0 C above the 1951-1980 average")
elif anomaly > 0:
    print("warmer than average")
else:
    print("at or below average")

at least 1.0 C above the 1951-1980 average

Example 2.2 Add the numbers 1 to 100 with a loop.

total = 0
for n in range(1, 101):
    total += n
print(total)

5050

Example 2.3 If CO₂ kept rising 2.5 ppm a year from the 2025 mean, when would it pass 450 ppm?

ppm = 427.35   # 2025 annual mean, Mauna Loa
year = 2025
while ppm < 450:
    ppm += 2.5
    year += 1
print(year, round(ppm, 2))

2035 452.35

Example 3.1 Write and test a Celsius-to-Fahrenheit function.

def c_to_f(c):
    """Convert Celsius to Fahrenheit."""
    return c * 9 / 5 + 32

for c in (-40, 0, 1.28, 100):
    print(c, "->", round(c_to_f(c), 2))

-40 -> -40.0
0 -> 32.0
1.28 -> 34.3
100 -> 212.0

Example 3.2 Turn the stall-speed formula into a function and use it at two densities.

def stall_speed(weight_n, rho, area, cl_max=1.6):
    return (2 * weight_n / (rho * area * cl_max)) ** 0.5

print(round(stall_speed(10791, 1.227, 16.2), 1))
print(round(stall_speed(10791, 0.738, 16.2), 1))

26.0
33.6

Example 3.3 Write a function that returns year-to-year changes of a list (Mauna Loa annual means 2022–2025).

def growth(values):
    """Year-to-year changes of a list."""
    return [b - a for a, b in zip(values, values[1:])]

co2 = [418.53, 421.08, 424.61, 427.35]   # 2022-2025
print([round(d, 2) for d in growth(co2)])

[2.55, 3.53, 2.74]

Example 4.1 Summarize ten days of USGS daily mean flow for the Potomac at Little Falls (June 1–10, 2025, in cubic feet per second).

flows = [44700, 40200, 29100, 22700, 17300, 14100, 12600, 12200, 14300, 22600]
print(len(flows), max(flows), min(flows))
print(sum(flows) / len(flows))
print(sorted(flows)[:3])

10 44700 12200
22980.0
[12200, 12600, 14100]

Example 4.2 Find the warmest year in a dictionary of NASA GISTEMP anomalies.

gistemp = {2016: 1.01, 2020: 1.01, 2023: 1.17, 2024: 1.28, 2025: 1.18}
warmest = max(gistemp, key=gistemp.get)
print(warmest, gistemp[warmest])
print([y for y, a in gistemp.items() if a > 1.1])

2024 1.28
[2023, 2024, 2025]

Example 4.3 Count words with a dictionary.

words = "reduce reuse recycle reduce repair reduce".split()
counts = {}
for w in words:
    counts[w] = counts.get(w, 0) + 1
print(counts)

{'reduce': 3, 'reuse': 1, 'recycle': 1, 'repair': 1}

Example 5.1 Parse the last five rows of NOAA’s Mauna Loa annual-mean CSV and find the average rate of rise.

import csv, io
text = """year,mean,unc
2021,416.41,0.12
2022,418.53,0.12
2023,421.08,0.12
2024,424.61,0.12
2025,427.35,0.12
"""
rows = list(csv.DictReader(io.StringIO(text)))
years = [int(r["year"]) for r in rows]
ppm = [float(r["mean"]) for r in rows]
print(years[0], years[-1], ppm[-1])
print(round((ppm[-1] - ppm[0]) / (years[-1] - years[0]), 2), "ppm per year")

2021 2025 427.35
2.73 ppm per year

Example 5.2 Write a small CSV of USGS flows, read it back and convert to m³/s.

with open("flows.csv", "w") as f:
    f.write("date,cfs\n2025-06-01,44700\n2025-06-02,40200\n2025-06-03,29100\n")
with open("flows.csv") as f:
    next(f)                      # skip the header
    for line in f:
        date, cfs = line.strip().split(",")
        print(date, round(int(cfs) * 0.0283168, 1), "m3/s")

2025-06-01 1265.8 m3/s
2025-06-02 1138.3 m3/s
2025-06-03 824.0 m3/s

Example 5.3 Read lines with comments and bad values without crashing.

lines = ["# comment", "year,mean", "2024,424.61", "2025,bad", "2026,"]
good = []
for line in lines:
    if line.startswith("#") or line.startswith("year"):
        continue
    year, value = line.split(",")
    try:
        good.append((int(year), float(value)))
    except ValueError:
        print("skipping", repr(line))
print(good)

skipping '2025,bad'
skipping '2026,'
[(2024, 424.61)]

Example 6.1 Use NumPy on Mauna Loa annual means 2021–2025: differences, mean and trend.

import numpy as np
years = np.array([2021, 2022, 2023, 2024, 2025])
ppm = np.array([416.41, 418.53, 421.08, 424.61, 427.35])
print(np.diff(ppm))
print(round(ppm.mean(), 3), round(np.diff(ppm).mean(), 3))
slope, intercept = np.polyfit(years, ppm, 1)
print(round(slope, 3))

[2.12 2.55 3.53 2.74]
421.596 2.735
2.796

Example 6.2 Convert ten days of Potomac flow to m³/s, count high-flow days and find the total volume.

import numpy as np
cfs = np.array([44700, 40200, 29100, 22700, 17300, 14100, 12600, 12200, 14300, 22600])
m3s = cfs * 0.0283168
print(np.round(m3s[:3], 1))
print(int((cfs > 20000).sum()), "days above 20,000 cfs")
print(round(m3s.sum() * 86400 / 1e9, 3), "km3 in 10 days")

[1265.8 1138.3  824. ]
5 days above 20,000 cfs
0.562 km3 in 10 days

Example 6.3 Plot selected NASA GISTEMP annual anomalies and save the figure.

import matplotlib
matplotlib.use("Agg")
import matplotlib.pyplot as plt
years = [1880, 1951, 1980, 2000, 2016, 2020, 2023, 2024, 2025]
anom = [-0.17, -0.07, 0.26, 0.39, 1.01, 1.01, 1.17, 1.28, 1.18]
fig, ax = plt.subplots()
ax.plot(years, anom, "o-")
ax.set_xlabel("Year"); ax.set_ylabel("Anomaly (C)")
ax.set_title("NASA GISTEMP v4, selected years")
fig.savefig("gistemp.png")
print(len(ax.lines[0].get_xdata()), "points saved")

9 points saved

Example 7.1 Find the bug in a mean function called on an empty list.

def mean(xs):
    return sum(xs) / len(xs)

try:
    print(mean([]))
except ZeroDivisionError as e:
    print("bug found:", e)

bug found: division by zero

Example 7.2 Fix it with a clear error and add assertion tests.

def mean(xs):
    if not xs:
        raise ValueError("mean() of empty list")
    return sum(xs) / len(xs)

assert mean([1, 2, 3]) == 2
assert abs(mean([0.1, 0.2]) - 0.15) < 1e-12
try:
    mean([])
except ValueError:
    print("all tests passed")

all tests passed

Example 7.3 Write a pytest-style test for f_to_c using known points.

def f_to_c(f):
    return (f - 32) * 5 / 9

def test_f_to_c():
    assert f_to_c(32) == 0
    assert f_to_c(212) == 100
    assert abs(f_to_c(-40) - (-40)) < 1e-12

test_f_to_c()
print("ok")

ok

Example 8.1 Search a sorted list of one million even numbers for the last one, both ways.

def linear_search(xs, target):
    for i, x in enumerate(xs):
        if x == target:
            return i
    return -1

def binary_search(xs, target):
    lo, hi, steps = 0, len(xs) - 1, 0
    while lo <= hi:
        steps += 1
        mid = (lo + hi) // 2
        if xs[mid] == target:
            return mid, steps
        if xs[mid] < target:
            lo = mid + 1
        else:
            hi = mid - 1
    return -1, steps

data = list(range(0, 2_000_000, 2))
print(linear_search(data, 1_999_998))
print(binary_search(data, 1_999_998))

999999
(999999, 20)

Example 8.2 Find √2 by bisection on x² − 2 from [0, 2] to 10⁻⁶.

def bisect_root(f, a, b, tol=1e-6):
    n = 0
    while b - a > tol:
        m = (a + b) / 2
        if f(a) * f(m) <= 0:
            b = m
        else:
            a = m
        n += 1
    return (a + b) / 2, n

root, n = bisect_root(lambda x: x**2 - 2, 0, 2)
print(round(root, 6), n)

1.414214 21

Example 8.3 Compare log₂ n and n² growth.

import math
for n in (10, 1_000, 1_000_000):
    print(n, math.ceil(math.log2(n)), n * n)

10 4 100
1000 10 1000000
1000000 20 1000000000000

Example 9.1 Simulate a 20 m/s projectile at 45° and compare with the exact range.

import math
v0, angle, g, dt = 20.0, math.radians(45), 9.81, 0.001
x, y = 0.0, 0.0
vx, vy = v0 * math.cos(angle), v0 * math.sin(angle)
while y >= 0:
    vy -= g * dt
    x += vx * dt
    y += vy * dt
print(round(x, 2), "m simulated;", round(v0**2 / g, 2), "m exact")

40.77 m simulated; 40.77 m exact

Example 9.2 Simulate one orbit 400 km above Earth with 1-second steps.

import math
GM = 3.986e14
r0 = 6.771e6
x, y = r0, 0.0
vx, vy = 0.0, math.sqrt(GM / r0)
dt, t = 1.0, 0
period = 2 * math.pi * math.sqrt(r0**3 / GM)
while t < period:
    r3 = (x * x + y * y) ** 1.5
    vx -= GM * x / r3 * dt          # Euler-Cromer: update velocity first
    vy -= GM * y / r3 * dt
    x += vx * dt
    y += vy * dt
    t += dt
print(round(period / 60, 1), "min;", round(math.hypot(x, y) / 1000), "km from Earth's center")

92.4 min; 6771 km from Earth's center

Example 9.3 Step a logistic fish population (r = 0.5, K = 1,000, P₀ = 100) for four years.

r, K, P, h = 0.5, 1000, 100.0, 1.0
for year in range(1, 5):
    P += h * r * P * (1 - P / K)
    print(year, round(P, 1))

1 145.0
2 207.0
3 289.1
4 391.8