Use NumPy arrays, plots, fitting, numerical methods, units and instrument data files to do real engineering and science calculations in Python.

Scientific Python

## Why this matters

A plot is often the only part of your work a reader looks at closely. A good scientific plot answers three questions at a glance: **what was measured** (axis labels with units), **what the values are** (points, not just a smooth line), and **how sure you are** (error bars). Reviewers, supervisors and clients will ask about all three. matplotlib is the standard Python plotting library, and its object-oriented interface gives you full control.

## From replicates to error bars

Most plotted points are a summary of several replicate measurements. Suppose you measure the extension of a spring at four loads, three times each:

```python
import numpy as np

load_N = np.array([1.0, 2.0, 3.0, 4.0])
extension_mm = np.array([
    [2.1, 2.0, 2.2],
    [4.0, 4.3, 4.1],
    [6.2, 5.9, 6.1],
    [8.3, 8.0, 8.2],
])                                         # rows = loads, columns = replicates

mean = extension_mm.mean(axis=1)           # one mean per load
sd = extension_mm.std(axis=1, ddof=1)      # sample SD per load
n = extension_mm.shape[1]
sem = sd / np.sqrt(n)                      # standard error of the mean
```

Decide, and **say in the caption**, what your error bars show:

- **SD**: the spread of individual measurements.
- **SEM** (standard error of the mean): how well you know the mean. It shrinks with more replicates.
- **95% confidence interval**: roughly `2 × SEM` for large `n`; for small `n` multiply SEM by the t-value (4.30 for n = 3).

Unlabelled error bars are a common reason papers get sent back.

## The matplotlib object interface

Use `fig, ax = plt.subplots()` and call methods on `ax`. This scales to multi-panel figures and avoids the confusing global state of `plt.plot`.

```python
import matplotlib.pyplot as plt

fig, ax = plt.subplots(figsize=(4.5, 3.2), dpi=150)

ax.errorbar(load_N, mean, yerr=sem, fmt="o", capsize=3,
            color="black", label="measured (mean ± SEM, n = 3)")

# Fitted line for comparison
k, c = np.polyfit(load_N, mean, 1)
x_fit = np.linspace(0, 4.5, 100)
ax.plot(x_fit, k * x_fit + c, "-", color="tab:blue",
        label=f"fit: {k:.2f} mm/N")

ax.set_xlabel("Load (N)")
ax.set_ylabel("Extension (mm)")
ax.set_xlim(0, 4.5)
ax.set_ylim(bottom=0)
ax.legend(frameon=False)
fig.tight_layout()
fig.savefig("spring.png")      # or .pdf / .svg for vector output
```

Matplotlib is not available in this course's in-browser exercises, so run plotting code in Jupyter, Google Colab or a local Python install. The calculations that feed the plot (means, SEMs, fits) are plain NumPy and work anywhere.

## What makes a plot "publication quality"

1. **Axis labels with units** in a consistent style, such as `Load (N)` or `Temperature / K`. Never leave the default blank axes.
2. **Markers for data, lines for models.** Measured points are discrete; connecting them with a line suggests you measured in between.
3. **Readable sizes.** Make the figure at its final printed size (a journal column is about 8.5 cm, or 3.3 inches, wide) so fonts are not shrunk to illegibility. Font size 8 to 10 pt at final size is typical.
4. **Vector formats** (`.pdf`, `.svg`) for line art; `.png` at 300 dpi or more for raster.
5. **Colour that survives.** Use colour-blind-safe palettes (matplotlib's default `tab10` and the `viridis` colormap are good) and vary marker shape as well as colour, so it still works in greyscale.
6. **Honest axes.** Start bar charts at zero. If you zoom a line plot, make it obvious.
7. **A caption** that says what the error bars are and what `n` is.

## Several panels

```python
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(7, 3), sharey=False)
ax1.errorbar(load_N, mean, yerr=sem, fmt="o", capsize=3)
residuals = mean - (k * load_N + c)
ax2.axhline(0, color="grey", lw=0.8)
ax2.plot(load_N, residuals, "s")
ax2.set_xlabel("Load (N)")
ax2.set_ylabel("Residual (mm)")
```

A **residual plot** next to a fit is one of the most useful panels you can add: if the residuals show a pattern (a curve, a funnel), the model is missing something.

## Common mistakes

- **Using `std()` without `ddof=1`** for three replicates: it underestimates the spread by about 18%.
- **Mixing up SD and SEM** or not saying which you plotted.
- **Plotting means with no indication of `n`.** Three replicates and three hundred look identical without it.
- **Default tiny figures saved at screen resolution**, then stretched in a report.
- **Rainbow (`jet`) colormaps**, which create false boundaries and fail for colour-blind readers.
- **Units in the title instead of on the axes.**

## Quick recap

- Summarise replicates with `mean(axis=...)`, `std(ddof=1)` and `SEM = SD / sqrt(n)`.
- `fig, ax = plt.subplots()`, then `ax.errorbar(x, y, yerr=..., fmt="o", capsize=3)`.
- Label every axis with units, use markers for data and lines for models, and state what the error bars mean.
- Save at final size in a vector format.

## Try it

In the **Do** tab you will compute the means and standard errors behind an error-bar plot and write a helper that formats axis labels with units.


Turn replicate measurements into means and standard errors and plot them with matplotlib, error bars, labelled axes and units.

Publication-quality plots

Matplotlib Tutorial: For Physicists, Engineers, and Mathematicians

Standard Deviation vs Standard Error, Clearly Explained!!!

## Why this matters

Many physical laws are not straight lines. A pendulum's period grows with the square root of its length. Radioactive activity, capacitor voltage during discharge and first-order reaction concentrations fall exponentially. Fatigue life, sound intensity and earthquake energy span many powers of ten. Plotted on ordinary axes, these curves hide their structure and squash small values into a corner. **Logarithmic axes** spread them out and, for the right kind of law, turn the curve into a straight line whose slope you can read off and fit.

## Two shapes, two kinds of axis

**Power law**: `y = a x^b`. Take logs of both sides:

```text
log y = log a + b log x
```

On **log-log** axes (both axes logarithmic) this is a straight line with slope `b` (the exponent) and intercept `log a`.

**Exponential**: `y = A e^(k t)` (decay when `k < 0`). Take the natural log of `y` only:

```text
ln y = ln A + k t
```

On **semilog** axes (log `y`, linear `t`) this is a straight line with slope `k`.

A quick diagnostic: plot your data both ways. Whichever plot is straight tells you the type of law.

## Plotting with log axes

```python
import numpy as np
import matplotlib.pyplot as plt

L = np.array([0.10, 0.20, 0.40, 0.80, 1.60])      # pendulum length, m
T = np.array([0.636, 0.898, 1.268, 1.795, 2.538])  # period, s

fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(7, 3))
ax1.plot(L, T, "o")
ax1.set_xlabel("Length (m)"); ax1.set_ylabel("Period (s)")

ax2.loglog(L, T, "o")          # same as ax.set_xscale("log"); ax.set_yscale("log")
ax2.set_xlabel("Length (m)"); ax2.set_ylabel("Period (s)")
```

`ax.semilogy(t, y)` sets only the y-axis to log. A log axis keeps the **original values** on its tick labels (0.1, 1, 10), which is easier to read than plotting `np.log10(y)` on a linear axis.

## Fitting on the linearised data

To get numbers, fit a straight line to the transformed data with `np.polyfit(x, y, 1)`, which returns `[slope, intercept]`.

### Worked example 1: pendulum power law

```python
b, log_a = np.polyfit(np.log(L), np.log(T), 1)
a = np.exp(log_a)
print(round(b, 3), round(a, 3))     # 0.5 2.006
```

The exponent is 0.5, matching `T = 2π sqrt(L / g)`. The prefactor `a = 2π / sqrt(g)` gives `g = (2π / a)^2 ≈ 9.81 m/s²`.

Which log base you use does not matter for the exponent, as long as you use the same base on both axes; the prefactor must be converted back with the matching inverse (`np.exp` for `np.log`, `10**` for `np.log10`).

### Worked example 2: capacitor discharge

A capacitor discharges through a resistor: `V(t) = V0 e^(-t / τ)`.

```python
t = np.array([0, 1, 2, 3, 4, 5], dtype=float)       # ms
V = np.array([9.00, 5.46, 3.31, 2.01, 1.22, 0.74])  # volts

k, ln_V0 = np.polyfit(t, np.log(V), 1)
tau = -1 / k
print(round(tau, 2), round(np.exp(ln_V0), 2))       # 2.0 ms, 9.0 V
```

The half-life is `τ ln 2 ≈ 1.39 ms`: the voltage halves every 1.39 ms, whatever its starting value. That constant-ratio behaviour is the signature of an exponential.

## Caveats of fitting in log space

Linearised fits are quick, robust and usually good enough. Know their limits:

- **Zero and negative values** have no logarithm. Remove them (and say so) or use a model that handles an offset, such as `y = A e^(kt) + C`, fitted with non-linear least squares (next module).
- **Noise is reshaped.** A least-squares fit to `log y` weights small values more heavily, because a fixed error is a bigger fraction of a small number. If your noise is constant in absolute terms, the log fit gives slightly biased parameters; a non-linear fit on the original scale is more correct.
- **Offsets break straightness.** A background level (dark current, baseline absorbance) makes a semilog plot curve at long times. Subtract the background first if you know it.
- **Range matters.** A power law that is straight over two decades may bend outside them. Do not extrapolate beyond the data.

## Common mistakes

- **Fitting `y` against `log x`** when you meant log-log, or forgetting to transform `y` for an exponential.
- **Reporting `log a` as `a`.** Convert the intercept back with `exp` or `10**`.
- **Mixing bases**: fitting with `np.log10` and converting back with `np.exp`.
- **Taking the log of data containing zeros**, which gives `-inf` and a `RuntimeWarning`, then a fit full of `nan`.
- **Reading slopes off a log axis by eye using tick distances** instead of fitting the transformed values.

## Quick recap

- Power law `y = a x^b`: straight on log-log; slope = `b`.
- Exponential `y = A e^(kt)`: straight on semilog-y; slope = `k`, time constant `τ = -1/k`.
- `ax.loglog`, `ax.semilogy` for plotting; `np.polyfit` on the logged data for numbers.
- Watch for zeros, offsets and noise reshaping; use a non-linear fit when they matter.

## Try it

In the **Do** tab you will write power-law and exponential fitting functions and test them on noisy, seeded data.


Use log-log and semilog axes to turn power laws and exponentials into straight lines, then read the exponent or rate from a linear fit.

Log scales and fits

Using a log-log graph to find a power law | A Level Physics

How to Linearize Data in Physics Lab

## Why this matters

A lot of scientific data is two-dimensional: a temperature map across a heated plate, a stress field in a beam from finite-element software, a detector image from a microscope or X-ray camera, a potential-energy surface, a gel scan. In NumPy these are just 2D arrays. The skills are the same each time: build the coordinate grid, compute or load the values, find the interesting features (peaks, hot spots, regions above a threshold), and show the result as a **heatmap** or **contour plot** that a reader can interpret.

## Core idea: a grid of coordinates

`np.meshgrid` turns two 1D coordinate arrays into two 2D arrays, one holding the x-coordinate of every grid point and one holding the y-coordinate:

```python
import numpy as np

x = np.linspace(-2, 2, 201)        # 201 columns
y = np.linspace(-1, 1, 101)        # 101 rows
X, Y = np.meshgrid(x, y)           # both have shape (101, 201)
```

Note the shape: **rows follow y, columns follow x**, so `Z[i, j]` is the value at `(x[j], y[i])`. This "row = y" convention matches how images are stored and is the single biggest source of confusion with 2D data.

With the grid in hand, any formula becomes a 2D field by broadcasting:

```python
# A Gaussian hot spot centred at (0.5, -0.2) on a 20 °C plate
Z = 20 + 60 * np.exp(-((X - 0.5)**2 + (Y + 0.2)**2) / (2 * 0.3**2))
```

## Finding features

`np.argmax` returns the position of the largest value in the **flattened** array. Convert it back to a row and column with `np.unravel_index`:

```python
i, j = np.unravel_index(np.argmax(Z), Z.shape)
print(x[j], y[i], Z[i, j])          # 0.5 -0.2 80.0
```

Thresholds work with Boolean masks, exactly as in 1D:

```python
hot = Z > 50                         # 2D Boolean array
print(hot.mean())                    # fraction of the plate above 50 °C
cell_area = (x[1] - x[0]) * (y[1] - y[0])
print(hot.sum() * cell_area)         # hot area in square units
```

Reductions take an axis. `Z.mean(axis=0)` averages over rows (giving a profile along x); `Z.max(axis=1)` gives the hottest value in each row. A **line profile** through the peak is just a slice: `Z[i, :]`.

## Heatmaps with imshow and pcolormesh

```python
import matplotlib.pyplot as plt

fig, ax = plt.subplots(figsize=(5, 3))
im = ax.imshow(Z, origin="lower", extent=[x[0], x[-1], y[0], y[-1]],
               cmap="viridis", aspect="auto")
fig.colorbar(im, ax=ax, label="Temperature (°C)")
ax.set_xlabel("x (m)")
ax.set_ylabel("y (m)")
```

Three arguments matter:

- `origin="lower"` puts row 0 at the bottom, as on a graph. The default (`"upper"`) is image convention, with row 0 at the top. Getting this wrong flips your plot vertically.
- `extent` sets the axis values, so ticks show metres rather than pixel indices.
- A **colorbar with a label and unit** is required; a heatmap without one is unreadable.

`ax.pcolormesh(X, Y, Z, shading="auto")` does the same job and also handles grids that are not evenly spaced.

## Contours

Contour lines join points of equal value, like a topographic map:

```python
cs = ax.contour(X, Y, Z, levels=[30, 40, 50, 60, 70], colors="white", linewidths=0.8)
ax.clabel(cs, fmt="%d °C", fontsize=7)
```

`contourf` fills between levels. Contours are better than colour when readers need to read actual values, and they print well in greyscale. Overlaying a few contour lines on a heatmap often gives the best of both.

## Colormaps

- **Sequential** maps (`viridis`, `magma`, `cividis`) for values that go from low to high, like temperature or intensity.
- **Diverging** maps (`RdBu_r`, `coolwarm`) for values with a meaningful centre, like residuals or tension versus compression. Centre them on zero with `vmin=-m, vmax=m`.
- Avoid `jet` and `rainbow`: they create bands that look like features and are hard for colour-blind readers.

Fix `vmin` and `vmax` when comparing several images side by side, or each panel gets its own scale and the comparison is meaningless.

## Common mistakes

- **Swapping x and y.** `Z[i, j]` is `(y[i], x[j])`. Check with an asymmetric test pattern.
- **Forgetting `origin="lower"`**, so the plot is upside down compared with your coordinates.
- **No `extent`**, so the axes show pixel numbers instead of physical units.
- **Comparing panels with different colour scales.**
- **Using `argmax` and treating the result as a column index**: it is a flat index until you unravel it.

## Quick recap

- `X, Y = np.meshgrid(x, y)` gives 2D coordinate arrays of shape `(len(y), len(x))`.
- Peaks: `np.unravel_index(np.argmax(Z), Z.shape)`. Thresholds: Boolean masks, `.mean()` for a fraction, `.sum() * cell_area` for an area.
- `imshow(..., origin="lower", extent=...)` plus a labelled colorbar; `contour` for lines of equal value.
- Sequential colormaps for magnitudes, diverging ones for signed data.

## Try it

In the **Do** tab you will build a temperature field on a grid, find its hottest point, and measure the area above a threshold.


Build 2D grids with meshgrid, find peaks and thresholds in 2D arrays, and show them as heatmaps and contour plots.

Images and 2D data

Heatmap Visualization with Matplotlib in Python | Complete Guide

Publication-quality matplotlib plots with error bars and units, log and semilog axes to linearise power laws and exponentials, and heatmaps and contours of 2D data.