Fyrefly
v0.12.0

Fyrefly

A quant-based Python library for flexible math, financial formulas, equation solving, SQL-powered data analysis, AI-driven insights, visualization, and API testing — every function takes any number of arguments.

Installation

Fyrefly is a single package — every function in this reference comes from one install.

$pip install fyrefly

Quickstart

Most real usage touches several parts of the library in one script. Here's a realistic end-to-end example — load a file, clean it, query it, chart it, and ask a question about it in plain English.

from fyrefly import load, clean, sql, viz, ask

# Load a dataset — CSV, Excel, or a Google Sheet all work
c = load("sales.csv")

# Clean it up: dedupe, strip whitespace, standardize column names
c = clean(c)

# Query it with real SQL, by the variable name you gave it
top_regions = sql("select region, sum(amount) as total from c group by region order by total desc")

# Visualize the result — one line, no separate charting library
viz.bar("region", "amount")

# Or just ask in plain English (needs your own API key)
ask(c, "which region grew the fastest last quarter?", api_key="AIza...")

Everything below is organized by what it does. Math, sequences, and geometry need nothing but Python. Data operations need a file to load. The AI functions are the only ones that need an API key of your own — everything else runs completely offline.

Math

The building blocks. Each function takes any number of arguments and applies the operation left to right — add(1, 2, 3, 4) sums all four; div(100, 5, 2) divides 100 by 5, then that result by 2.

add(*nos)

Adds any count of numbers.

add(5, 4, 6, 10)
# 25

add(1, 2, 3, 4, 5)
# 15

mul(*nos)

Multiplies any count of numbers.

mul(5, 4, 6, 10)
# 1200

div(*nos)

Divides numbers in sequence: first ÷ second ÷ third ...

div(100, 5, 2)
# 10.0

mod(*nos)

Finds the remainder in sequence: first % second % third ...

mod(17, 5)
# 2

diff(*nos)

Subtracts numbers in sequence: first − second − third ...

diff(10, 3, 2)
# 5

Quant & Finance

Interest and profit/loss formulas for everyday financial math. profit() and loss() are deliberately strict: call the wrong one — say profit() on numbers that are actually a loss — and you get a clear error telling you which function to use instead, rather than a silently wrong (negative) answer.

si(p, n, r)

Simple Interest: (p × n × r) / 100.

Formulas and arguments
CallFormulaWhat each argument means
si(p, n, r)
  • p () — the principal (money you start with)
  • n — time, in years
  • r — rate of interest, % per year
Returns the simple interest earned
si(1000, 2, 5)
# 100.0

ci(p, n, r, freq="annual")

Compound interest earned (principal excluded). freq is "annual" or "semi-annual".

Formulas and arguments
CallFormulaWhat each argument means
ci(p, n, r)
  • p () — the principal
  • n — time, in years
  • r — rate, % per year
Returns the compound interest earned (the principal is not included)
ci(p, n, r, "semi-annual")
  • freq — "annual" (default) or "semi-annual": half the rate, twice as many periods
Returns the compound interest earned
ci(1000, 2, 10)
# 210.0

ci(1000, 2, 10, "semi-annual")
# 215.5

si_ci_diff(p, r, n)

Difference between Compound and Simple Interest. Supports n = 2 or 3 years only.

Formulas and arguments
CallFormulaWhat each argument means
si_ci_diff(p, r, n=2)
  • p () — the principal
  • r — rate, % per year
  • n — years, 2 or 3 only
Returns the difference between compound and simple interest
si_ci_diff(p, r, n=3)
  • n — 3 years
Returns
si_ci_diff(1000, 10, 2)
# 10.0

si_ci_diff(1000, 10, 3)
# 31.0

profit(sp, cp, pct=None)

Profit = sp − cp. Pass "%" as the third argument for (value, percent). Raises a clear error if sp < cp — that's a loss, not a profit.

Formulas and arguments
CallFormulaWhat each argument means
profit(sp, cp)
  • sp () — selling price
  • cp () — cost price
Returns the profit
profit(sp, cp, "%")
  • pct — pass "%" as the third argument
Returns the pair (profit, profit percentage)
profit(120, 100)
# 20

profit(120, 100, "%")
# (20, 20.0)

profit(80, 100)
# raises: this is a LOSS. Use loss(80, 100) instead.

loss(sp, cp, pct=None)

Loss = cp − sp. Pass "%" as the third argument for (value, percent). Raises a clear error if sp > cp — that's a profit, not a loss.

Formulas and arguments
CallFormulaWhat each argument means
loss(sp, cp)
  • sp () — selling price
  • cp () — cost price
Returns the loss
loss(sp, cp, "%")
  • pct — pass "%" as the third argument
Returns the pair (loss, loss percentage)
loss(80, 100)
# 20

loss(80, 100, "%")
# (20, 20.0)

Sequences

Each progression type works the same way: call it directly for the full list of terms, or use .tn() and .sn() when you only need one term or a running total — no need to generate the whole sequence first.

ap(a, d, n) · ap.tn(a, d, n) · ap.sn(a, d, n)

Arithmetic Progression: full sequence, the nth term, or the sum of the first n terms.

Formulas and arguments
CallFormulaWhat each argument means
ap.tn(a, d, n)
  • a — the first term
  • d — the common difference
  • n — which term you want
Returns the th term
ap.sn(a, d, n)
  • a — the first term
  • d — the common difference
  • n — how many terms to add
Returns the sum of the first terms
ap(2, 3, 4)
# [2, 5, 8, 11]

ap.tn(2, 3, 4)
# 11

ap.sn(2, 3, 4)
# 26.0

gp(a, r, n) · gp.tn(a, r, n) · gp.sn(a, r, n)

Geometric Progression: full sequence, the nth term, or the sum of the first n terms.

Formulas and arguments
CallFormulaWhat each argument means
gp.tn(a, r, n)
  • a — the first term
  • r — the common ratio
  • n — which term you want
Returns the th term
gp.sn(a, r, n)
  • a — the first term
  • r — the common ratio (if the sum is )
  • n — how many terms to add
Returns the sum of the first terms
gp(2, 3, 4)
# [2, 6, 18, 54]

gp.tn(2, 3, 4)
# 54

gp.sn(2, 3, 4)
# 80.0

hp(a, d, n) · hp.tn(a, d, n) · hp.sn(a, d, n)

Harmonic Progression — reciprocals of an AP with first term 1/a. No closed-form sum exists, so .sn() totals the actual terms.

Formulas and arguments
CallFormulaWhat each argument means
hp.tn(a, d, n)
  • a — the first term of the progression
  • d — the common difference of the reciprocals
  • n — which term you want
Returns the th term
hp.sn(a, d, n)
  • a — the first term
  • d — the common difference of the reciprocals
  • n — how many terms to add
Returns the sum (there is no simple closed form, so the terms are added)
hp(2, 3, 4)
# [2.0, 0.2857..., 0.1538..., 0.1052...]

hp.tn(2, 3, 4)
# 0.10526315789473684

Averages & Statistics

avg() takes raw numbers directly; avgc() is for when you already know a value's frequency rather than listing it out repeatedly — e.g. avgc(5, 3) for "the value 5, occurring 3 times", instead of avg(5, 5, 5).

avg(*nos)

Plain average of individual numbers.

Formulas and arguments
CallFormulaWhat each argument means
avg(*nos)
  • nos () — the numbers, separated by commas
Returns the mean
avg(10, 12, 13)
# 11.666666666666666

avgc(*value_count_pairs)

Weighted average from (value, count) pairs — e.g. "this value occurred this many times."

Formulas and arguments
CallFormulaWhat each argument means
avgc(*pairs)
  • pairs () — (value, count) pairs: each value and how many times it occurred
Returns the weighted mean
avgc(5, 3)
# 5.0   (three 5's)

avgc(5, 3, 10, 2)
# 7.0   (three 5's and two 10's)

mode(*nos)

Most frequently occurring value(s). Returns a list — may contain more than one value if there's a tie.

Formulas and arguments
CallFormulaWhat each argument means
mode(*nos)
  • nos () — the numbers, separated by commas
Returns a list of the most frequent value(s)
mode(1, 2, 2, 3)
# [2]

mode(1, 1, 2, 2)
# [1, 2]

median(*nos)

Median of any count of numbers.

Formulas and arguments
CallFormulaWhat each argument means
median(*nos)
  • nos () — the numbers, separated by commas
Returns the median
median(1, 3, 2)
# 2

median(1, 2, 3, 4)
# 2.5

Geometry

Coordinate geometry formulas you'd otherwise look up every time. Points are plain (x, y) tuples throughout, so these compose easily with your own data.

hyp(a, b)

Hypotenuse of a right triangle given the other two sides.

Formulas and arguments
CallFormulaWhat each argument means
hyp(a, b)
  • a — one of the two shorter sides
  • b — the other shorter side
Returns the hypotenuse (Pythagoras' theorem)
hyp(3, 4)
# 5.0

slope(point1, point2)

Slope between two (x, y) tuples: (y2 − y1) / (x2 − x1). Raises a clear error for a vertical line.

Formulas and arguments
CallFormulaWhat each argument means
slope(point1, point2)
  • point1 () — the first point, as an (x, y) tuple
  • point2 () — the second point
Returns the slope (an error for a vertical line)
slope((1, 2), (3, 4))
# 1.0

centroid(point1, point2, point3)

Centroid of a triangle given three (x, y) tuples.

Formulas and arguments
CallFormulaWhat each argument means
centroid(point1, point2, point3)
  • point1 () — first corner as an (x, y) tuple
  • point2 () — second corner
  • point3 () — third corner
Returns the centroid as an (x, y) tuple
centroid((1, 2), (3, 2), (4, 5))
# (2.6666666666666665, 3.0)

dist(point1, point2)

Distance between two (x, y) tuples: √((x2−x1)² + (y2−y1)²).

Formulas and arguments
CallFormulaWhat each argument means
dist(point1, point2)
  • point1 () — the first point, as an (x, y) tuple
  • point2 () — the second point
Returns the distance
dist((6, 2), (10, 5))
# 5.0

Mensuration — 2D Shapes

Area and perimeter for every common 2D shape. Each shape is its own object with .area() and .perimeter() (or close variants) as methods.

circle.area(r) · circle.perimeter(r)

Area and circumference of a circle.

Formulas and arguments
CallFormulaWhat each argument means
circle.area(r)
  • r — radius
Returns the area
circle.perimeter(r)
  • r — radius
Returns the circumference
circle.area(7)
# 153.93804002589985

circle.perimeter(7)
# 43.982297150257104

semicircle.area(r) · semicircle.perimeter(r)

Area and perimeter (curved edge + diameter) of a semicircle.

Formulas and arguments
CallFormulaWhat each argument means
semicircle.area(r)
  • r — radius
Returns the area
semicircle.perimeter(r)
  • r — radius
Returns the curved edge plus the diameter
semicircle.area(7)
# 76.96902001294993

semicircle.perimeter(7)
# 35.99114857512855

square.area(s) · square.perimeter(s)

Area and perimeter of a square.

Formulas and arguments
CallFormulaWhat each argument means
square.area(s)
  • s — side
Returns the area
square.perimeter(s)
  • s — side
Returns the perimeter
square.area(5)
# 25

square.perimeter(5)
# 20

rectangle.area(l, w) · rectangle.perimeter(l, w)

Area and perimeter of a rectangle.

Formulas and arguments
CallFormulaWhat each argument means
rectangle.area(l, w)
  • l — length
  • w — width
Returns the area
rectangle.perimeter(l, w)
  • l — length
  • w — width
Returns the perimeter
rectangle.area(4, 6)
# 24

rectangle.perimeter(4, 6)
# 20

triangle.area(base, height) · triangle.area_sides(a, b, c) · triangle.perimeter(a, b, c)

Area from base/height, area from all three sides via Heron's formula, and perimeter. Raises a clear error if the three sides can't form a real triangle.

Formulas and arguments
CallFormulaWhat each argument means
triangle.area(base, height)
  • base () — the base
  • height () — the height, measured at right angles to the base
Returns the area
triangle.area_sides(a, b, c)
  • a — first side
  • b — second side
  • c — third side
Returns the area (Heron's formula); an error if they cannot form a triangle
triangle.perimeter(a, b, c)
  • a — first side
  • b — second side
  • c — third side
Returns the perimeter
triangle.area(10, 5)
# 25.0

triangle.area_sides(3, 4, 5)
# 6.0

triangle.perimeter(3, 4, 5)
# 12

equilateral_triangle.area(s) · equilateral_triangle.perimeter(s)

Area and perimeter of an equilateral triangle from its side length.

Formulas and arguments
CallFormulaWhat each argument means
equilateral_triangle.area(s)
  • s — side
Returns the area
equilateral_triangle.perimeter(s)
  • s — side
Returns the perimeter
equilateral_triangle.area(4)
# 6.928203230275509

equilateral_triangle.perimeter(4)
# 12

parallelogram.area(base, height) · parallelogram.perimeter(a, b)

Area and perimeter of a parallelogram.

Formulas and arguments
CallFormulaWhat each argument means
parallelogram.area(base, height)
  • base () — the base
  • height () — the perpendicular height (not the slanted side)
Returns the area
parallelogram.perimeter(a, b)
  • a — one side
  • b — the adjacent side
Returns the perimeter
parallelogram.area(10, 5)
# 50

parallelogram.perimeter(10, 5)
# 30

rhombus.area(d1, d2) · rhombus.perimeter(s)

Area from the two diagonals, and perimeter from the side length.

Formulas and arguments
CallFormulaWhat each argument means
rhombus.area(d1, d2)
  • d1 () — the first diagonal
  • d2 () — the second diagonal
Returns the area
rhombus.perimeter(s)
  • s — side
Returns the perimeter
rhombus.area(6, 8)
# 24.0

rhombus.perimeter(5)
# 20

kite.area(d1, d2) · kite.perimeter(a, b)

Area from the two diagonals, and perimeter from the two distinct side lengths.

Formulas and arguments
CallFormulaWhat each argument means
kite.area(d1, d2)
  • d1 () — the first diagonal
  • d2 () — the second diagonal
Returns the area
kite.perimeter(a, b)
  • a — one of the two side lengths
  • b — the other
Returns the perimeter
kite.area(6, 8)
# 24.0

kite.perimeter(5, 7)
# 24

trapezoid.area(a, b, height) · trapezoid.perimeter(a, b, c, d)

Area from the two parallel sides and height, and perimeter from all four sides.

Formulas and arguments
CallFormulaWhat each argument means
trapezoid.area(a, b, height)
  • a — one parallel side
  • b — the other parallel side
  • height () — the distance between them
Returns the area
trapezoid.perimeter(a, b, c, d)
  • a, b, c, d — the four sides
Returns the perimeter
trapezoid.area(8, 5, 4)
# 26.0

trapezoid.perimeter(8, 5, 3, 4)
# 20

ellipse.area(a, b) · ellipse.perimeter(a, b)

Area is exact. Perimeter uses Ramanujan's approximation — extremely accurate, since no exact closed-form formula exists for an ellipse's perimeter.

Formulas and arguments
CallFormulaWhat each argument means
ellipse.area(a, b)
  • a — the semi-major axis (half the long width)
  • b — the semi-minor axis (half the short width)
Returns the area
ellipse.perimeter(a, b)
  • a — the semi-major axis
  • b — the semi-minor axis
Returns the perimeter (Ramanujan's approximation)
ellipse.area(5, 3)
# 47.12388980384689

ellipse.perimeter(5, 3)
# ~25.5 (approximate)

polygon.area(n_sides, length) · polygon.perimeter(n_sides, length)

Area and perimeter of any regular polygon, given its number of sides and side length.

Formulas and arguments
CallFormulaWhat each argument means
polygon.area(n_sides, length)
  • n_sides () — the number of sides
  • length () — the length of each side
Returns the area of a regular polygon
polygon.perimeter(n_sides, length)
  • n_sides () — the number of sides
  • length () — the length of each side
Returns the perimeter
polygon.area(6, 4)
# area of a regular hexagon

polygon.perimeter(6, 4)
# 24

sector.area(r, angle) · sector.arc_length(r, angle)

Area and arc length of a circular sector. Angle is in degrees.

Formulas and arguments
CallFormulaWhat each argument means
sector.area(r, angle)
  • r — radius
  • angle () — the angle at the centre, in degrees
Returns the area of the sector
sector.arc_length(r, angle)
  • r — radius
  • angle () — the angle at the centre, in degrees
Returns the length of the arc
sector.area(7, 90)
# 38.48451000647496

sector.arc_length(7, 90)
# 10.995574287564276

annulus.area(outer_r, inner_r)

Area of a ring shape — the area between two concentric circles.

Formulas and arguments
CallFormulaWhat each argument means
annulus.area(outer_r, inner_r)
  • outer_r () — the outer radius
  • inner_r () — the inner radius
Returns the area of the ring
annulus.area(10, 6)
# 201.06192982974676

Mensuration — 3D Shapes

Volume and surface area for every common 3D shape. Surface area is split into .tsa() (Total Surface Area, every face combined), .csa() (Curved Surface Area, curved part only — for shapes with a genuine curve), and .lsa() (Lateral Surface Area, side faces only — for flat-sided shapes). Sphere and Torus use plain .surface_area() instead, since there's nothing to disambiguate on those two.

cube.volume(s) · cube.lsa(s) · cube.tsa(s)

Volume, lateral surface area (4 side faces), and total surface area (all 6 faces) of a cube.

Formulas and arguments
CallFormulaWhat each argument means
cube.volume(s)
  • s — edge length
Returns the volume
cube.lsa(s)
  • s — edge length
Returns the four side faces
cube.tsa(s)
  • s — edge length
Returns all six faces
cube.volume(3)
# 27

cube.lsa(3)
# 36

cube.tsa(3)
# 54

cuboid.volume(l, w, h) · cuboid.lsa(l, w, h) · cuboid.tsa(l, w, h)

Volume, lateral surface area, and total surface area of a cuboid.

Formulas and arguments
CallFormulaWhat each argument means
cuboid.volume(l, w, h)
  • l — length
  • w — width
  • h — height
Returns the volume
cuboid.lsa(l, w, h)
  • l — length
  • w — width
  • h — height
Returns the four side faces
cuboid.tsa(l, w, h)
  • l — length
  • w — width
  • h — height
Returns all six faces
cuboid.volume(2, 3, 4)
# 24

cuboid.lsa(2, 3, 4)
# 40

cuboid.tsa(2, 3, 4)
# 52

sphere.volume(r) · sphere.surface_area(r)

Volume and surface area of a sphere. No TSA/CSA split — the whole surface is uniform.

Formulas and arguments
CallFormulaWhat each argument means
sphere.volume(r)
  • r — radius
Returns the volume
sphere.surface_area(r)
  • r — radius
Returns the surface area
sphere.volume(3)
# 113.09733552923254

sphere.surface_area(3)
# 113.09733552923255

hemisphere.volume(r) · hemisphere.csa(r) · hemisphere.tsa(r)

Volume, curved surface area (dome only), and total surface area (dome + flat base) of a hemisphere.

Formulas and arguments
CallFormulaWhat each argument means
hemisphere.volume(r)
  • r — radius
Returns the volume
hemisphere.csa(r)
  • r — radius
Returns the curved dome only
hemisphere.tsa(r)
  • r — radius
Returns the dome plus the flat base ()
hemisphere.volume(3)
# 56.54866776461627

hemisphere.csa(3)
# 56.548667764616276

hemisphere.tsa(3)
# 84.82300164692441

cylinder.volume(r, h) · cylinder.csa(r, h) · cylinder.tsa(r, h)

Volume, curved surface area (side only), and total surface area (side + both circular ends) of a cylinder.

Formulas and arguments
CallFormulaWhat each argument means
cylinder.volume(r, h)
  • r — radius of the base
  • h — height
Returns the volume
cylinder.csa(r, h)
  • r — radius of the base
  • h — height
Returns the side only
cylinder.tsa(r, h)
  • r — radius of the base
  • h — height
Returns the side plus both circular ends
cylinder.volume(3, 7)
# 197.92033717615698

cylinder.csa(3, 7)
# 131.94689145077132

cylinder.tsa(3, 7)
# 188.49555921538757

cone.volume(r, h) · cone.csa(r, h) · cone.tsa(r, h)

Volume, curved surface area, and total surface area of a cone. Slant height is computed internally from r and h.

Formulas and arguments
CallFormulaWhat each argument means
cone.volume(r, h)
  • r — radius of the base
  • h — vertical height
Returns the volume
cone.csa(r, h)
  • r — radius of the base
  • h — vertical height (, the slant height, is worked out for you)
Returns the curved surface
cone.tsa(r, h)
  • r — radius of the base
  • h — vertical height
Returns the curved surface plus the base
cone.volume(3, 4)
# 37.69911184307752

cone.csa(3, 4)
# 47.12388980384689 (slant height = 5)

cone.tsa(3, 4)
# 75.39822368615503

frustum.volume(r1, r2, h) · frustum.csa(r1, r2, h) · frustum.tsa(r1, r2, h)

Volume, curved surface area, and total surface area of a frustum (a cone with the tip cut off). Slant height is computed internally.

Formulas and arguments
CallFormulaWhat each argument means
frustum.volume(r1, r2, h)
  • r1 () — radius of one end
  • r2 () — radius of the other end
  • h — vertical height between them
Returns the volume
frustum.csa(r1, r2, h)
  • r1 () — radius of one end
  • r2 () — radius of the other end
  • h — vertical height
Returns the curved surface
frustum.tsa(r1, r2, h)
  • r1 () — radius of one end
  • r2 () — radius of the other end
  • h — vertical height
Returns the curved surface plus both ends
frustum.volume(3, 5, 6)
# 307.8760800517997

frustum.csa(3, 5, 6)
# 158.95341225273762

frustum.tsa(3, 5, 6)
# 265.76756247479057

pyramid.volume(s, h) · pyramid.lsa(s, h) · pyramid.tsa(s, h)

Square-based pyramid. h is the vertical height — slant height for the triangular faces is computed internally.

Formulas and arguments
CallFormulaWhat each argument means
pyramid.volume(s, h)
  • s — side of the square base
  • h — vertical height
Returns the volume
pyramid.lsa(s, h)
  • s — side of the square base
  • h — vertical height (, the slant height, is worked out for you)
Returns the four triangular faces
pyramid.tsa(s, h)
  • s — side of the square base
  • h — vertical height
Returns the four triangles plus the base
pyramid.volume(6, 8)
# 96.0

pyramid.lsa(6, 8)
# 102.52804494381036

pyramid.tsa(6, 8)
# 138.52804494381036

prism.volume(base_area, h) · prism.lsa(base_perimeter, h) · prism.tsa(base_area, base_perimeter, h)

General prism of any base shape — described by its base area and base perimeter, rather than a specific shape.

Formulas and arguments
CallFormulaWhat each argument means
prism.volume(base_area, h)
  • base_area () — area of the base
  • h — height (the length of the prism)
Returns the volume
prism.lsa(base_perimeter, h)
  • base_perimeter () — perimeter of the base
  • h — height
Returns the side faces
prism.tsa(base_area, base_perimeter, h)
  • base_area () — area of the base
  • base_perimeter () — perimeter of the base
  • h — height
Returns the side faces plus both bases
prism.volume(20, 10)
# 200

prism.lsa(18, 10)
# 180

prism.tsa(20, 18, 10)
# 220

torus.volume(R, r) · torus.surface_area(R, r)

Volume and surface area of a torus (donut shape). R is the distance from the tube's center to the torus's center; r is the tube's own radius. No TSA/CSA split — the whole surface is uniform.

Formulas and arguments
CallFormulaWhat each argument means
torus.volume(R, r)
  • R — distance from the centre of the torus to the centre of the tube
  • r — radius of the tube
Returns the volume
torus.surface_area(R, r)
  • R — distance from the centre of the torus to the centre of the tube
  • r — radius of the tube
Returns the surface area
torus.volume(5, 2)
# 394.78417604357435

torus.surface_area(5, 2)
# 394.78417604357435

Equations

eqn.l() solves a system of simultaneous equations — give it one tuple per equation (its coefficients, then the constant) and it works out every unknown. eqn.q() finds the roots of a polynomial of any degree, including complex roots when there's no real solution.

eqn.l(*equations)

Solves n linear equations in n unknowns. Each equation is (coefficients..., constant) — n+1 numbers for n unknowns.

Formulas and arguments
CallFormulaWhat each argument means
eqn.l(*equations)
  • equations () — one tuple per equation: the coefficients of the unknowns, then the constant on the right
Returns the unknowns, in order (, , …)
eqn.l((2, 3, 5), (7, 6, 10))
# (x, y) solving 2x+3y=5, 7x+6y=10

eqn.l((1,1,1,6), (2,-1,1,3), (1,2,-1,2))
# (a, b, c) — scales to any n

eqn.q(*coeffs) · eqn.q.sum(*coeffs) · eqn.q.mul(*coeffs)

Roots of a polynomial of any degree, plus the sum and product of all roots (Vieta's formulas — valid at any degree).

Formulas and arguments
CallFormulaWhat each argument means
eqn.q(*coeffs)
  • coeffs () — the coefficients of , highest power first (higher degrees are solved numerically)
Returns the roots
eqn.q.sum(*coeffs)
  • coeffs () — the coefficients, highest power first
Returns the sum of the roots (Vieta's formulas)
eqn.q.mul(*coeffs)
  • coeffs () — the coefficients, highest power first
Returns the product of the roots (Vieta's formulas)
eqn.q(1, 3, 4)
# roots of x^2 + 3x + 4 = 0

eqn.q(3, 4, 6, 9, 10)
# roots of the quartic

eqn.q.sum(5, 2, 3)
# -0.4

eqn.q.mul(5, 2, 3)
# 0.6

Logs, Exponents & Roots

Everyday algebra helpers — logarithms with any base, exponents in either direction, and roots of any degree, including correct handling of negative numbers for odd roots.

log(x, base=e)

Logarithm of x. Natural log by default; pass a second argument for a different base.

Formulas and arguments
CallFormulaWhat each argument means
log(x, base)
  • x — the number (greater than 0)
  • base () — the base (default , the natural logarithm)
Returns the logarithm
log(math.e)
# 1.0

log(8, 2)
# 3.0

exp(a, b)

a raised to the power b. Works for positive and negative exponents.

Formulas and arguments
CallFormulaWhat each argument means
exp(a, b)
  • a — the base
  • b — the exponent (can be negative: )
Returns raised to the power
exp(2, 5)
# 32

exp(2, -1)
# 0.5

sqrt(no) · curt(no) · nroot(no, n)

Square root, cube root, and the generic nth root. Odd roots of negative numbers work correctly.

Formulas and arguments
CallFormulaWhat each argument means
sqrt(no)
  • no () — the number
Returns the square root
curt(no)
  • no () — the number
Returns the cube root (negative numbers allowed)
nroot(no, n)
  • no () — the number
  • n — which root
Returns the th root (odd roots of negatives are negative)
sqrt(16)
# 4.0

curt(27)
# 3.0

curt(-27)
# -3.0

nroot(16, 2)
# 4.0

Trigonometry

Degrees by default — pass mode="rad" for radians.

sin(angle, mode="deg") · cos(...) · tan(...)

The three primary trigonometric ratios.

Formulas and arguments
CallFormulaWhat each argument means
sin(angle)
  • angle () — the angle, in degrees (or radians with "rad")
  • mode — "deg" (default) or "rad"
Returns the sine
cos(angle)
  • angle () — the angle
  • mode — "deg" or "rad"
Returns the cosine
tan(angle)
  • angle () — the angle
  • mode — "deg" or "rad"
Returns the tangent
(identity)
Returns the Pythagorean identity, true for every angle
sin(30)
# 0.5

cos(60)
# 0.5

tan(45)
# 1.0

cosec(angle, mode="deg") · sec(...) · cot(...)

The three reciprocal trigonometric ratios. Raise a clear error at an undefined point.

Formulas and arguments
CallFormulaWhat each argument means
cosec(angle)
  • angle () — the angle
  • mode — "deg" or "rad"
Returns undefined where
sec(angle)
  • angle () — the angle
  • mode — "deg" or "rad"
Returns undefined where
cot(angle)
  • angle () — the angle
  • mode — "deg" or "rad"
Returns undefined where
cosec(30)
# 2.0

sec(60)
# 2.0

cot(45)
# 1.0

Classical Mechanics

Newton's laws, motion, projectiles, gravity, energy, momentum, circular motion, friction, rotation, oscillations, orbits, pulleys, fluids and elasticity. All SI units and degrees. Most laws solve for any missing variable, show their working with explain=True, and the motion ones can draw charts.

newton.second(F=10, m=5) motion.suvat(u=0, a=9.81, t=3) ...(explain=True)

How the physics functions work. Everything uses SI units (kilogram, metre, second, newton, joule, watt) and angles in degrees, so the numbers you pass in and get back never need converting. Use convert() for other units. Each family has short functions for the common question, such as newton.force(m, a).

Solve for any missing variable. The main laws also take their variables by name: give all but one, and the law finds the last. That is the same law solved three ways. Add explain=True to any of them to print the equation, what you gave, and the result with its unit. Give the wrong number of variables, a negative mass, or an unknown name, and you get a plain sentence saying what to fix.

newton.second(m=5, a=2)                       # 10   (force)
newton.second(F=10, m=5)                      # 2   (acceleration)
newton.second(F=10, a=2)                      # 5   (mass)
newton.second(F=10, m=5, explain=True)        # 2

# Newton's second law:  F = m x a
#   Given:   F = 10 N,  m = 5 kg
#   Result:  a = 2 m/s^2

newton.force(m, a) newton.mass(F, a) newton.accel(F, m) newton.weight(m, g=9.81)

Newton's laws of motion. The second law has three short forms, one for each unknown. weight() multiplies by gravity (9.81 unless you pass g=). net() adds forces on a line, or as (x, y) pairs on a plane; equilibrium() is Newton's first law: it is True when the forces cancel. third() gives the reaction force, components() splits a force at an angle, and resultant() combines two forces with the angle between them.

Every variable is checked: a mass must be greater than zero, and asking for a mass from a zero acceleration says why that cannot be done.

Formulas and arguments
CallFormulaWhat each argument means
newton.force(m, a)
  • m — mass of the object, kg
  • a — acceleration, m/s²
Returns the force , in N
newton.mass(F, a)
  • F — net force, N
  • a — acceleration, m/s² (not zero)
Returns the mass , in kg
newton.accel(F, m)
  • F — net force, N
  • m — mass, kg
Returns the acceleration , in m/s²
newton.second(F=, m=, a=)
  • F — net force, N
  • m — mass, kg
  • a — acceleration, m/s²
Returns the one variable you leave out (give the other two, by name). Add explain=True to print the working
newton.weight(m, g=9.81)
  • m — mass, kg
  • g — gravitational field strength, m/s² (default 9.81, Earth)
Returns the weight , in N
newton.net(*forces)
  • forces () — the forces to add: numbers on a line (positive one way, negative the other) or pairs
Returns the net force, a number or an pair
newton.equilibrium(*forces, tol)
  • forces () — the forces, as for net
  • tol — how close to zero still counts as zero (default )
Returns True if the forces cancel (Newton's first law)
newton.third(force)
  • force () — the force one body exerts on another
Returns the equal and opposite force
newton.components(force, angle, mode)
  • force () — size of the force, N
  • angle () — angle from the axis, in degrees
  • mode — "deg" (default) or "rad"
Returns the pair
newton.resultant(f1, f2, angle=90)
  • f1 () — first force, N
  • f2 () — second force, N
  • angle () — angle between the two forces, degrees (default 90, where this is Pythagoras)
Returns the size of the combined force
newton.lift(m, a, g=9.81)
  • m — mass of the person, kg
  • a — acceleration of the lift, m/s², positive when going up
  • g — gravitational field strength (default 9.81)
Returns the force a scale reads, in N
newton.force(5, 2)                            # 10
newton.mass(10, 2)                            # 5
newton.accel(10, 5)                           # 2
newton.weight(10)                             # 98.1
newton.weight(10, g=1.62)                     # 16.2   (on the Moon)
newton.net(5, -2, 1)                          # 4
newton.net((3, 4), (1, -4))                   # (4, 0)
newton.equilibrium(10, -4, -6)                # True
newton.third(7)                               # -7
newton.components(10, 30)                     # (8.66025, 5)
newton.resultant(3, 4)                        # 5   (at right angles)

motion.suvat(s=, u=, v=, a=, t=) motion.fall(height) motion.stopping_distance(speed, decel)

Straight-line motion with steady acceleration. Five quantities describe it: distance s, start speed u, end speed v, acceleration a and time t. Give any three to suvat() and it returns all five, choosing the right equation for you (all ten combinations work). If the numbers cannot describe real motion, such as stopping before the distance is reached, it says so.

final_speed(u, a, t), distance(u, a, t), time(u, v, a) and accel(u, v, t) are the short forms. fall(height) returns the time and speed of a drop from rest, and stopping_distance() gives the distance to stop at a steady deceleration.

Formulas and arguments
CallFormulaWhat each argument means
motion.suvat(s=, u=, v=, a=, t=)
  • s — distance travelled, m
  • u — starting speed, m/s
  • v — final speed, m/s
  • a — acceleration, m/s² (steady)
  • t — time taken, s
Returns all five values as a dictionary, from any three you give (by name). Fyrefly picks the right equation
motion.final_speed(u, a, t)
  • u — starting speed, m/s
  • a — acceleration, m/s²
  • t — time, s
Returns the final speed
motion.distance(u, a, t)
  • u — starting speed, m/s
  • a — acceleration, m/s²
  • t — time, s
Returns the distance
motion.time(u, v, a)
  • u — starting speed, m/s
  • v — final speed, m/s
  • a — acceleration, m/s²
Returns the time
motion.accel(u, v, t)
  • u — starting speed, m/s
  • v — final speed, m/s
  • t — time, s
Returns the acceleration
motion.speed(distance, time)
  • distance () — distance travelled, m
  • time () — time taken, s
Returns the average speed
motion.fall(height, g=9.81)
  • height () — height dropped from rest, m
  • g — gravitational field strength (default 9.81)
Returns the pair (time, speed at the bottom)
motion.stopping_distance(speed, decel)
  • speed () — starting speed, m/s
  • decel () — deceleration, m/s², as a positive number
Returns the distance to stop
motion.suvat(u=0, a=9.81, t=3)                # {'s': 44.145, 'u': 0, 'v': 29.43, 'a': 9.81, 't': 3}   (a dropped stone)
motion.suvat(u=0, v=20, s=100)                # {'s': 100, 'u': 0, 'v': 20, 'a': 2, 't': 10}   (a car speeding up)
motion.suvat(u=10, v=0, a=-2)                 # {'s': 25, 'u': 10, 'v': 0, 'a': -2, 't': 5}   (braking)
motion.final_speed(0, 2, 5)                   # 10
motion.distance(0, 2, 5)                      # 25
motion.fall(20)                               # (2.01928, 19.8091)   ((time, speed) from 20 m)
motion.stopping_distance(20, 5)               # 40

projectile.range(v, angle, h0=0) .height .time .at(t, ...) .best_angle(v) .path(...) .plot(...)

Projectile motion with no air resistance. Give the launch speed and the angle above the horizontal (and, if it starts above the ground, h0=). range(), height() and time() give the distance, peak and time of flight, at(t, ...) gives the position at any moment, and velocity_at() gives its speed components. best_angle() finds the angle for the longest throw: 45° from the ground, and lower from a height.

path() returns the whole flight as a table with columns t, x and y, so you can use sql() or viz on it. plot() draws the flight with the peak and landing point marked, and compare=[15, 30, 45, 60, 75] draws several angles together. save= writes the chart to a file, and mode="rad" accepts radians.

Formulas and arguments
CallFormulaWhat each argument means
projectile.time(v, angle, h0, g)
  • v — launch speed, m/s
  • angle () — launch angle above the horizontal, degrees
  • h0 () — starting height above the ground, m (default 0)
  • g — gravitational field strength (default 9.81)
Returns the time of flight , in s (from the ground, )
projectile.range(v, angle, h0, g)
  • v — launch speed, m/s
  • angle () — launch angle above the horizontal, degrees
  • h0 () — starting height above the ground, m (default 0)
  • g — gravitational field strength (default 9.81)
Returns the horizontal distance to landing, in m (from the ground, )
projectile.height(v, angle, h0, g)
  • v — launch speed, m/s
  • angle () — launch angle above the horizontal, degrees
  • h0 () — starting height above the ground, m (default 0)
  • g — gravitational field strength (default 9.81)
Returns the greatest height , in m
projectile.at(t, v, angle, h0, g)
  • t — time since launch, s
  • v — launch speed, m/s
  • angle () — launch angle above the horizontal, degrees
  • h0 () — starting height above the ground, m (default 0)
  • g — gravitational field strength (default 9.81)
Returns the position
projectile.velocity_at(t, v, angle, g)
  • t — time since launch, s
  • v — launch speed, m/s
  • angle () — launch angle, degrees
  • g — gravitational field strength
Returns the velocity
projectile.best_angle(v, h0, g)
  • v — launch speed, m/s
  • h0 () — starting height, m (0 gives )
  • g — gravitational field strength
Returns the angle in degrees that gives the longest range
projectile.path(v, angle, h0, g, steps)
  • v — launch speed
  • angle () — launch angle
  • h0 () — starting height
  • g — gravity
  • steps — how many points to calculate (default 60)
Returns a table with columns t, x, y
projectile.plot(v, angle, h0, g, title, save, compare)
  • v — launch speed
  • angle () — launch angle
  • h0 () — starting height
  • g — gravity
  • title — heading for the chart
  • save — file name to also write the chart to
  • compare — a list of other angles to draw together
Returns draws the flight (peak and landing marked)
projectile.range(20, 45)                      # 40.7747
projectile.height(20, 45)                     # 10.1937
projectile.time(20, 45)                       # 2.88321
projectile.range(20, 30)                      # 35.3119
projectile.range(20, 60)                      # 35.3119   (the same: complementary angles)
projectile.best_angle(20, h0=10)              # 39.3207   (from a 10 m cliff)
projectile.at(1, 20, 45)                      # (14.1421, 9.23714)   ((x, y) after 1 s)
df = projectile.path(20, 45)
projectile.plot(20, 45)
projectile.plot(20, 45, compare=[15, 30, 45, 60, 75])
projectile.plot(20, 45): the peak and the landing point are marked
projectile.plot(20, 45): the peak and the landing point are marked
compare=[15, 30, 45, 60, 75]: 30° and 60° land together, 45° goes farthest
compare=[15, 30, 45, 60, 75]: 30° and 60° land together, 45° goes farthest

gravity.force(m1, m2, r) gravity.field(body) gravity.orbit_speed(body, r) gravity.escape_speed(body)

Gravitation and orbits. Wherever a body is needed you can pass its name ("earth", "moon", "mars", "venus", "jupiter", "sun") instead of a mass, and the distance then defaults to its surface. Pass a mass in kilograms and a distance r (from the centre, in metres) for any other body. force() is Newton's law of gravitation, and gravity.law(F=, m1=, m2=, r=) solves it for any missing variable. orbital_period() gives the time for one orbit, and kepler3() scales a known period to a new distance (T² grows as a³). weight_on(body, m) gives your weight there.

Constants: gravity.G is 6.6743e-11, and gravity.bodies lists the built-in masses and radii.

Formulas and arguments
CallFormulaWhat each argument means
gravity.force(m1, m2, r)
  • m1 () — first mass, kg
  • m2 () — second mass, kg
  • r — distance between their centres, m
Returns the pull , in N ( N m²/kg²)
gravity.law(F=, m1=, m2=, r=)
  • F — force, N
  • m1 () — first mass
  • m2 () — second mass
  • r — distance between centres
Returns the variable you leave out
gravity.field(M, r)
  • M — the body's mass in kg, or a name such as "earth"
  • r — distance from its centre, m (for a named body, the default is its surface)
Returns the field strength , in m/s²
gravity.orbit_speed(M, r)
  • M — mass of the body being orbited, or its name
  • r — orbit radius from the centre, m
Returns the circular-orbit speed
gravity.escape_speed(M, r)
  • M — mass of the body, or its name
  • r — distance from the centre (default: its surface)
Returns the escape speed
gravity.orbital_period(M, a)
  • M — mass of the body being orbited, or its name
  • a — orbit radius, m
Returns the period , in s
gravity.kepler3(T1, a1, a2)
  • T1 () — a known period
  • a1 () — the radius that goes with
  • a2 () — the new radius
Returns the new period (same time unit as )
gravity.weight_on(body, m)
  • body — name of the body, e.g. "moon"
  • m — your mass, kg
Returns your weight there, in N
gravity.force(5.972e24, 1, 6.371e6)           # 9.81997   (pull on 1 kg at Earth's surface)
gravity.field('earth')                        # 9.81997   (g at the surface)
gravity.field('moon')                         # 1.62338
gravity.escape_speed('earth')                 # 11186   (m/s)
gravity.orbit_speed('earth')                  # 7909.68   (lowest orbit, m/s)
gravity.orbital_period('earth', 4.2164e7)     # 86164.8   (seconds: geostationary)
gravity.kepler3(1, 1, 1.524)                  # 1.88138   (Mars year, in Earth years)
gravity.weight_on('moon', 60)                 # 97.4029   (newtons)

energy.kinetic(m, v) energy.potential(m, h) energy.work(F, d, angle) energy.power(work, time) energy.spring(k, x)

Work, energy and power. kinetic() is ½mv², potential() is mgh, spring() is the energy stored in a stretched spring, work() is F·d·cos(angle) with the angle between the force and the movement, and power() is work over time (power_fv() is force times speed). speed_after_drop() and height_reached() use conservation of energy, and efficiency() returns a percentage and refuses a useful output larger than the input.

energy.kinetic_law, potential_law, spring_law and power_law solve their equation for any missing variable by name.

Formulas and arguments
CallFormulaWhat each argument means
energy.kinetic(m, v)
  • m — mass, kg
  • v — speed, m/s
Returns the kinetic energy, in J
energy.potential(m, h, g=9.81)
  • m — mass, kg
  • h — height above your zero level, m
  • g — gravitational field strength
Returns the potential energy, in J
energy.spring(k, x)
  • k — spring constant, N/m
  • x — stretch or squeeze, m
Returns the stored energy, in J
energy.work(force, distance, angle=0)
  • force () — size of the force, N
  • distance () — distance moved, m
  • angle () — angle between the force and the movement, degrees (default 0)
Returns the work , in J
energy.power(work, time)
  • work () — work done, J
  • time () — time taken, s
Returns the power , in W
energy.power_fv(force, speed)
  • force () — steady force, N
  • speed () — speed, m/s
Returns the power , in W
energy.speed_after_drop(height, v0, g)
  • height () — height fallen, m
  • v0 () — starting speed, m/s (default 0)
  • g — gravitational field strength
Returns the speed at the bottom
energy.height_reached(speed, g)
  • speed () — speed thrown straight up, m/s
  • g — gravitational field strength
Returns the greatest height
energy.efficiency(useful, supplied)
  • useful — useful energy or work out
  • supplied — energy or work put in (same unit)
Returns the efficiency , as a percentage
energy.kinetic_law(KE=, m=, v=)
  • KE — kinetic energy, J
  • m — mass, kg
  • v — speed, m/s
Returns the variable you leave out
energy.potential_law(PE=, m=, g=, h=)
  • PE — potential energy, J
  • m — mass, kg
  • g — field strength
  • h — height, m
Returns the variable you leave out
energy.spring_law(E=, k=, x=)
  • E — stored energy, J
  • k — spring constant, N/m
  • x — stretch, m
Returns the variable you leave out
energy.power_law(P=, W=, t=)
  • P — power, W
  • W — work, J
  • t — time, s
Returns the variable you leave out
energy.kinetic(2, 3)                          # 9
energy.potential(2, 10)                       # 196.2
energy.spring(200, 0.1)                       # 1
energy.work(10, 5, 60)                        # 25
energy.power(100, 4)                          # 25
energy.speed_after_drop(20)                   # 19.8091   (m/s after falling 20 m)
energy.height_reached(20)                     # 20.3874   (thrown up at 20 m/s)
energy.efficiency(30, 120)                    # 25   (percent)
energy.kinetic_law(KE=9, m=2)                 # 3   (the speed)

momentum.p(m, v) momentum.collide(m1, u1, m2, u2, e=1) momentum.elastic momentum.inelastic momentum.recoil

Momentum and collisions in a line. Speeds are signed: right is positive, left is negative. collide() returns the two speeds after the collision. e is the coefficient of restitution: 1 is perfectly elastic (bouncy), 0 is perfectly inelastic (they stick together). elastic() and inelastic() are the two ends of that range, and energy_lost() tells you how much kinetic energy went missing. Momentum is always conserved, and the chosen e sets how fast they separate.

impulse(), impulse_change(), recoil() (a gun firing a bullet), com() (centre of mass, on a line or in a plane) and rocket() (the speed a rocket gains from burning fuel) round out the family.

Formulas and arguments
CallFormulaWhat each argument means
momentum.p(m, v)
  • m — mass, kg
  • v — velocity, m/s (signed)
Returns the momentum , in kg m/s
momentum.impulse(force, time)
  • force () — steady force, N
  • time () — how long it acts, s
Returns the impulse , in N s
momentum.impulse_change(m, u, v)
  • m — mass, kg
  • u — speed before, m/s
  • v — speed after, m/s
Returns the impulse
momentum.collide(m1, u1, m2, u2, e)
  • m1 () — mass of the first object, kg
  • u1 () — its speed before, m/s (signed: right is positive)
  • m2 () — mass of the second object, kg
  • u2 () — its speed before, m/s
  • e — coefficient of restitution, from 0 (they stick) to 1 (perfectly bouncy, the default)
Returns the speeds afterwards
momentum.elastic(m1, u1, m2, u2)
  • m1 () — mass of the first object, kg
  • u1 () — its speed before, m/s (signed: right is positive)
  • m2 () — mass of the second object, kg
  • u2 () — its speed before, m/s
Returns the speeds afterwards ()
momentum.inelastic(m1, u1, m2, u2)
  • m1 () — mass of the first object, kg
  • u1 () — its speed before, m/s (signed: right is positive)
  • m2 () — mass of the second object, kg
  • u2 () — its speed before, m/s
Returns the one common speed (, they stick)
momentum.energy_lost(m1, u1, m2, u2, e)
  • m1 () — mass of the first object, kg
  • u1 () — its speed before, m/s (signed: right is positive)
  • m2 () — mass of the second object, kg
  • u2 () — its speed before, m/s
  • e — coefficient of restitution (default 1)
Returns the kinetic energy lost, in J
momentum.recoil(m_gun, m_bullet, v_bullet)
  • m_gun () — mass of the gun, kg
  • m_bullet () — mass of the bullet, kg
  • v_bullet () — bullet speed, m/s
Returns the recoil speed of the gun (negative: backwards)
momentum.com(masses, positions)
  • masses () — a list of masses, kg
  • positions () — a list of positions (numbers on a line, or pairs), in the same order
Returns the centre of mass
momentum.rocket(exhaust_speed, m_start, m_end)
  • exhaust_speed () — speed of the exhaust relative to the rocket, m/s
  • m_start () — mass with fuel, kg
  • m_end () — mass after burning, kg
Returns the speed gained
momentum.law(p=, m=, v=)
  • p — momentum
  • m — mass
  • v — velocity
Returns the variable you leave out
momentum.p(2, 3)                              # 6
momentum.collide(2, 3, 1, 0)                  # (1, 4)   (a 2 kg ball hits a 1 kg ball at rest)
momentum.elastic(1, 5, 1, 0)                  # (0, 5)   (equal masses swap speeds)
momentum.inelastic(1, 5, 1, 0)                # 2.5   (they stick: one common speed)
momentum.energy_lost(1, 5, 1, 0, e=0)         # 6.25   (joules)
momentum.recoil(5, 0.01, 400)                 # -0.8
momentum.com([1, 3], [0, 4])                  # 3
momentum.rocket(3000, 1000, 500)              # 2079.44   (exhaust 3000 m/s, half the mass burnt)

circular.accel(v, r) circular.force(m, v, r) circular.banked_speed(r, angle, mu=0) circular.vertical_circle(r)

Circular motion. accel() and force() are the centripetal acceleration v²/r and force mv²/r. speed(), period() and omega() convert between speed, the time for one lap, and angular speed. banked_speed() gives the fastest steady speed round a banked curve (with optional friction mu=), and vertical_circle() returns the smallest speeds that keep an object moving round a vertical loop: at the top, and at the bottom.

Formulas and arguments
CallFormulaWhat each argument means
circular.accel(v, r)
  • v — speed, m/s
  • r — radius of the circle, m
Returns the centripetal acceleration , in m/s²
circular.force(m, v, r)
  • m — mass, kg
  • v — speed, m/s
  • r — radius, m
Returns the inward force , in N
circular.speed(r, period)
  • r — radius, m
  • period () — time for one lap, s
Returns the speed
circular.period(r, v)
  • r — radius, m
  • v — speed, m/s
Returns the time for one lap
circular.omega(period)
  • period () — time for one lap, s
Returns the angular speed , in rad/s
circular.banked_speed(r, angle, mu, g)
  • r — radius of the curve, m
  • angle () — banking angle, degrees
  • mu () — friction coefficient (default 0)
  • g — gravitational field strength
Returns the fastest steady speed on the curve
circular.vertical_circle(r, g)
  • r — radius of the loop, m
  • g — gravitational field strength
Returns the smallest speeds (top, bottom) that stay on the loop
circular.accel(10, 5)                         # 20
circular.force(2, 10, 5)                      # 40
circular.period(5, 5)                         # 6.28319
circular.omega(2)                             # 3.14159   (rad/s for a 2 s lap)
circular.banked_speed(50, 20)                 # 13.3614   (50 m radius, banked 20 degrees)
circular.banked_speed(50, 20, mu=0.3)         # 19.1206
circular.vertical_circle(2)                   # (4.42945, 9.90454)   ((top, bottom))

friction.force(mu, N) friction.incline(angle, mu=0) friction.repose(mu) friction.stopping_distance(v, mu)

Friction and slopes. incline() gives the acceleration down a slope and returns 0 when friction is enough to hold the object still. repose() is the steepest angle at which an object stays put, stopping_distance() is how far something slides before stopping, and slide_time() is the time to slide down a slope of a given length (it tells you if the object would not move at all).

Formulas and arguments
CallFormulaWhat each argument means
friction.force(mu, normal)
  • mu () — coefficient of friction
  • normal () — normal force pressing the surfaces together, N
Returns the friction force , in N
friction.incline(angle, mu, g)
  • angle () — angle of the slope, degrees
  • mu () — friction coefficient (default 0)
  • g — gravitational field strength
Returns the acceleration down the slope (0 if friction holds it)
friction.repose(mu)
  • mu () — coefficient of friction
Returns the steepest angle in degrees at which it stays put
friction.stopping_distance(v, mu, g)
  • v — starting speed, m/s
  • mu () — friction coefficient
  • g — gravitational field strength
Returns the sliding distance
friction.slide_time(length, angle, mu, g)
  • length () — length of the slope, m
  • angle () — angle of the slope, degrees
  • mu () — friction coefficient
  • g — gravitational field strength
Returns the time to slide down ( is from the incline formula)
friction.force(0.5, 20)                       # 10
friction.incline(30)                          # 4.905   (a smooth slope)
friction.incline(30, 0.2)                     # 3.20586
friction.incline(10, 0.5)                     # 0   (friction wins: 0)
friction.repose(1)                            # 45   (degrees)
friction.stopping_distance(20, 0.8)           # 25.4842
friction.slide_time(5, 30)                    # 1.42784

inertia.sphere(m, r) rotation.torque(F, arm, angle=90) rotation.rolling_speed(height, shape)

Rotation. inertia holds the moment of inertia of common shapes: point, ring, disc, sphere, hollow_sphere, rod_center, rod_end and plate, plus shift() for the parallel axis theorem. rotation has torque(), angular_accel(), angular_momentum(), kinetic() and lever() (the effort that balances a load). rolling_speed() gives the speed at the bottom of a slope for something that rolls without slipping: a sphere beats a disc, which beats a ring, and all of them lose to sliding.

Formulas and arguments
CallFormulaWhat each argument means
inertia.point(m, r)
  • m — mass, kg
  • r — distance from the axis, m
Returns the moment of inertia , in kg m²
inertia.ring(m, r)
  • m — mass, kg
  • r — radius of the ring, m
Returns about its axis
inertia.disc(m, r)
  • m — mass, kg
  • r — radius, m
Returns of a solid disc or cylinder
inertia.sphere(m, r)
  • m — mass, kg
  • r — radius, m
Returns of a solid sphere about a diameter
inertia.hollow_sphere(m, r)
  • m — mass, kg
  • r — radius, m
Returns of a thin shell
inertia.rod_center(m, length)
  • m — mass, kg
  • length () — length of the rod, m
Returns about the middle of the rod
inertia.rod_end(m, length)
  • m — mass, kg
  • length () — length of the rod, m
Returns about one end
inertia.plate(m, a, b)
  • m — mass, kg
  • a — one side, m
  • b — the other side, m
Returns about the centre, axis perpendicular to the plate
inertia.shift(I_cm, m, d)
  • I_cm () — moment of inertia about the centre of mass
  • m — mass, kg
  • d — distance between the two parallel axes, m
Returns about the new axis (parallel axis theorem)
rotation.torque(force, arm, angle=90)
  • force () — the force, N
  • arm () — distance from the pivot to where the force acts, m
  • angle () — angle between the arm and the force, degrees (default 90)
Returns the torque , in N m
rotation.angular_accel(torque, I)
  • torque () — torque, N m
  • I — moment of inertia, kg m²
Returns the angular acceleration , in rad/s²
rotation.torque_law(torque=, I=, alpha=)
  • torque () — torque
  • I — moment of inertia
  • alpha () — angular acceleration
Returns the variable you leave out
rotation.angular_momentum(I, omega)
  • I — moment of inertia, kg m²
  • omega () — angular speed, rad/s
Returns the angular momentum
rotation.kinetic(I, omega)
  • I — moment of inertia
  • omega () — angular speed, rad/s
Returns the rotational kinetic energy, in J
rotation.rolling_speed(height, shape, g)
  • height () — height of the drop, m
  • shape () — "sphere", "hollow_sphere", "disc", "cylinder" or "ring"; is
  • g — gravitational field strength
Returns the speed at the bottom
rotation.lever(load, load_arm, effort_arm)
  • load — the load, N
  • load_arm — its distance from the pivot, m
  • effort_arm — distance from the pivot to where you push, m
Returns the effort force needed, in N
inertia.sphere(5, 2)                          # 8
inertia.rod_center(6, 2)                      # 2
inertia.shift(inertia.rod_center(6, 2), 6, 1) # 8   (= rod_end(6, 2))
rotation.torque(10, 0.5)                      # 5
rotation.torque(10, 0.5, 30)                  # 2.5
rotation.angular_accel(5, 2.5)                # 2
rotation.rolling_speed(1, 'sphere')           # 3.74357   (m/s after a 1 m drop)
rotation.rolling_speed(1, 'ring')             # 3.13209
rotation.lever(100, 0.5, 2)                   # 25

shm.spring_period(m, k) shm.pendulum_period(L) shm.position(A, period, t) shm.plot(A, period)

Springs, pendulums and simple harmonic motion. spring_period() and pendulum_period() give the time for one swing; frequency() and omega() convert it. position() and velocity() give the motion at any time, max_speed() the fastest point, energy() the total energy, and hooke() the restoring force of a spring. path() returns the motion as a table of t and x, and plot() draws it.

Formulas and arguments
CallFormulaWhat each argument means
shm.spring_period(m, k)
  • m — mass, kg
  • k — spring constant, N/m
Returns the period , in s
shm.pendulum_period(length, g)
  • length () — length of the string, m (small swings)
  • g — gravitational field strength
Returns the period
shm.hooke(k, x)
  • k — spring constant, N/m
  • x — displacement from the middle, m
Returns the restoring force (negative: it pulls back)
shm.frequency(period)
  • period () — time for one swing, s
Returns the frequency , in Hz
shm.omega(period)
  • period () — time for one swing, s
Returns the angular frequency , in rad/s
shm.position(A, period, t, phase)
  • A — amplitude, m
  • period () — period, s
  • t — time, s
  • phase () — starting phase, radians (default 0)
Returns the displacement
shm.velocity(A, period, t, phase)
  • A — amplitude, m
  • period () — period, s
  • t — time, s
  • phase () — starting phase (default 0)
Returns the velocity
shm.max_speed(A, period)
  • A — amplitude, m
  • period () — period, s
Returns the greatest speed, at the middle
shm.energy(k, A)
  • k — spring constant, N/m
  • A — amplitude, m
Returns the total energy , in J
shm.path(A, period, cycles, steps) / shm.plot(...)
  • A — amplitude
  • period () — period
  • cycles — how many cycles to cover (default 2)
  • steps — points to calculate
  • title — chart heading
  • save — file name to write the chart to
Returns a table of t, x; or a chart
shm.spring_period(0.5, 200)                   # 0.314159
shm.pendulum_period(1)                        # 2.00607
shm.frequency(2)                              # 0.5
shm.position(0.1, 2, 1)                       # -0.1   (half a period in)
shm.max_speed(0.1, 2)                         # 0.314159
shm.energy(100, 0.1)                          # 0.5
shm.hooke(10, 0.3)                            # -3
shm.plot(0.2, 2)
shm.plot(0.2, 2): a 0.2 m swing with a 2 second period, over two cycles
shm.plot(0.2, 2): a 0.2 m swing with a 2 second period, over two cycles

damping.position(x0, m, k, c, t) damping.regime(m, k, c) damping.driven_amplitude(F0, m, k, c, w) damping.plot(...)

Damped and driven oscillations. A mass on a spring with friction c. damping_ratio() tells you which of three behaviours you get, and regime() names it: underdamped (it rings and fades), critically damped (back to rest fastest, no overshoot) or overdamped (it creeps back). position() gives the displacement at any time for all three, starting from any position and speed. omega_damped(), period(), q_factor() and decay_time() describe the ringing.

Push it at a steady rhythm and you have a driven oscillator. driven_amplitude() is the swing size at a driving frequency w, resonance_frequency() is the frequency that swings it most, and peak_amplitude() the biggest swing. plot() draws the motion with its decaying envelope, plot_resonance() draws the resonance curve, and path() returns the motion as a table.

Formulas and arguments
CallFormulaWhat each argument means
damping.omega0(m, k)
  • m — mass, kg
  • k — spring constant, N/m
Returns the natural frequency , rad/s
damping.gamma(m, c)
  • m — mass, kg
  • c — damping constant, N s/m
Returns the decay rate , 1/s
damping.damping_ratio(m, k, c)
  • m — mass, kg
  • k — spring constant, N/m
  • c — damping constant, N s/m
Returns the damping ratio
damping.critical(m, k)
  • m — mass, kg
  • k — spring constant, N/m
Returns the critical damping constant, N s/m
damping.regime(m, k, c)
  • m — mass
  • k — spring constant
  • c — damping constant
Returns "underdamped", "critically damped", "overdamped", or "undamped"
damping.omega_damped(m, k, c)
  • m — mass
  • k — spring constant
  • c — damping constant
Returns the ringing frequency (underdamped only)
damping.period(m, k, c)
  • m — mass
  • k — spring constant
  • c — damping constant
Returns the time of one ring
damping.q_factor(m, k, c)
  • m — mass
  • k — spring constant
  • c — damping constant
Returns the quality factor
damping.decay_time(m, c)
  • m — mass, kg
  • c — damping constant, N s/m
Returns the time for the amplitude to fall to
damping.amplitude(A, m, c, t)
  • A () — starting amplitude, m
  • m — mass
  • c — damping constant
  • t — time, s
Returns the envelope size at time
damping.position(x0, m, k, c, t, v0)
  • x0 () — starting displacement, m
  • m — mass
  • k — spring constant
  • c — damping constant
  • t — time, s
  • v0 () — starting speed, m/s (default 0)
Returns the displacement at time (underdamped, critical, overdamped, in that order)
damping.driven_amplitude(F0, m, k, c, w)
  • F0 () — size of the driving force, N
  • m — mass
  • k — spring constant
  • c — damping constant
  • w () — driving angular frequency, rad/s
Returns the steady swing size , in m
damping.resonance_frequency(m, k, c)
  • m — mass
  • k — spring constant
  • c — damping constant
Returns the driving frequency of the biggest swing (0 if there is no peak)
damping.peak_amplitude(F0, m, k, c)
  • F0 () — driving force, N
  • m — mass
  • k — spring constant
  • c — damping constant
Returns the biggest swing
damping.path(x0, m, k, c, v0, duration, steps) / plot(...) / plot_resonance(F0, m, k, c)
  • x0 () — start displacement
  • m,k,c — mass, spring constant, damping constant
  • v0 () — start speed
  • duration — how long to cover, s
  • F0 () — driving force (for the resonance curve)
Returns a table of t, x; or a chart
damping.damping_ratio(1, 100, 4)              # 0.2
damping.regime(1, 100, 4)                     # 'underdamped'
damping.regime(1, 100, 20)                    # 'critically damped'
damping.critical(1, 100)                      # 20
damping.omega_damped(1, 100, 4)               # 9.79796   (rad/s)
damping.q_factor(1, 100, 4)                   # 2.5
damping.decay_time(1, 4)                      # 0.5   (seconds)
damping.position(1, 1, 100, 4, 0.5)           # -0.00554445   (starts 1 m out)
damping.driven_amplitude(10, 1, 100, 4, 10)   # 0.25
damping.resonance_frequency(1, 100, 4)        # 9.59166
damping.peak_amplitude(10, 1, 100, 4)         # 0.255155
damping.plot(1, 1, 100, 4)
damping.plot_resonance(10, 1, 100, 4)
damping.plot(1, 1, 100, 4): the swings shrink inside their envelope
damping.plot(1, 1, 100, 4): the swings shrink inside their envelope
damping.plot_resonance(10, 1, 100, 4): the biggest swing is just below the natural frequency
damping.plot_resonance(10, 1, 100, 4): the biggest swing is just below the natural frequency

orbit.period(M, a) orbit.speed(M, r, a) orbit.hohmann(M, r1, r2) orbit.path(M, a, e) orbit.plot(M, a, e)

Orbits and Kepler's laws. Give the body you orbit ("earth", "sun", or a mass in kg), the semi-major axis a in metres, and the eccentricity e (0 is a circle, close to 1 is a long thin ellipse). periapsis() and apoapsis() give the closest and farthest distances, period() is Kepler's third law, and speed() is the speed at any distance (vis-viva), fastest when closest. areal_velocity() is the constant area swept per second (Kepler's second law), and total_energy() is negative for any bound orbit. orbit.AU is 1.496×10¹¹ m.

hohmann() plans the cheapest transfer between two circular orbits: the two speed changes, their total, and the coasting time. path() returns one orbit as a table at equal time steps, and plot() draws it with the body at a focus. The dots are equally spaced in time, so they bunch up where the body is slow (Kepler's second law).

Formulas and arguments
CallFormulaWhat each argument means
orbit.periapsis(a, e)
  • a — semi-major axis (half the long width of the ellipse), m
  • e — eccentricity, 0 (circle) to just below 1
Returns the closest distance
orbit.apoapsis(a, e)
  • a — semi-major axis, m
  • e — eccentricity
Returns the farthest distance
orbit.semi_major(r_peri, r_apo)
  • r_peri () — closest distance, m
  • r_apo () — farthest distance, m
Returns the semi-major axis
orbit.eccentricity(r_peri, r_apo)
  • r_peri () — closest distance
  • r_apo () — farthest distance
Returns the eccentricity
orbit.period(M, a)
  • M — mass of the body orbited, kg, or its name (e.g. "sun")
  • a — semi-major axis, m
Returns the period , in s
orbit.speed(M, r, a)
  • M — mass of the body orbited, or its name
  • r — the distance you are at, m
  • a — semi-major axis, m
Returns the speed there (vis-viva)
orbit.total_energy(M, m, a)
  • M — mass of the body orbited
  • m — mass of the orbiting object, kg
  • a — semi-major axis, m
Returns the total energy , in J (negative: bound)
orbit.areal_velocity(M, a, e)
  • M — mass of the body orbited
  • a — semi-major axis, m
  • e — eccentricity (default 0)
Returns the area swept per second (Kepler's second law)
orbit.hohmann(M, r1, r2)
  • M () — the body orbited (mass or name);
  • r1 () — radius of the starting circular orbit, m
  • r2 () — radius of the final circular orbit, m
Returns a dictionary with dv1, dv2, total (m/s) and time (s, the coasting time)
orbit.path(M, a, e, steps) / plot(M, a, e, steps)
  • M — body orbited
  • a — semi-major axis
  • e — eccentricity
  • steps — points (or dots) at equal time steps
Returns a table of t, x, y, speed; or a chart
orbit.period('sun', orbit.AU) / 86400         # 365.203   (days in Earth's year)
orbit.periapsis(1e7, 0.5)                     # 5e+06
orbit.apoapsis(1e7, 0.5)                      # 1.5e+07
orbit.eccentricity(5e6, 1.5e7)                # 0.5
orbit.speed('earth', 5e6, 1e7)                # 10935.1   (at the closest point)
orbit.speed('earth', 1.5e7, 1e7)              # 3645.04   (at the farthest point)
orbit.total_energy('earth', 1000, 1e7)        # -1.99295e+10   (joules)
orbit.hohmann('earth', 6.771e6, 4.2164e7)     # {'dv1': 2399.43, 'dv2': 1457.2, 'total': 3856.63, 'time': 19044.6}   (low orbit to geostationary)
orbit.path('earth', 1e7, 0.5, steps=8)
             t             x             y         speed
0     0.000000  5.000000e+06  0.000000e+00  10935.115857
1  1244.019305 -1.958050e+06  8.249843e+06   7359.264654
2  2488.038611 -9.351309e+06  7.797409e+06   5061.049936
3  3732.057916 -1.361876e+07  4.391779e+06   3981.366172
4  4976.077221 -1.500000e+07  1.060575e-09   3645.038619
5  6220.096527 -1.361876e+07 -4.391779e+06   3981.366172
6  7464.115832 -9.351309e+06 -7.797409e+06   5061.049936
7  8708.135137 -1.958050e+06 -8.249843e+06   7359.264654
8  9952.154443  5.000000e+06 -2.121150e-09  10935.115857
orbit.plot('earth', 1e7, 0.5)
orbit.plot('earth', 1e7, 0.5): equal time steps, so the dots bunch up far from the centre
orbit.plot('earth', 1e7, 0.5): equal time steps, so the dots bunch up far from the centre

pulley.atwood(m1, m2) pulley.on_table(m_table, m_hang, mu) pulley.on_incline(m_slope, angle, m_hang, mu)

Pulleys and connected masses. A light string over a smooth pulley. Each function returns (acceleration, tension), and a positive acceleration means the hanging mass goes down. atwood() is two hanging masses, on_table() is a block on a table pulled by a hanging mass, and on_incline() is a block on a slope. Friction mu is optional, and the answer is 0 when friction is enough to hold everything still. On a slope the acceleration comes out negative if the block slides down and drags the hanging mass up.

newton.lift(m, a) (and pulley.lift_tension()) give the apparent weight in a lift and the cable tension; a is positive when it accelerates upwards.

Formulas and arguments
CallFormulaWhat each argument means
pulley.atwood(m1, m2, g)
  • m1 () — one hanging mass, kg
  • m2 () — the other hanging mass, kg
  • g — gravitational field strength
Returns (acceleration, tension); positive means goes down
pulley.on_incline(m_slope, angle, m_hang, mu, g)
  • m_slope () — mass of the block on the slope, kg
  • angle () — slope angle, degrees
  • m_hang () — the hanging mass, kg
  • mu () — friction coefficient (default 0)
  • g — gravitational field strength
Returns (acceleration, tension); positive means the hanging mass goes down, 0 if friction holds
pulley.on_table(m_table, m_hang, mu, g)
  • m_table () — mass of the block on the table, kg
  • m_hang () — the hanging mass, kg
  • mu () — friction coefficient (default 0)
  • g — gravitational field strength
Returns (acceleration, tension)
pulley.lift_tension(m, a, g)
  • m — mass of the lift and load, kg
  • a — acceleration, m/s², positive going up
  • g — gravitational field strength
Returns the cable tension , in N
pulley.atwood(3, 2)                           # (1.962, 23.544)   (acceleration, then tension)
pulley.atwood(2, 2)                           # (0, 19.62)   (balanced)
pulley.on_table(2, 1)                         # (3.27, 6.54)
pulley.on_table(2, 1, mu=1)                   # (0, 9.81)   (friction wins)
pulley.on_incline(5, 30, 2)                   # (-0.700714, 21.0214)   (the block slides down)
pulley.on_incline(5, 30, 2, mu=0.5)           # (0, 19.62)   (held by friction)
newton.lift(70, 2)                            # 826.7   (a 70 kg person, lift speeding up)
newton.lift(70, -2)                           # 546.7   (slowing)

fluid.pressure(F, A) fluid.hydrostatic(depth) fluid.buoyancy(volume) fluid.continuity(A1, v1, A2) fluid.bernoulli(...)

Fluids. Water (1000 kg/m³) is the default fluid; pass rho= for another. pressure() is force over area, hydrostatic() and absolute_pressure() give the pressure at a depth, and hydraulic() is a hydraulic press. buoyancy(), floats() and submerged_fraction() cover floating and sinking. continuity(), flow_rate(), bernoulli() and torricelli() describe flowing liquid, reynolds() tells smooth flow from turbulent, and terminal_velocity() is a small sphere sinking through a thick liquid. fluid.pressure_law solves P = F/A for any missing variable.

Formulas and arguments
CallFormulaWhat each argument means
fluid.pressure(F, A)
  • F — force, N
  • A — area it is spread over, m²
Returns the pressure , in Pa
fluid.pressure_law(P=, F=, A=)
  • P — pressure
  • F — force
  • A — area
Returns the variable you leave out
fluid.hydrostatic(depth, rho, g)
  • depth () — depth below the surface, m
  • rho () — density of the fluid, kg/m³ (default 1000, water)
  • g — gravitational field strength
Returns the extra pressure
fluid.absolute_pressure(depth, rho, g)
  • depth () — depth, m
  • rho () — fluid density
  • g — gravity
Returns the total pressure ( Pa)
fluid.hydraulic(F1, A1, A2)
  • F1 () — force on the first piston, N
  • A1 () — area of the first piston
  • A2 () — area of the second piston
Returns the force on the second piston
fluid.buoyancy(volume, rho, g)
  • volume () — volume under the surface, m³
  • rho () — density of the fluid
  • g — gravity
Returns the upthrust , in N
fluid.floats(rho_object, rho_fluid)
  • rho_object () — density of the object
  • rho_fluid () — density of the fluid (default 1000)
Returns True if it floats
fluid.submerged_fraction(rho_object, rho_fluid)
  • rho_object () — density of the object
  • rho_fluid () — density of the fluid
Returns the fraction under the surface (1 if it sinks)
fluid.continuity(A1, v1, A2)
  • A1 () — cross-section of the first pipe, m²
  • v1 () — speed in it, m/s
  • A2 () — cross-section of the second pipe, m²
Returns the speed in the second pipe
fluid.flow_rate(A, v)
  • A — cross-section area, m²
  • v — speed of the fluid, m/s
Returns the volume per second , in m³/s
fluid.bernoulli(P1, v1, h1, v2, h2, rho, g)
  • P1 () — pressure at the first point, Pa
  • v1 () — speed there, m/s
  • h1 () — height there, m
  • v2 () — speed at the second point
  • h2 () — height at the second point
  • rho () — fluid density (default 1000)
  • g — gravity
Returns the pressure at the second point
fluid.torricelli(depth, g)
  • depth () — depth of the hole below the surface, m
  • g — gravity
Returns the jet speed
fluid.reynolds(rho, v, length, viscosity)
  • rho () — fluid density, kg/m³
  • v — flow speed, m/s
  • length () — pipe diameter or object size, m
  • viscosity () — viscosity, Pa s
Returns the Reynolds number (smooth below ~2000, turbulent above ~4000)
fluid.terminal_velocity(radius, rho_sphere, rho_fluid, viscosity, g)
  • radius () — radius of the sphere, m
  • rho_sphere () — its density
  • rho_fluid () — density of the fluid
  • viscosity () — viscosity, Pa s
  • g — gravity
Returns the steady sinking speed
fluid.pressure(100, 0.5)                      # 200
fluid.hydrostatic(10)                         # 98100   (10 m of water)
fluid.absolute_pressure(10)                   # 199425
fluid.hydraulic(10, 0.01, 1)                  # 1000   (a small push, a big force)
fluid.buoyancy(0.01)                          # 98.1   (a 10 litre volume)
fluid.floats(900)                             # True   (oil)
fluid.submerged_fraction(900)                 # 0.9
fluid.continuity(2, 3, 1)                     # 6   (the pipe narrows)
fluid.bernoulli(2e5, 2, 0, 4, 0)              # 194000   (pressure falls as it speeds up)
fluid.torricelli(5)                           # 9.90454
fluid.reynolds(1000, 1, 0.1, 1e-3)            # 100000

elasticity.stress(F, A) elasticity.strain(extension, length) elasticity.youngs(F, A, L, dL) elasticity.series(*ks)

Elasticity. How solids stretch. stress() is force per area, strain() is extension per length, and youngs() is their ratio, the stiffness of the material. extension() tells you how much a rod stretches under a load, elastic_energy() the energy stored, and series() and parallel() combine springs. shear_modulus() and bulk_modulus() cover twisting and squeezing. elasticity.youngs_law solves E = F L / (A dL) for any one missing variable.

Formulas and arguments
CallFormulaWhat each argument means
elasticity.stress(F, A)
  • F — force pulling, N
  • A — cross-section area, m²
Returns the stress , in Pa
elasticity.strain(extension, length)
  • extension () — how much it stretches, m
  • length () — its original length, m
Returns the strain (no unit)
elasticity.youngs(F, A, L, dL)
  • F — force, N
  • A — cross-section area, m²
  • L — original length, m
  • dL () — extension, m
Returns Young's modulus , in Pa
elasticity.youngs_law(E=, F=, A=, L=, dL=)
  • E — Young's modulus
  • F — force
  • A — area
  • L — length
  • dL () — extension
Returns the variable you leave out
elasticity.extension(F, L, A, E)
  • F — force, N
  • L — original length, m
  • A — cross-section area, m²
  • E — Young's modulus, Pa
Returns the extension , in m
elasticity.elastic_energy(F, extension)
  • F — the final force, N
  • extension () — the extension, m
Returns the stored energy , in J
elasticity.series(*ks)
  • ks () — the spring constants, N/m, separated by commas
Returns the combined of springs end to end
elasticity.parallel(*ks)
  • ks () — the spring constants, N/m
Returns the combined of springs side by side
elasticity.shear_modulus(F, A, shear_strain)
  • F — sideways force, N
  • A — area, m²
  • shear_strain () — the shear strain
Returns the shear modulus
elasticity.bulk_modulus(dP, dV, V)
  • dP () — change in pressure, Pa
  • dV () — change in volume (negative when squeezed)
  • V — original volume
Returns the bulk modulus
elasticity.stress(100, 1e-6)                  # 1e+08
elasticity.strain(1e-3, 2)                    # 0.0005
elasticity.youngs(100, 1e-6, 2, 1e-3)         # 2e+11   (steel, about 2e11)
elasticity.extension(100, 2, 1e-6, 2e11)      # 0.001
elasticity.youngs_law(E=2e11, F=100, A=1e-6, L=2) # 0.001   (the extension, in metres)
elasticity.elastic_energy(100, 1e-3)          # 0.05
elasticity.series(10, 10)                     # 5
elasticity.parallel(10, 10)                   # 20

Quantum & Atomic Physics

Light as photons, matter waves, the Bohr atom, particles in boxes, tunnelling, the nucleus, radioactive decay and X-rays. Energies and sizes accept unit= and length= so you can work in eV and nm, and every function is shown with the formula it uses.

photon.energy(wavelength=500, length="nm", unit="eV") constants.h

How the quantum functions work. Like the mechanics functions, everything is in SI units by default: joules, metres, kilograms, seconds, kelvin, and angles in degrees. Atomic physics is usually quoted in electron-volts and nanometres, so any function that takes or returns an energy accepts unit= ("J", "eV", "keV", "MeV"), and any function that takes or returns a wavelength or size accepts length= ("m", "cm", "mm", "um", "nm", "pm", "angstrom", "fm"). The unit applies to the numbers you give and the number you get back. A wrong unit name says which ones are allowed.

constants holds the fundamental constants: constants.h, hbar, c, e, m_e, m_p, m_n, u, k_B, a0, R_inf, alpha, sigma, b_wien, Ry, and others. constants.names() lists them all.

photon.energy(wavelength=500, length="nm", unit="eV") # 2.47968   (green light, in eV)
photon.energy(wavelength=500e-9)              # 3.97289e-19   (the same, in joules)
constants.h                                   # 6.62607e-34
constants.c                                   # 2.99792e+08
constants.names()                             # ('N_A', 'R', 'R_inf', 'Ry', 'a0', 'alpha', 'atm', 'b_wien', 'c', 'e', 'eV', 'eps0', 'h', 'hbar', 'k_B', 'k_e', 'm_H', 'm_e', 'm_n', 'm_neutron_u', 'm_p', 'mu0', 'sigma', 'u', 'u_MeV')

photon.energy(wavelength=, frequency=) photon.wavelength(energy=) photon.photoelectric(work_function, wavelength=)

Light as photons. Light comes in packets, each carrying energy proportional to its frequency. Give energy() a wavelength or a frequency, frequency() a wavelength or an energy, and wavelength() an energy or a frequency. Give exactly one; two is an error. momentum() is the momentum of a photon.

photoelectric() is Einstein's photoelectric effect: the greatest kinetic energy of electrons knocked out of a metal with a given work function. If the photon does not have enough energy, no electrons come out and the answer is 0, however bright the light. threshold_wavelength() and threshold_frequency() give the cut-off, and stopping_voltage() the voltage that just stops the fastest electrons. compton_shift() is how much a photon's wavelength grows when it scatters off an electron.

Formulas and arguments
CallFormulaWhat each argument means
photon.energy(wavelength=, frequency=, unit, length)
  • wavelength () — wavelength of the light
  • frequency () — frequency, Hz (give exactly one of wavelength or frequency)
  • unit — energy unit for what you give and what you get back: "J" (default), "eV", "keV" or "MeV"
  • length — length unit for wavelengths and sizes you give and get back: "m" (default), "cm", "mm", "um", "nm", "pm", "angstrom", "fm"
Returns the energy of one photon ( J s)
photon.frequency(wavelength=, energy=, unit, length)
  • wavelength () — wavelength
  • energy () — energy of the photon (give exactly one)
  • unit — energy unit for what you give and what you get back: "J" (default), "eV", "keV" or "MeV"
  • length — length unit for wavelengths and sizes you give and get back: "m" (default), "cm", "mm", "um", "nm", "pm", "angstrom", "fm"
Returns the frequency , in Hz
photon.wavelength(energy=, frequency=, unit, length)
  • energy () — energy of the photon
  • frequency () — frequency, Hz (give exactly one)
  • unit — energy unit for what you give and what you get back: "J" (default), "eV", "keV" or "MeV"
  • length — length unit for wavelengths and sizes you give and get back: "m" (default), "cm", "mm", "um", "nm", "pm", "angstrom", "fm"
Returns the wavelength
photon.momentum(wavelength, length)
  • wavelength () — wavelength of the light
  • length — length unit for wavelengths and sizes you give and get back: "m" (default), "cm", "mm", "um", "nm", "pm", "angstrom", "fm"
Returns the photon's momentum , in kg m/s
photon.photoelectric(work_function, wavelength=, frequency=, unit, length)
  • work_function () — work function of the metal (energy needed to free an electron), in the unit energy unit
  • wavelength () — wavelength of the light (or give the frequency)
  • frequency () — frequency, Hz
  • unit — energy unit for what you give and what you get back: "J" (default), "eV", "keV" or "MeV"
  • length — length unit for wavelengths and sizes you give and get back: "m" (default), "cm", "mm", "um", "nm", "pm", "angstrom", "fm"
Returns the greatest kinetic energy of the electrons; 0 if no electrons come out
photon.threshold_frequency(work_function, unit)
  • work_function () — work function of the metal (energy needed to free an electron), in the unit energy unit
  • unit — energy unit for what you give and what you get back: "J" (default), "eV", "keV" or "MeV"
Returns the lowest frequency that frees electrons
photon.threshold_wavelength(work_function, unit, length)
  • work_function () — work function of the metal (energy needed to free an electron), in the unit energy unit
  • unit — energy unit for what you give and what you get back: "J" (default), "eV", "keV" or "MeV"
  • length — length unit for wavelengths and sizes you give and get back: "m" (default), "cm", "mm", "um", "nm", "pm", "angstrom", "fm"
Returns the longest wavelength that frees electrons
photon.stopping_voltage(work_function, wavelength=, frequency=, unit, length)
  • work_function () — work function of the metal (energy needed to free an electron), in the unit energy unit
  • wavelength () — wavelength of the light (or the frequency)
  • frequency () — frequency
  • unit — energy unit for what you give and what you get back: "J" (default), "eV", "keV" or "MeV"
  • length — length unit for wavelengths and sizes you give and get back: "m" (default), "cm", "mm", "um", "nm", "pm", "angstrom", "fm"
Returns the voltage that just stops the fastest electrons, in volts
photon.compton_wavelength(mass, length)
  • mass () — mass of the particle, kg (default: an electron)
  • length — length unit for wavelengths and sizes you give and get back: "m" (default), "cm", "mm", "um", "nm", "pm", "angstrom", "fm"
Returns the Compton wavelength
photon.compton_shift(angle, mass, length)
  • angle () — angle the photon scatters through, degrees
  • mass () — mass of the particle hit, kg (default: an electron)
  • length — length unit for wavelengths and sizes you give and get back: "m" (default), "cm", "mm", "um", "nm", "pm", "angstrom", "fm"
Returns how much the wavelength grows,
photon.energy(wavelength=500, length="nm", unit="eV") # 2.47968
photon.wavelength(energy=2.5, unit="eV", length="nm") # 495.937
photon.frequency(wavelength=500e-9)           # 5.99585e+14
photon.momentum(500, length="nm")             # 1.32521e-27
photon.photoelectric(2.28, wavelength=400, unit="eV", length="nm") # 0.819605   (sodium, violet light)
photon.photoelectric(2.28, wavelength=700, unit="eV", length="nm") # 0   (red: nothing comes out)
photon.threshold_wavelength(2.28, unit="eV", length="nm") # 543.79
photon.stopping_voltage(2.28, wavelength=400, unit="eV", length="nm") # 0.819605
photon.compton_wavelength(length="pm")        # 2.42631
photon.compton_shift(90, length="pm")         # 2.42631

photon.wien(T) photon.stefan(T, area) photon.planck(wavelength, T) photon.plot_blackbody(temperatures)

Glowing bodies. Hot objects glow, and the colour depends on the temperature. wien() gives the wavelength where the glow is brightest, stefan() the total power radiated (it grows as the fourth power of the temperature), and planck() the brightness at one wavelength. plot_blackbody() draws Planck's curves for several temperatures, with the peaks marked.

Formulas and arguments
CallFormulaWhat each argument means
photon.wien(temperature, length)
  • temperature () — temperature, kelvin
  • length — length unit for wavelengths and sizes you give and get back: "m" (default), "cm", "mm", "um", "nm", "pm", "angstrom", "fm"
Returns the wavelength of the brightest glow ( m K)
photon.stefan(temperature, area, emissivity)
  • temperature () — temperature, kelvin
  • area () — radiating area, m² (default 1)
  • emissivity () — 0 to 1 (default 1, a perfect black body)
Returns the power radiated, in W ( W/(m² K⁴))
photon.planck(wavelength, temperature, length)
  • wavelength () — the wavelength you ask about
  • temperature () — temperature, kelvin
  • length — length unit for wavelengths and sizes you give and get back: "m" (default), "cm", "mm", "um", "nm", "pm", "angstrom", "fm"
Returns the brightness at that wavelength, in W/(m² m sr) ( J/K)
photon.plot_blackbody(temperatures, title, save)
  • temperatures () — a list of temperatures, kelvin (default 3000, 4500, 6000)
  • title — chart heading
  • save — file name to write the chart to
Returns draws the Planck curves with the Wien peaks marked
photon.wien(5772, length="nm")                # 502.039   (the Sun's surface: green-yellow)
photon.wien(300, length="um")                 # 9.65924   (a person: infrared)
photon.stefan(5772, 6.09e18)                  # 3.83296e+26   (watts from the Sun's area)
photon.stefan(300, area=1.8, emissivity=0.97) # 801.938   (a person)
photon.plot_blackbody([3000, 4500, 6000])
photon.plot_blackbody([3000, 4500, 6000]): hotter bodies are brighter, and their peak moves towards blue
photon.plot_blackbody([3000, 4500, 6000]): hotter bodies are brighter, and their peak moves towards blue

matter.de_broglie(mass, speed) matter.electron_wavelength(volts) matter.momentum_spread(dx)

Matter waves and uncertainty. Particles behave like waves, with a wavelength of h divided by momentum. A moving electron has a wavelength about the size of an atom, which is why electrons can show diffraction, and a baseball's is far too small to detect. electron_wavelength() gives it directly for an electron accelerated through a voltage (accurate up to about 10 kV).

The uncertainty principle says a particle cannot have a sharply defined position and momentum at the same time. Give one uncertainty and the function returns the smallest the other can be.

Formulas and arguments
CallFormulaWhat each argument means
matter.de_broglie(mass, speed, length)
  • mass () — mass of the particle, kg
  • speed () — its speed, m/s (well below the speed of light)
  • length — length unit for wavelengths and sizes you give and get back: "m" (default), "cm", "mm", "um", "nm", "pm", "angstrom", "fm"
Returns the wavelength
matter.de_broglie_energy(mass, energy, unit, length)
  • mass () — mass of the particle, kg
  • energy () — its kinetic energy
  • unit — energy unit for what you give and what you get back: "J" (default), "eV", "keV" or "MeV"
  • length — length unit for wavelengths and sizes you give and get back: "m" (default), "cm", "mm", "um", "nm", "pm", "angstrom", "fm"
Returns the wavelength
matter.electron_wavelength(volts, length)
  • volts () — accelerating voltage, V (up to about 10 kV)
  • length — length unit for wavelengths and sizes you give and get back: "m" (default), "cm", "mm", "um", "nm", "pm", "angstrom", "fm"
Returns the wavelength of an electron
matter.momentum_spread(position_spread, length)
  • position_spread () — the uncertainty in position
  • length — length unit for wavelengths and sizes you give and get back: "m" (default), "cm", "mm", "um", "nm", "pm", "angstrom", "fm"
Returns the smallest possible , in kg m/s ()
matter.position_spread(momentum_spread, length)
  • momentum_spread () — the uncertainty in momentum, kg m/s
  • length — length unit for wavelengths and sizes you give and get back: "m" (default), "cm", "mm", "um", "nm", "pm", "angstrom", "fm"
Returns the smallest possible
matter.speed_spread(mass, position_spread, length)
  • mass () — mass of the particle, kg
  • position_spread () — uncertainty in position
  • length — length unit for wavelengths and sizes you give and get back: "m" (default), "cm", "mm", "um", "nm", "pm", "angstrom", "fm"
Returns the smallest possible , in m/s
matter.energy_spread(lifetime)
  • lifetime () — how long the state lasts, s
Returns the smallest possible , in J
matter.de_broglie(constants.m_e, 1e6, length="nm") # 0.72739   (an electron at a thousand km/s)
matter.electron_wavelength(100, length="nm")  # 0.122643   (100 volts)
matter.de_broglie(0.145, 40)                  # 1.14243e-34   (a baseball, in metres)
matter.de_broglie_energy(constants.m_e, 100, unit="eV", length="nm") # 0.122643
matter.momentum_spread(1e-10)                 # 5.27286e-25   (an electron held to one angstrom)
matter.speed_spread(constants.m_e, 1e-10)     # 578838   (m/s)
matter.energy_spread(1e-9)                    # 5.27286e-26   (joules, for a state that lasts a nanosecond)

bohr.energy(n, Z=1) bohr.wavelength(n_from, n_to) bohr.spectrum(n_to) bohr.plot_levels()

The Bohr atom. Hydrogen, and any atom with just one electron (Z is the nuclear charge: 2 for He⁺, 3 for Li²⁺). The electron can only have certain energies, labelled by the whole number n, and the light an atom gives out or takes in comes from jumps between them. energy(), radius(), speed() and ionization() describe a level. transition() is the energy of the photon for a jump (positive when emitted), and wavelength() is its wavelength from the Rydberg formula. series_name() names the family of lines that end on a level: Lyman, Balmer, Paschen and so on.

spectrum(n_to) returns the first lines of a series as a table, and plot_levels() draws the energy levels with the Balmer jumps as arrows.

Formulas and arguments
CallFormulaWhat each argument means
bohr.energy(n, Z, unit)
  • n — the level, a whole number from 1
  • Z — nuclear charge, a whole number (default 1, hydrogen; 2 is He⁺)
  • unit — energy unit for what you give and what you get back: "J" (default), "eV", "keV" or "MeV"
Returns the energy of level (negative: bound)
bohr.radius(n, Z, length)
  • n — the level, a whole number from 1
  • Z — nuclear charge, a whole number (default 1, hydrogen; 2 is He⁺)
  • length — length unit for wavelengths and sizes you give and get back: "m" (default), "cm", "mm", "um", "nm", "pm", "angstrom", "fm"
Returns the orbit radius ( pm)
bohr.speed(n, Z)
  • n — the level, a whole number from 1
  • Z — nuclear charge, a whole number (default 1, hydrogen; 2 is He⁺)
Returns the electron's speed , in m/s ()
bohr.ionization(n, Z, unit)
  • n — the level it starts in (default 1)
  • Z — nuclear charge (default 1)
  • unit — energy unit for what you give and what you get back: "J" (default), "eV", "keV" or "MeV"
Returns the energy to free the electron
bohr.transition(n_from, n_to, Z, unit)
  • n_from () — the level it jumps from
  • n_to () — the level it jumps to
  • Z — nuclear charge
  • unit — energy unit for what you give and what you get back: "J" (default), "eV", "keV" or "MeV"
Returns the photon's energy (positive when given out, negative when taken in)
bohr.wavelength(n_from, n_to, Z, length)
  • n_from () — the level it jumps from
  • n_to () — the level it jumps to
  • Z — nuclear charge
  • length — length unit for wavelengths and sizes you give and get back: "m" (default), "cm", "mm", "um", "nm", "pm", "angstrom", "fm"
Returns the wavelength of the light (Rydberg formula, m⁻¹)
bohr.series_name(n_to)
  • n_to () — the level the jumps end on
Returns "Lyman", "Balmer", "Paschen", "Brackett", "Pfund"
bohr.spectrum(n_to, lines, Z, length, unit)
  • n_to () — the level the series ends on
  • lines — how many lines to list (default 5)
  • Z — nuclear charge
  • length — wavelength unit (default "nm")
  • unit — energy unit (default "eV")
Returns a table of n_from, n_to, wavelength and energy
bohr.plot_levels(Z, n_max, jumps, title, save)
  • Z — nuclear charge
  • n_max — highest level to draw (default 6)
  • jumps — list of (n_from, n_to) arrows to draw
  • title — heading
  • save — file name
Returns draws the energy levels
bohr.energy(1, unit="eV")                     # -13.6057
bohr.energy(2, unit="eV")                     # -3.40142
bohr.energy(1, Z=2, unit="eV")                # -54.4228   (He+: four times deeper)
bohr.ionization(1, unit="eV")                 # 13.6057
bohr.radius(1, length="pm")                   # 52.9177
bohr.speed(1)                                 # 2.18769e+06   (m/s)
bohr.wavelength(3, 2, length="nm")            # 656.112   (hydrogen alpha: red)
bohr.transition(3, 2, unit="eV")              # 1.88968
bohr.series_name(2)                           # 'Balmer'
bohr.spectrum(2, lines=4)   # the Balmer series
   n_from  n_to  wavelength_nm  energy_eV
0       3     2     656.112276   1.889680
1       4     2     486.009094   2.551067
2       5     2     433.936691   2.857196
3       6     2     410.070173   3.023487
bohr.plot_levels()
bohr.plot_levels(): the four visible Balmer lines are jumps down to n = 2
bohr.plot_levels(): the four visible Balmer lines are jumps down to n = 2

box.energy(n, mass, width) box.probability(n, width, x1, x2) box.photon_wavelength(n_from, n_to, mass, width) box.plot(width)

A particle in a box. A particle trapped between two walls can only have certain energies, and they grow as n². A narrower box means higher energies. levels() lists the lowest few, wavefunction() gives the wave function, and probability() gives the exact chance of finding the particle between two points. Notice that for n = 2 there is a node in the middle, where the chance of finding it is zero. photon_wavelength() gives the light for a jump between levels, and plot() draws the wave functions at the height of their energies. The width and the answers share the length= unit.

Formulas and arguments
CallFormulaWhat each argument means
box.energy(n, mass, width, unit, length)
  • n — the level, a whole number from 1
  • mass () — mass of the particle, kg
  • width () — width of the box (same length unit as length)
  • unit — energy unit for what you give and what you get back: "J" (default), "eV", "keV" or "MeV"
  • length — length unit for wavelengths and sizes you give and get back: "m" (default), "cm", "mm", "um", "nm", "pm", "angstrom", "fm"
Returns the energy of level
box.levels(mass, width, count, unit, length)
  • mass () — mass of the particle, kg
  • width () — width of the box (same length unit as length)
  • count — how many levels (default 4)
  • unit — energy unit for what you give and what you get back: "J" (default), "eV", "keV" or "MeV"
  • length — length unit for wavelengths and sizes you give and get back: "m" (default), "cm", "mm", "um", "nm", "pm", "angstrom", "fm"
Returns a list of the lowest levels
box.wavefunction(n, width, x, length)
  • n — the level
  • width () — width of the box
  • x — position, between 0 and
  • length — length unit for wavelengths and sizes you give and get back: "m" (default), "cm", "mm", "um", "nm", "pm", "angstrom", "fm"
Returns the value of the wave function
box.probability(n, width, x1, x2, length)
  • n — the level
  • width () — width of the box
  • x1 () — start of the region
  • x2 () — end of the region
  • length — length unit for wavelengths and sizes you give and get back: "m" (default), "cm", "mm", "um", "nm", "pm", "angstrom", "fm"
Returns the chance of finding the particle between and
box.transition(n_from, n_to, mass, width, unit, length)
  • n_from () — level it jumps from
  • n_to () — level it jumps to
  • mass () — mass of the particle, kg
  • width () — width of the box (same length unit as length)
  • unit — energy unit for what you give and what you get back: "J" (default), "eV", "keV" or "MeV"
  • length — length unit for wavelengths and sizes you give and get back: "m" (default), "cm", "mm", "um", "nm", "pm", "angstrom", "fm"
Returns the photon's energy
box.photon_wavelength(n_from, n_to, mass, width, length)
  • n_from () — level it jumps from
  • n_to () — level it jumps to
  • mass () — mass of the particle, kg
  • width () — width of the box (same length unit as length)
  • length — length unit for wavelengths and sizes you give and get back: "m" (default), "cm", "mm", "um", "nm", "pm", "angstrom", "fm" (the width and the answer share it)
Returns the wavelength of the light
box.plot(width, levels, mass, length, title, save)
  • width () — width of the box
  • levels — which levels to draw (default 1, 2, 3)
  • mass () — mass (default an electron)
  • length — unit of the width (default "nm")
Returns draws the wave functions at the height of their energies
box.energy(1, constants.m_e, 1, unit="eV", length="nm") # 0.37603   (an electron in a 1 nm box)
box.levels(constants.m_e, 1, 4, unit="eV", length="nm") # (0.37603, 1.50412, 3.38427, 6.01648)
box.energy(1, constants.m_e, 0.5, unit="eV", length="nm") # 1.50412   (half the width, four times the energy)
box.wavefunction(1, 1, 0.5)                   # 1.41421
box.probability(1, 1, 0, 0.5)                 # 0.5
box.probability(2, 1, 0.45, 0.55)             # 0.00645107   (around the node)
box.photon_wavelength(2, 1, constants.m_e, 1, length="nm") # 1099.06   (nanometres)
box.plot(1)
box.plot(1): each wave function sits at its energy, and n has n − 1 nodes
box.plot(1): each wave function sits at its energy, and n has n − 1 nodes

qho.energy(n, omega) qho.zero_point(omega) qho.spacing(omega)

The quantum harmonic oscillator. A particle on a spring has evenly spaced levels, n = 0, 1, 2, and so on, a gap of ħω between each. Even at the lowest level it has some energy left, the zero_point() energy, which is why nothing is ever perfectly still. omega() gives the angular frequency from a mass and spring constant, and turning_point() is how far the classical particle with that energy would swing.

Formulas and arguments
CallFormulaWhat each argument means
qho.omega(mass, k)
  • mass () — mass of the particle, kg
  • k — spring constant, N/m
Returns the angular frequency , rad/s
qho.energy(n, omega, unit)
  • n — the level, 0, 1, 2, …
  • omega () — angular frequency, rad/s
  • unit — energy unit for what you give and what you get back: "J" (default), "eV", "keV" or "MeV"
Returns the energy of level
qho.zero_point(omega, unit)
  • omega () — angular frequency, rad/s
  • unit — energy unit for what you give and what you get back: "J" (default), "eV", "keV" or "MeV"
Returns the lowest possible energy
qho.spacing(omega, unit)
  • omega () — angular frequency, rad/s
  • unit — energy unit for what you give and what you get back: "J" (default), "eV", "keV" or "MeV"
Returns the gap between neighbouring levels
qho.turning_point(n, mass, omega, length)
  • n — the level
  • mass () — mass, kg
  • omega () — angular frequency, rad/s
  • length — length unit for wavelengths and sizes you give and get back: "m" (default), "cm", "mm", "um", "nm", "pm", "angstrom", "fm"
Returns how far the classical particle with that energy would swing
qho.omega(constants.m_e, 1.0)                 # 1.04774e+15   (an electron on a 1 N/m spring)
qho.zero_point(1e15, unit="eV")               # 0.329106
qho.energy(0, 1e15, unit="eV")                # 0.329106
qho.energy(2, 1e15, unit="eV")                # 1.64553
qho.spacing(1e15, unit="eV")                  # 0.658212
qho.turning_point(2, constants.m_e, 1e15, length="nm") # 0.760814

tunnel.probability(mass, barrier, energy, width) tunnel.approx(...) tunnel.decay_constant(...)

Quantum tunnelling. A particle can get through a barrier taller than its energy, because its wave does not stop dead at the wall but dies away inside it. The thicker the barrier, or the greater the shortfall in energy, the smaller the chance. probability() is the exact chance for a rectangular barrier, approx() the simpler thick-barrier estimate e−2κL, and decay_constant() how quickly the wave fades. The energy must be less than the barrier height.

Formulas and arguments
CallFormulaWhat each argument means
tunnel.decay_constant(mass, barrier, energy, unit, length)
  • mass () — mass of the particle, kg
  • barrier () — height of the barrier
  • energy () — energy of the particle (less than ), both in the unit energy unit
  • unit — energy unit for what you give and what you get back: "J" (default), "eV", "keV" or "MeV"
  • length — the answer is per this length unit (default per metre)
Returns how fast the wave dies away in the barrier,
tunnel.approx(mass, barrier, energy, width, unit, length)
  • mass () — mass of the particle, kg
  • barrier () — height of the barrier
  • energy () — energy of the particle (less than ), both in the unit energy unit
  • width () — thickness of the barrier
  • unit — energy unit for what you give and what you get back: "J" (default), "eV", "keV" or "MeV"
  • length — length unit for wavelengths and sizes you give and get back: "m" (default), "cm", "mm", "um", "nm", "pm", "angstrom", "fm"
Returns the chance of getting through a thick barrier
tunnel.probability(mass, barrier, energy, width, unit, length)
  • mass () — mass of the particle, kg
  • barrier () — height of the barrier
  • energy () — energy of the particle (less than ), both in the unit energy unit
  • width () — thickness of the barrier
  • unit — energy unit for what you give and what you get back: "J" (default), "eV", "keV" or "MeV"
  • length — length unit for wavelengths and sizes you give and get back: "m" (default), "cm", "mm", "um", "nm", "pm", "angstrom", "fm"
Returns the exact chance of getting through
tunnel.probability(constants.m_e, 5, 1, 0.5, unit="eV", length="nm") # 9.08458e-05   (a 1 eV electron, 5 eV barrier, 0.5 nm)
tunnel.approx(constants.m_e, 5, 1, 0.5, unit="eV", length="nm") # 3.54873e-05
tunnel.probability(constants.m_e, 5, 1, 1.0, unit="eV", length="nm") # 3.22394e-09   (twice as thick)
tunnel.probability(constants.m_e, 5, 4, 0.5, unit="eV", length="nm") # 0.015199   (nearly over the top)
tunnel.decay_constant(constants.m_e, 5, 1, unit="eV", length="nm") # 10.2463   (per nm)

nuclear.binding_energy(Z, A, atomic_mass) nuclear.remaining(N0, t, half_life) nuclear.age(fraction_left, half_life)

The nucleus. radius() gives the size of a nucleus from its mass number. A nucleus is lighter than the particles it is made of, and the missing mass is the energy that binds it: mass_defect() (in atomic mass units), binding_energy() and binding_per_nucleon(), which peaks near iron. Give the atom's atomic mass in u. mass_energy() is Einstein's E = mc², and u_to_energy() converts atomic mass units to energy (1 u is 931.494 MeV).

Radioactive decay. Each nucleus has a fixed chance per unit time of decaying, so the number left falls by half every half-life. remaining() gives the number left, activity() the decays per unit time, age() the age of a sample from the fraction of the isotope still present (carbon dating), and decay_law solves N = N₀e−λt for any missing variable. Use the same time unit for the time and the half-life. plot() draws the decay curve.

Formulas and arguments
CallFormulaWhat each argument means
nuclear.radius(mass_number, length)
  • mass_number () — the number of protons plus neutrons
  • length — length unit for wavelengths and sizes you give and get back: "m" (default), "cm", "mm", "um", "nm", "pm", "angstrom", "fm"
Returns the nuclear radius ( fm)
nuclear.mass_energy(mass, unit)
  • mass () — mass, kg
  • unit — energy unit for what you give and what you get back: "J" (default), "eV", "keV" or "MeV"
Returns the energy locked up in that mass
nuclear.u_to_energy(mass_u, unit)
  • mass_u () — mass in atomic mass units
  • unit — energy unit (default "MeV")
Returns the equivalent energy
nuclear.mass_defect(Z, A, atomic_mass)
  • Z — number of protons
  • A — mass number (protons plus neutrons)
  • atomic_mass () — mass of the atom, in u
Returns the missing mass , in u
nuclear.binding_energy(Z, A, atomic_mass, unit)
  • Z — number of protons
  • A — mass number
  • atomic_mass () — atomic mass, u
  • unit — energy unit (default "MeV")
Returns the binding energy
nuclear.binding_per_nucleon(Z, A, atomic_mass, unit)
  • Z — protons
  • A — mass number
  • atomic_mass () — atomic mass, u
  • unit — energy unit (default "MeV")
Returns the binding energy for each particle
nuclear.decay_constant(half_life)
  • half_life () — time for half to decay (in any time unit)
Returns the decay constant , per that time unit
nuclear.half_life(decay_constant)
  • decay_constant () — the decay constant
Returns the half-life
nuclear.mean_life(half_life)
  • half_life () — the half-life
Returns the mean life
nuclear.remaining(N0, t, half_life)
  • N0 () — the number at the start
  • t — time passed
  • half_life () — the half-life (the same time unit as )
Returns the number left
nuclear.decay_law(N=, N0=, lam=, t=)
  • N — number left
  • N0 () — number at the start
  • lam () — decay constant
  • t — time passed
Returns the variable you leave out
nuclear.activity(N, half_life)
  • N — the number of nuclei present
  • half_life () — the half-life
Returns decays per unit time (becquerels if the half-life is in seconds)
nuclear.age(fraction_left, half_life)
  • fraction_left () — fraction of the isotope still present, 0 to 1
  • half_life () — the half-life
Returns the age (carbon dating)
nuclear.time_to_fraction(fraction_left, half_life)
  • fraction_left () — the fraction left, 0 to 1
  • half_life () — the half-life
Returns the time until that fraction is left
nuclear.path(N0, half_life, half_lives, steps) / plot(...)
  • N0 () — starting number
  • half_life () — half-life
  • half_lives — how many half-lives to cover (default 6)
  • steps — points
  • time_label — label for the time axis
Returns a table of t, N; or a chart
nuclear.radius(27, length="fm")               # 3.6   (aluminium)
nuclear.mass_defect(2, 4, 4.002603254)        # 0.0303766   (helium-4, in u)
nuclear.binding_energy(2, 4, 4.002603254)     # 28.2957   (MeV)
nuclear.binding_per_nucleon(2, 4, 4.002603254) # 7.07392
nuclear.binding_per_nucleon(26, 56, 55.934936) # 8.79035   (iron-56 holds tighter)
nuclear.mass_energy(1e-3)                     # 8.98755e+13   (joules in one gram)
nuclear.u_to_energy(1)                        # 931.494
nuclear.decay_constant(5730)                  # 0.000120968   (carbon-14, per year)
nuclear.remaining(1000, 11460, 5730)          # 250   (two half-lives)
nuclear.age(0.25, 5730)                       # 11460   (years)
nuclear.activity(1e6, 3600)                   # 192.541   (per second, half-life of an hour)
nuclear.decay_law(N0=1000, lam=0.1, t=5)      # 606.531
nuclear.decay_law(N=500, N0=1000, lam=0.1)    # 6.93147   (the time to halve)
nuclear.plot(1000, 5730, half_lives=5, time_label="time (years)")
nuclear.plot(1000, 5730, half_lives=5): each half-life halves what is left
nuclear.plot(1000, 5730, half_lives=5): each half-life halves what is left

xray.bragg_angle(spacing, wavelength) xray.bragg_spacing(angle, wavelength) xray.min_wavelength(volts)

X-rays. A crystal reflects X-rays strongly only at certain angles, set by the spacing of its planes of atoms (Bragg's law), which is how crystal structures are measured. bragg_angle() and bragg_spacing() convert between the two, with an optional order=. If no reflection is possible at all, you are told. min_wavelength() and max_energy() give the shortest wavelength and the biggest photon energy an X-ray tube can make at a given voltage.

Formulas and arguments
CallFormulaWhat each argument means
xray.bragg_angle(spacing, wavelength, order, length, mode)
  • spacing () — distance between the crystal planes
  • wavelength () — wavelength of the X-rays
  • order () — order of the reflection (default 1)
  • length — length unit for wavelengths and sizes you give and get back: "m" (default), "cm", "mm", "um", "nm", "pm", "angstrom", "fm"
  • mode — "deg" (default) or "rad" for the angle
Returns the glancing angle
xray.bragg_spacing(angle, wavelength, order, length, mode)
  • angle () — the glancing angle (degrees)
  • wavelength () — wavelength
  • order () — order (default 1)
  • length — length unit for wavelengths and sizes you give and get back: "m" (default), "cm", "mm", "um", "nm", "pm", "angstrom", "fm"
Returns the plane spacing
xray.min_wavelength(volts, length)
  • volts () — voltage across the X-ray tube, V
  • length — length unit for wavelengths and sizes you give and get back: "m" (default), "cm", "mm", "um", "nm", "pm", "angstrom", "fm"
Returns the shortest wavelength the tube makes
xray.max_energy(volts, unit)
  • volts () — voltage across the tube, V
  • unit — energy unit for what you give and what you get back: "J" (default), "eV", "keV" or "MeV"
Returns the biggest photon energy
xray.bragg_angle(0.282, 0.154, length="nm")   # 15.8458   (degrees, rock salt with copper X-rays)
xray.bragg_angle(0.282, 0.154, order=2, length="nm") # 33.0998
xray.bragg_spacing(15.8, 0.154, length="nm")  # 0.282797
xray.min_wavelength(30000, length="pm")       # 41.3281
xray.max_energy(30000, unit="keV")            # 30

Thermodynamics & Statistical Physics

Gases, heat and its transfer, entropy and the laws of thermodynamics, engines, and the statistics of many particles. Every function is shown with its formula, and each table says which argument is which letter of it.

gas.pressure(1, 300, 0.0248) gas.law(P=, V=, n=, T=) engine.carnot(600, 300)

How the thermodynamics functions work. Everything is in SI units: pressure in pascals, volume in cubic metres, temperature in kelvin, energy in joules, amount of gas in moles, and molar mass in kg/mol. Use heat.to_kelvin(25, "C") to change a thermometer reading, or convert() for other units. The statistical functions accept unit= ("J", "eV", "keV", "MeV") for energies.

Solve for any missing variable. The main laws take their variables by name: give all but one and the law finds the last. gas.law(P=, V=, n=, T=), gas.combined(...), heat.law(Q=, m=, c=, dT=) and engine.first_law(dU=, Q=, W=) work this way, and explain=True prints the equation, what you gave, and the result. Every function below is shown with its formula, and the table says which argument is which letter of the formula.

gas.pressure(1, 300, 0.0248)                  # 100578   (1 mol at 300 K in 24.8 litres, in pascals)
gas.law(P=101325, V=0.0248, T=300)            # 1.00743   (solves for n, in mol)
gas.law(P=101325, V=0.0248, T=300, explain=True) # 1.00743
heat.to_kelvin(25)                            # 298.15   (25 degrees C)
constants.R                                   # 8.31446

gas.law(P=, V=, n=, T=) gas.pressure(n, T, V) gas.volume(n, T, P) gas.moles(P, V, T) gas.temperature(P, V, n) gas.combined(...) gas.density(...) gas.van_der_waals(...)

The gas law. An ideal gas obeys P V = n R T. pressure(), volume(), moles() and temperature() each solve it for one quantity, and law() solves for whichever one you leave out. combined() compares two states of the same gas (P1 V1 / T1 = P2 V2 / T2) and finds any one of the six. density() gives the density of a gas of a given molar mass. van_der_waals() corrects the pressure for real gases: a is the pull between molecules and b the volume they fill themselves.

Formulas and arguments
CallFormulaWhat each argument means
gas.law(P=, V=, n=, T=)
  • P — pressure, Pa
  • V — volume, m³
  • n — amount of gas, mol
  • T — temperature, kelvin (give three, it finds the fourth)
Returns the missing one of , , , ( J/(mol K))
gas.pressure(n, T, V)
  • n — amount of gas, mol
  • T — temperature, kelvin
  • V — volume, m³
Returns the pressure , in Pa
gas.volume(n, T, P)
  • n — amount of gas, mol
  • T — temperature, kelvin
  • P — pressure, Pa
Returns the volume , in m³
gas.moles(P, V, T)
  • P — pressure, Pa
  • V — volume, m³
  • T — temperature, kelvin
Returns the amount of gas , in mol
gas.temperature(P, V, n)
  • P — pressure, Pa
  • V — volume, m³
  • n — amount of gas, mol
Returns the temperature , in K
gas.combined(P1=, V1=, T1=, P2=, V2=, T2=)
  • P1 () — starting pressure
  • V1 () — starting volume
  • T1 () — starting temperature, kelvin
  • P2 () — final pressure
  • V2 () — final volume
  • T2 () — final temperature, kelvin (give five, it finds the sixth)
Returns the missing one of the six
gas.density(P, T, molar_mass)
  • P — pressure, Pa
  • T — temperature, kelvin
  • molar_mass () — molar mass of the gas, kg/mol (nitrogen is 0.028)
Returns the density , in kg/m³
gas.van_der_waals(n, V, T, a, b)
  • n — amount of gas, mol
  • V — volume, m³
  • T — temperature, kelvin
  • a — strength of the attraction between molecules, Pa m⁶/mol²
  • b — volume of a mole of molecules themselves, m³/mol
Returns the pressure of a real gas, in Pa
gas.volume(1, 273.15, 101325)                 # 0.022414   (molar volume at 0 °C and 1 atm, m³)
gas.moles(101325, 0.0224, 273.15)             # 0.999377
gas.temperature(101325, 0.0248, 1)            # 302.228
gas.combined(P1=1e5, V1=1, T1=300, P2=2e5, T2=600) # 1   (V2 in m³)
gas.density(101325, 300, 0.029)               # 1.17804   (air, kg/m³)
gas.pressure(1, 300, 0.0248)                  # 100578   (ideal)
gas.van_der_waals(1, 0.0248, 300, 0.1378, 3.18e-5) # 100483   (CO₂, real)

gas.isothermal_work(n, T, V1, V2) gas.isobaric_work(P, V1, V2) gas.adiabatic_pressure(...) gas.adiabatic_temperature(...) gas.adiabatic_work(...)

Work done by a gas. A gas that expands pushes on its surroundings and does work. isothermal_work() is for a change at constant temperature, isobaric_work() for constant pressure. In an adiabatic change no heat flows in or out, so squeezing the gas heats it: adiabatic_pressure() and adiabatic_temperature() give the end state, and adiabatic_work() the work. The work is positive when the gas does work on the surroundings, negative when it is squeezed. gamma defaults to 5/3 for a single-atom gas; use 1.4 for air.

Formulas and arguments
CallFormulaWhat each argument means
gas.isothermal_work(n, T, V1, V2)
  • n — amount of gas, mol
  • T — temperature, kelvin, held constant
  • V1 () — starting volume
  • V2 () — final volume
Returns the work done BY the gas, in J (negative when squeezed)
gas.isobaric_work(P, V1, V2)
  • P — the constant pressure, Pa
  • V1 () — starting volume
  • V2 () — final volume
Returns the work done BY the gas, in J
gas.adiabatic_pressure(P1, V1, V2, gamma)
  • P1 () — starting pressure
  • V1 () — starting volume
  • V2 () — final volume
  • gamma () — ratio (default 5/3 for a single-atom gas; 1.4 for air)
Returns the final pressure , in Pa
gas.adiabatic_temperature(T1, V1, V2, gamma)
  • T1 () — starting temperature, kelvin
  • V1 () — starting volume
  • V2 () — final volume
  • gamma () — ratio (default 5/3 for a single-atom gas; 1.4 for air)
Returns the final temperature , in K
gas.adiabatic_work(P1, V1, P2, V2, gamma)
  • P1 () — starting pressure
  • V1 () — starting volume
  • P2 () — final pressure
  • V2 () — final volume
  • gamma () — ratio (default 5/3 for a single-atom gas; 1.4 for air)
Returns the work done BY the gas, in J
gas.isothermal_work(1, 300, 1, 2)             # 1728.94   (doubling the volume of 1 mol at 300 K, J)
gas.isobaric_work(1e5, 1e-3, 3e-3)            # 200
gas.adiabatic_pressure(1e5, 1e-3, 5e-4, 1.4)  # 263902   (air squeezed to half the volume)
gas.adiabatic_temperature(300, 1e-3, 5e-4, 1.4) # 395.852   (and heated, K)
gas.adiabatic_work(1e5, 1e-3, 2.639e5, 5e-4, 1.4) # -79.875   (negative: work was done on the gas)

gas.cv(dof) gas.cp(dof) gas.gamma(dof) gas.internal_energy(n, T, dof) gas.kinetic_energy(T)

Heat capacity and energy. A molecule can store energy in each of its "degrees of freedom": 3 for a single atom (it can only move in three directions), 5 for a diatomic gas such as nitrogen (add two rotations), 6 for a bent or larger molecule. The equipartition rule gives each one ½ k T. From that, cv(), cp(), gamma() and internal_energy() follow. kinetic_energy() is the average energy of one molecule's motion, 3/2 k T, in whatever energy unit you choose.

Formulas and arguments
CallFormulaWhat each argument means
gas.cv(dof)
  • dof () — degrees of freedom per molecule: 3 for a single atom (default), 5 for a diatomic gas, 6 for a polyatomic gas
Returns the molar heat capacity at constant volume, in J/(mol K)
gas.cp(dof)
  • dof () — degrees of freedom per molecule: 3 for a single atom (default), 5 for a diatomic gas, 6 for a polyatomic gas
Returns the molar heat capacity at constant pressure, in J/(mol K)
gas.gamma(dof)
  • dof () — degrees of freedom per molecule: 3 for a single atom (default), 5 for a diatomic gas, 6 for a polyatomic gas
Returns the ratio (no unit)
gas.internal_energy(n, T, dof)
  • n — amount of gas, mol
  • T — temperature, kelvin
  • dof () — degrees of freedom per molecule: 3 for a single atom (default), 5 for a diatomic gas, 6 for a polyatomic gas
Returns the energy stored in the gas, in J
gas.kinetic_energy(T, unit)
  • T — temperature, kelvin
  • unit — energy unit for what you give and what you get back: "J" (default), "eV", "keV" or "MeV"
Returns the average translational energy of one molecule
gas.cv()                                      # 12.4717   (single atom, J/(mol K))
gas.cp(5)                                     # 29.1006   (diatomic)
gas.gamma(5)                                  # 1.4
gas.internal_energy(1, 300, 3)                # 3741.51   (J)
gas.kinetic_energy(300, unit="eV")            # 0.038778   (one molecule at room temperature)

gas.rms_speed(T, M) gas.mean_speed(T, M) gas.most_probable_speed(T, M) gas.maxwell(v, T, M) gas.speed_fraction(v1, v2, T, M) gas.mean_free_path(T, P, d) gas.plot_maxwell(M, temperatures)

How fast the molecules move. The molecules of a gas move at a spread of speeds, described by the Maxwell distribution. most_probable_speed() is the peak of the curve, mean_speed() the average, and rms_speed() the root-mean-square, which sets the energy. All three grow as the square root of the temperature. maxwell() gives the height of the curve, and speed_fraction() the share of molecules between two speeds. mean_free_path() is how far a molecule goes between collisions. plot_maxwell() draws the curves for several temperatures. Molar mass is in kg/mol (nitrogen: 0.028).

Formulas and arguments
CallFormulaWhat each argument means
gas.rms_speed(T, molar_mass)
  • T — temperature, kelvin
  • molar_mass () — molar mass of the gas, kg/mol (nitrogen is 0.028)
Returns the root-mean-square speed, in m/s
gas.mean_speed(T, molar_mass)
  • T — temperature, kelvin
  • molar_mass () — molar mass of the gas, kg/mol (nitrogen is 0.028)
Returns the average speed, in m/s
gas.most_probable_speed(T, molar_mass)
  • T — temperature, kelvin
  • molar_mass () — molar mass of the gas, kg/mol (nitrogen is 0.028)
Returns the speed most molecules have, in m/s
gas.maxwell(v, T, molar_mass)
  • v — the speed you ask about, m/s
  • T — temperature, kelvin
  • molar_mass () — molar mass of the gas, kg/mol (nitrogen is 0.028)
Returns the share of molecules per unit of speed, , in s/m
gas.speed_fraction(v1, v2, T, molar_mass)
  • v1 () — lower speed, m/s
  • v2 () — upper speed, m/s
  • T — temperature, kelvin
  • molar_mass () — molar mass of the gas, kg/mol (nitrogen is 0.028)
Returns the fraction of molecules with speeds between and (0 to 1)
gas.mean_free_path(T, P, diameter)
  • T — temperature, kelvin
  • P — pressure, Pa
  • diameter () — diameter of a molecule, m
Returns the average distance between collisions , in m
gas.plot_maxwell(molar_mass, temperatures, title, save)
  • molar_mass () — molar mass of the gas, kg/mol (nitrogen is 0.028)
  • temperatures () — list of temperatures, kelvin (default 100, 300, 1000)
  • title — chart heading
  • save — file name to write the chart to
Returns draws the Maxwell curves with the most probable speeds marked
gas.rms_speed(300, 0.028)                     # 516.963   (nitrogen at room temperature, m/s)
gas.mean_speed(300, 0.028)                    # 476.287
gas.most_probable_speed(300, 0.028)           # 422.098
gas.maxwell(400, 300, 0.028)                  # 0.00195593
gas.speed_fraction(0, 422.1, 300, 0.028)      # 0.427596   (slower than the peak)
gas.mean_free_path(300, 101325, 3.7e-10)      # 6.72078e-08   (metres: about 67 nm)
gas.plot_maxwell(0.028, temperatures=(100, 300, 1000))
gas.plot_maxwell(): hotter gas is faster and more spread out; the dots mark the most probable speed
gas.plot_maxwell(): hotter gas is faster and more spread out; the dots mark the most probable speed

heat.energy(m, c, dT) heat.law(Q=, m=, c=, dT=) heat.latent(m, L) heat.mix(...) heat.to_kelvin(value, scale) heat.linear_expansion(...) heat.newton_cooling(...)

Heat and temperature. energy() is the heat needed to change the temperature of a mass, and law() solves the same equation for any of its four variables. latent() is the heat to melt or boil with no change of temperature. mix() gives the final temperature when two bodies are mixed (use the same scale for both and the answer comes back in it). to_kelvin() and from_kelvin() change between °C, °F and K. linear_expansion() and volume_expansion() give how much a heated body grows, and newton_cooling() the temperature of something cooling in a room.

Formulas and arguments
CallFormulaWhat each argument means
heat.energy(m, c, dT)
  • m — mass, kg
  • c — specific heat capacity, J/(kg K) (water is 4186)
  • dT () — temperature change, K
Returns the heat needed, in J
heat.law(Q=, m=, c=, dT=)
  • Q — heat added, J
  • m — mass, kg
  • c — specific heat capacity, J/(kg K)
  • dT () — temperature change, K (give three, it finds the fourth)
Returns the missing one of the four
heat.latent(m, L)
  • m — mass that melts or boils, kg
  • L — latent heat, J/kg (ice melting is 334000)
Returns the heat needed, in J
heat.mix(m1, c1, T1, m2, c2, T2)
  • m1 () — mass of the first body, kg
  • c1 () — its specific heat capacity
  • T1 () — its temperature
  • m2 () — mass of the second body
  • c2 () — its specific heat capacity
  • T2 () — its temperature (use one scale for both, say °C)
Returns the final temperature , in the scale you used
heat.to_kelvin(value, scale)
  • value — the temperature you have
  • scale — "C" (default), "F" or "K"
Returns the temperature in kelvin
heat.from_kelvin(kelvin, scale)
  • kelvin () — the temperature in kelvin
  • scale — "C" (default), "F" or "K"
Returns the temperature in that scale
heat.linear_expansion(L0, alpha, dT)
  • L0 () — original length
  • alpha () — coefficient of linear expansion, 1/K
  • dT () — temperature change, K
Returns the change in length (same unit as )
heat.volume_expansion(V0, beta, dT)
  • V0 () — original volume
  • beta () — coefficient of volume expansion, 1/K (about 3 α for a solid)
  • dT () — temperature change, K
Returns the change in volume
heat.newton_cooling(T0, T_env, k, t)
  • T0 () — starting temperature
  • T_env () — temperature of the surroundings
  • k — cooling rate constant, 1/s
  • t — time that has passed, s
Returns the temperature after time
heat.energy(2, 4186, 30)                      # 251160   (2 kg of water warmed by 30 K, J)
heat.law(Q=251160, m=2, c=4186)               # 30   (solves for the temperature change)
heat.latent(1, 334000)                        # 334000   (melting 1 kg of ice, J)
heat.mix(1, 4186, 20, 0.5, 900, 100)          # 27.7653   (water and hot aluminium, °C)
heat.to_kelvin(98.6, "F")                     # 310.15
heat.from_kelvin(300, "F")                    # 80.33
heat.linear_expansion(10, 12e-6, 50)          # 0.006   (a 10 m steel rod heated by 50 K, m)
heat.newton_cooling(90, 20, 0.05, 10)         # 62.4571   (a hot drink after 10 s, °C)

heat.conduction(k, A, dT, L) heat.resistance(L, k, A) heat.convection(h, A, dT) heat.radiation(emissivity, area, T, T_surroundings)

How heat travels. Heat moves by conduction through solids, by convection into a moving fluid, and by radiation as light. conduction() gives the power through a slab, and resistance() its thermal resistance (resistances of layers in series add). convection() is Newton's law of cooling for a surface in a fluid. radiation() is the net power a surface gives out: what it radiates minus what it takes in from the surroundings.

Formulas and arguments
CallFormulaWhat each argument means
heat.conduction(k, A, dT, L)
  • k — thermal conductivity, W/(m K)
  • A — area of the slab, m²
  • dT () — temperature difference across it, K
  • L — thickness, m
Returns the rate of heat flow , in W
heat.resistance(L, k, A)
  • L — thickness, m
  • k — thermal conductivity, W/(m K)
  • A — area, m²
Returns the thermal resistance, in K/W
heat.convection(h, A, dT)
  • h — convection coefficient, W/(m² K)
  • A — surface area, m²
  • dT () — temperature difference between surface and fluid, K
Returns the rate of heat flow , in W
heat.radiation(emissivity, area, T, T_surroundings)
  • emissivity () — 0 to 1 (1 is a perfect black body)
  • area () — surface area, m²
  • T — temperature of the body, kelvin
  • T_surroundings () — temperature of the surroundings, kelvin (default 0)
Returns the net power radiated, in W ( W/(m² K⁴))
heat.conduction(0.8, 10, 20, 0.2)             # 800   (a 10 m² brick wall, 20 K difference, W)
heat.resistance(0.2, 0.8, 10)                 # 0.025   (K/W)
heat.convection(10, 2, 15)                    # 300   (W)
heat.radiation(1, 1, 300)                     # 459.3   (a black body at 300 K, W)
heat.radiation(0.9, 1.5, 310, 293)            # 142.778   (a person in a room, W)

entropy.heat(Q, T) entropy.phase_change(m, L, T) entropy.heating(m, c, T1, T2) entropy.ideal_gas(...) entropy.universe(Q, T_hot, T_cold) entropy.gibbs(H, T, S)

Entropy. Entropy measures how spread out the energy is, and the total of it never falls in a real process. heat() is the entropy gained when heat flows in at a steady temperature, phase_change() the same for melting or boiling, and heating() when the temperature itself changes. ideal_gas() gives the change for a gas that changes temperature and volume. universe() shows why heat flows from hot to cold: the total entropy change is positive in that direction, and Fyrefly refuses the reverse. helmholtz() and gibbs() give the free energies that tell you whether a change will happen by itself.

Formulas and arguments
CallFormulaWhat each argument means
entropy.heat(Q, T)
  • Q — heat flowing in, J (negative if flowing out)
  • T — the constant temperature, kelvin
Returns the entropy change , in J/K
entropy.phase_change(m, L, T)
  • m — mass that melts or boils, kg
  • L — latent heat, J/kg
  • T — temperature of the change, kelvin
Returns the entropy change , in J/K
entropy.heating(m, c, T1, T2)
  • m — mass, kg
  • c — specific heat capacity, J/(kg K)
  • T1 () — starting temperature, kelvin
  • T2 () — final temperature, kelvin
Returns the entropy change , in J/K
entropy.ideal_gas(n, T1, T2, V1, V2, dof)
  • n — amount of gas, mol
  • T1 () — starting temperature
  • T2 () — final temperature
  • V1 () — starting volume
  • V2 () — final volume
  • dof () — degrees of freedom per molecule: 3 for a single atom (default), 5 for a diatomic gas, 6 for a polyatomic gas
Returns the entropy change of the gas, in J/K
entropy.universe(Q, T_hot, T_cold)
  • Q — heat that flows, J
  • T_hot () — temperature of the hot body, kelvin
  • T_cold () — temperature of the cold body, kelvin
Returns the total entropy change, which is never negative
entropy.helmholtz(U, T, S)
  • U — internal energy, J
  • T — temperature, kelvin
  • S — entropy, J/K
Returns the Helmholtz free energy , in J
entropy.gibbs(H, T, S)
  • H — enthalpy, J
  • T — temperature, kelvin
  • S — entropy, J/K
Returns the Gibbs free energy , in J (a change happens by itself when falls)
entropy.heat(1000, 300)                       # 3.33333   (J/K)
entropy.phase_change(1, 334000, 273.15)       # 1222.77   (melting 1 kg of ice)
entropy.heating(1, 4186, 293.15, 373.15)      # 1010.07   (heating 1 kg of water from 20 to 100 °C)
entropy.ideal_gas(1, 300, 300, 1, 2)          # 5.76315   (doubling the volume at constant temperature)
entropy.universe(1000, 500, 300)              # 1.33333   (1000 J flowing from 500 K to 300 K)
entropy.gibbs(-92200, 298, -198.7)            # -32987.4   (ammonia synthesis, J)

engine.first_law(dU=, Q=, W=) engine.efficiency(work, heat_in) engine.carnot(T_hot, T_cold) engine.refrigerator_cop(...) engine.heat_pump_cop(...) engine.otto(r) engine.diesel(r, rc) engine.plot_carnot()

The laws of thermodynamics and engines. The first law says energy is conserved: the energy in a gas rises by the heat put in minus the work it does (first_law() solves for any one of the three). The second law limits engines: no engine can be more efficient than a Carnot engine, whose efficiency carnot() depends only on the two temperatures (in kelvin). refrigerator_cop() and heat_pump_cop() give the best "coefficient of performance" for cooling and heating. max_work() and exhaust() split the heat taken in into work and waste. otto() and diesel() are the efficiencies of the ideal petrol and diesel cycles. plot_carnot() draws the four steps of the Carnot cycle on a pressure-volume diagram.

Formulas and arguments
CallFormulaWhat each argument means
engine.first_law(dU=, Q=, W=)
  • dU () — change in internal energy, J
  • Q — heat added to the gas, J
  • W — work done BY the gas, J (give two, it finds the third)
Returns the missing one of , ,
engine.efficiency(work, heat_in)
  • work () — useful work out, J
  • heat_in () — heat taken in, J
Returns the efficiency (0 to 1)
engine.carnot(T_hot, T_cold)
  • T_hot () — temperature of the hot side, kelvin
  • T_cold () — temperature of the cold side, kelvin
Returns the greatest efficiency any engine can have between these temperatures
engine.refrigerator_cop(T_hot, T_cold)
  • T_hot () — temperature of the hot side, kelvin
  • T_cold () — temperature of the cold side, kelvin
Returns the best heat removed per unit of work
engine.heat_pump_cop(T_hot, T_cold)
  • T_hot () — temperature of the hot side, kelvin
  • T_cold () — temperature of the cold side, kelvin
Returns the best heat delivered per unit of work
engine.max_work(Q_hot, T_hot, T_cold)
  • Q_hot () — heat taken from the hot side, J
  • T_hot () — temperature of the hot side, kelvin
  • T_cold () — temperature of the cold side, kelvin
Returns the most work possible, in J
engine.exhaust(Q_hot, work)
  • Q_hot () — heat taken in, J
  • work () — work given out, J
Returns the heat thrown away to the cold side, in J
engine.otto(compression_ratio, gamma)
  • compression_ratio () — volume before the squeeze divided by volume after (greater than 1)
  • gamma () — ratio of the gas (default 1.4, air)
Returns the efficiency of the ideal petrol-engine cycle
engine.diesel(compression_ratio, cutoff_ratio, gamma)
  • compression_ratio () — compression ratio (greater than 1)
  • cutoff_ratio () — volume at the end of burning divided by volume at the start (between 1 and r)
  • gamma () — ratio (default 1.4)
Returns the efficiency of the ideal diesel cycle
engine.plot_carnot(T_hot, T_cold, V1, V2, n, gamma, title, save)
  • T_hot () — hot temperature, K (default 600)
  • T_cold () — cold temperature, K (default 300)
  • V1 () — smallest volume, m³
  • V2 () — volume at the end of the hot expansion, m³
  • n — moles of gas
  • gamma () — heat capacity ratio
  • title — chart heading
  • save — file name to write the chart to
Returns draws the four steps of the Carnot cycle on a pressure-volume diagram
engine.first_law(Q=500, W=200)                # 300   (dU in J)
engine.first_law(dU=300, W=200)               # 500   (Q in J)
engine.efficiency(300, 1000)                  # 0.3
engine.carnot(600, 300)                       # 0.5
engine.refrigerator_cop(300, 270)             # 9
engine.heat_pump_cop(300, 270)                # 10
engine.max_work(1000, 600, 300)               # 500   (J)
engine.exhaust(1000, 400)                     # 600
engine.otto(10)                               # 0.601893   (compression ratio 10)
engine.diesel(18, 2)                          # 0.631578
engine.plot_carnot(600, 300)
engine.plot_carnot(600, 300): two isotherms and two adiabats; the area inside is the work done each cycle
engine.plot_carnot(600, 300): two isotherms and two adiabats; the area inside is the work done each cycle

stat.boltzmann_factor(E, T) stat.population_ratio(E_upper, E_lower, T) stat.partition(energies, T) stat.probabilities(energies, T) stat.average_energy(energies, T) stat.equipartition(T, dof)

Boltzmann statistics. A system at temperature T is likelier to be found in a low-energy state than a high one, by the Boltzmann factor e^(-E / kT). population_ratio() gives the ratio of two levels (with their degeneracies, the number of states that share an energy), such as the atoms in the n = 2 and n = 1 levels of hydrogen in a star. partition() adds up the factors for a list of levels, and probabilities() and average_energy() follow from it. equipartition() is the average energy of one particle with dof quadratic terms. Energies take unit=.

Formulas and arguments
CallFormulaWhat each argument means
stat.boltzmann_factor(energy, T, unit)
  • energy () — energy of the state
  • T — temperature, kelvin
  • unit — energy unit for what you give and what you get back: "J" (default), "eV", "keV" or "MeV"
Returns the Boltzmann factor: how likely the state is compared with a zero-energy state
stat.population_ratio(E_upper, E_lower, T, g_upper, g_lower, unit)
  • E_upper () — energy of the upper level
  • E_lower () — energy of the lower level
  • T — temperature, kelvin
  • g_upper () — number of states in the upper level (default 1)
  • g_lower () — number of states in the lower level (default 1)
  • unit — energy unit for what you give and what you get back: "J" (default), "eV", "keV" or "MeV"
Returns the ratio of the populations
stat.partition(energies, T, degeneracies, unit)
  • energies () — list of the energies of the levels (measured from any zero)
  • T — temperature, kelvin
  • degeneracies () — how many states share each energy (default 1 each)
  • unit — energy unit for what you give and what you get back: "J" (default), "eV", "keV" or "MeV"
Returns the partition function (energies measured from the lowest level)
stat.probabilities(energies, T, degeneracies, unit)
  • energies () — list of the energies of the levels (measured from any zero)
  • T — temperature, kelvin
  • degeneracies () — how many states share each energy (default 1 each)
  • unit — energy unit for what you give and what you get back: "J" (default), "eV", "keV" or "MeV"
Returns a list: the chance of finding the system in each level
stat.average_energy(energies, T, degeneracies, unit)
  • energies () — list of the energies of the levels (measured from any zero)
  • T — temperature, kelvin
  • degeneracies () — how many states share each energy (default 1 each)
  • unit — energy unit for what you give and what you get back: "J" (default), "eV", "keV" or "MeV"
Returns the average energy over the levels
stat.equipartition(T, dof, unit)
  • T — temperature, kelvin
  • dof () — number of quadratic energy terms (3 for the motion of a free particle)
  • unit — energy unit for what you give and what you get back: "J" (default), "eV", "keV" or "MeV"
Returns the average energy of one particle
stat.boltzmann_factor(0.025, 300, unit="eV")  # 0.380206   (about one kT)
stat.population_ratio(10.2, 0, 5800, unit="eV", g_upper=8, g_lower=2) # 5.48279e-09   (hydrogen n=2 to n=1 in the Sun)
stat.partition([0, 1e-21], 300)               # 1.7855   (a two-level system, energies in J)
stat.probabilities([0, 1e-21], 300)           # (0.560067, 0.439933)
stat.average_energy([0, 1e-21], 300)          # 4.39933e-22
stat.probabilities([0, 0], 300, degeneracies=[1, 3]) # (0.25, 0.75)
stat.equipartition(300, 3, unit="eV")         # 0.038778

stat.entropy(W) stat.gibbs_entropy(probabilities) stat.microstates(N, n) stat.einstein_solid(N, q)

Counting microstates. Boltzmann's great idea is that entropy counts the number of ways to arrange things: S = k ln W. microstates() counts the ways to choose n of N (heads in a run of coin flips), and einstein_solid() counts the ways for q energy units to be shared among N oscillators. entropy() turns a count into an entropy, and gibbs_entropy() does it from a list of probabilities. Counts are returned as exact whole numbers, however large.

Formulas and arguments
CallFormulaWhat each argument means
stat.entropy(W)
  • W — the number of microstates (ways to arrange things)
Returns the entropy , in J/K
stat.gibbs_entropy(probabilities)
  • probabilities () — list of the chance of each state (adding to 1)
Returns the entropy , in J/K
stat.microstates(N, n)
  • N — the number of particles or trials
  • n — how many are chosen (heads, excited atoms)
Returns the number of ways
stat.einstein_solid(N, q)
  • N — the number of oscillators
  • q — the number of energy units shared among them
Returns the number of microstates
stat.microstates(10, 5)                       # 252   (ways to get 5 heads in 10 flips)
stat.microstates(100, 50)                     # 100891344545564193334812497256
stat.einstein_solid(4, 3)                     # 20
stat.entropy(2**100)                          # 9.56993e-22   (J/K)
stat.gibbs_entropy([0.5, 0.5])                # 9.56993e-24   (J/K, one coin: k ln 2)

stat.fermi_dirac(E, mu, T) stat.bose_einstein(E, mu, T) stat.fermi_energy(n) stat.fermi_temperature(n) stat.plot_occupancy()

Quantum statistics. Identical quantum particles are not independent. Fermions (electrons, protons) obey the Pauli rule that no two can share a state, so the chance a state is filled, fermi_dirac(), never exceeds 1. Bosons (photons) can pile into one state, so bose_einstein() can be larger than 1. Both approach the classical Boltzmann curve when the states are sparsely filled. mu is the chemical potential, which sets where the filling changes over. fermi_energy() is the energy of the highest filled level in a metal at absolute zero, and fermi_temperature() writes it as a temperature, which shows why the electrons in a metal behave quantum mechanically even at room temperature.

Formulas and arguments
CallFormulaWhat each argument means
stat.fermi_dirac(energy, mu, T, unit)
  • energy () — energy of the state
  • mu () — chemical potential (the Fermi energy for electrons in a metal)
  • T — temperature, kelvin
  • unit — energy unit for what you give and what you get back: "J" (default), "eV", "keV" or "MeV"
Returns the chance the state is filled, between 0 and 1
stat.bose_einstein(energy, mu, T, unit)
  • energy () — energy of the state (above mu)
  • mu () — chemical potential
  • T — temperature, kelvin
  • unit — energy unit for what you give and what you get back: "J" (default), "eV", "keV" or "MeV"
Returns the average number of bosons in the state
stat.fermi_energy(density, mass, unit)
  • density () — free electrons per m³ (copper: 8.49×10²⁸)
  • mass () — particle mass, kg (default: an electron)
  • unit — energy unit for what you give and what you get back: "J" (default), "eV", "keV" or "MeV"
Returns the Fermi energy
stat.fermi_temperature(density, mass)
  • density () — free electrons per m³
  • mass () — particle mass, kg (default: an electron)
Returns the Fermi temperature , in K
stat.plot_occupancy(mu, T, unit, title, save)
  • mu () — chemical potential (default 0)
  • T — temperature, kelvin
  • unit — energy unit for the axis
  • title — chart heading
  • save — file name to write the chart to
Returns draws the Fermi-Dirac, Bose-Einstein, and Boltzmann curves
stat.fermi_dirac(0, 0, 300)                   # 0.5   (at mu: exactly half)
stat.fermi_dirac(-0.1, 0, 300, unit="eV")     # 0.979531   (0.1 eV below mu: nearly full)
stat.fermi_dirac(0.1, 0, 300, unit="eV")      # 0.0204688
stat.bose_einstein(0.1, 0, 300, unit="eV")    # 0.0213425
stat.fermi_energy(8.49e28, unit="eV")         # 7.04383   (copper)
stat.fermi_temperature(8.49e28)               # 81740.2   (kelvin)
stat.plot_occupancy(0, 300, unit="eV")
stat.plot_occupancy(): fermions fill up to 1, bosons pile up, and both match Boltzmann when the states are sparse
stat.plot_occupancy(): fermions fill up to 1, bosons pile up, and both match Boltzmann when the states are sparse

stat.einstein_heat_capacity(T, theta) stat.debye_heat_capacity(T, theta)

The heat capacity of solids. Classically every solid should store 3R (about 24.9 J/(mol K)) per mole at any temperature, and at room temperature most do. But at low temperature quantum effects freeze out the vibrations and the heat capacity drops towards zero. Einstein's model of identical vibrating atoms explains the drop, and Debye's model of waves in the solid explains it better. theta is the characteristic temperature of the solid (copper's Debye temperature is 343 K).

Formulas and arguments
CallFormulaWhat each argument means
stat.einstein_heat_capacity(T, theta)
  • T — temperature, kelvin
  • theta () — Einstein temperature of the solid, kelvin
Returns the molar heat capacity, in J/(mol K) (approaches when hot)
stat.debye_heat_capacity(T, theta)
  • T — temperature, kelvin
  • theta () — Debye temperature of the solid, kelvin (copper: 343)
Returns the molar heat capacity, in J/(mol K)
stat.debye_heat_capacity(300, 343)            # 23.3862   (copper at room temperature, J/(mol K))
stat.debye_heat_capacity(20, 343)             # 0.385289   (copper at 20 K)
stat.einstein_heat_capacity(300, 240)         # 23.6546
stat.einstein_heat_capacity(1e6, 100)         # 24.9434   (very hot: tends to 3R)
3 * constants.R                               # 24.9434

Fluid Mechanics

Flow in pipes, drag and lift, surface tension, liquids at rest, and rivers. Every function is shown with its formula, and each table says which argument is which letter of it.

pipe.poiseuille(radius, length, pressure_drop, viscosity) drag.terminal_speed(mass, Cd, rho, area) surface.capillary_rise(surface_tension, radius)

Fluid mechanics, beyond the basics. The fluid namespace (in Classical Mechanics) covers pressure, buoyancy, Bernoulli's equation and the Reynolds number. These five namespaces go further: pipe for flow in pipes, drag for the force of air and water on moving bodies, surface for surface tension, hydro for liquids at rest, and channel for rivers and weirs. Everything is in SI units: metres, kilograms, seconds and pascals. Water is the default liquid (1000 kg/m³), viscosity is the dynamic viscosity in Pa s (water about 0.001, air 0.000018), and the drag and sound functions default to air.

Every function is shown with its formula, and each table says which argument is which letter of the formula, and what comes back. Wrong inputs, such as a negative radius or a wide section smaller than the narrow one, give a plain sentence saying what to fix.

pipe.poiseuille(0.01, 2, 5000, 0.001)         # 0.00981748   (water through a 2 m pipe of radius 1 cm, m³/s)
drag.terminal_speed(80, 1.0, 1.2, 0.7)        # 43.227   (a skydiver, m/s)
surface.capillary_rise(0.0728, 0.0005)        # 0.029684   (water in a tube of radius 0.5 mm, m)
fluid.reynolds(1000, 1, 0.05, 0.001)          # 50000

pipe.poiseuille(radius, length, pressure_drop, viscosity) pipe.pressure_drop(...) pipe.resistance(...) pipe.regime(reynolds)

Smooth flow in a pipe. At a low Reynolds number the flow is smooth, or laminar, and Poiseuille's law gives the flow rate: it grows as the fourth power of the radius, so doubling the radius of a pipe multiplies the flow by sixteen. pressure_drop() turns it round to give the pressure needed for a given flow, and resistance() gives the hydraulic resistance (so that ΔP = R Q). regime() says whether a Reynolds number (from fluid.reynolds()) means laminar, transitional or turbulent flow.

Formulas and arguments
CallFormulaWhat each argument means
pipe.poiseuille(radius, length, pressure_drop, viscosity)
  • radius () — inner radius of the pipe, m
  • length () — length of the pipe, m
  • pressure_drop () — pressure difference between the ends, Pa
  • viscosity () — dynamic viscosity, Pa s (water is about 0.001)
Returns the flow rate , in m³/s (valid for smooth, laminar flow)
pipe.pressure_drop(radius, length, flow_rate, viscosity)
  • radius () — inner radius, m
  • length () — length, m
  • flow_rate () — flow rate, m³/s
  • viscosity () — dynamic viscosity, Pa s (water is about 0.001)
Returns the pressure difference needed, in Pa
pipe.resistance(radius, length, viscosity)
  • radius () — inner radius, m
  • length () — length, m
  • viscosity () — dynamic viscosity, Pa s (water is about 0.001)
Returns the hydraulic resistance , in Pa s/m³ (so that )
pipe.regime(reynolds)
  • reynolds () — the Reynolds number from fluid.reynolds()
Returns the word "laminar", "transitional" or "turbulent"
pipe.poiseuille(0.01, 2, 5000, 0.001)         # 0.00981748   (m³/s)
pipe.poiseuille(0.02, 2, 5000, 0.001)         # 0.15708   (radius doubled: 16 times as much)
pipe.pressure_drop(0.01, 2, 0.0098175, 0.001) # 5000.01   (Pa)
pipe.resistance(0.01, 2, 0.001)               # 509296   (Pa s/m³)
fluid.reynolds(1000, 1, 0.05, 0.001)          # 50000
pipe.regime(1000)                             # 'laminar'
pipe.regime(3000)                             # 'transitional'
pipe.regime(50000)                            # 'turbulent'

pipe.friction_factor(reynolds, diameter, roughness) pipe.darcy_weisbach(f, length, diameter, rho, v) pipe.head_loss(...) pipe.pump_power(...) pipe.plot_friction()

Friction and pumping. When the flow is turbulent, the pressure lost along a pipe follows the Darcy–Weisbach equation, with a friction factor f that depends on the Reynolds number and on the roughness of the pipe wall. friction_factor() gives 64/Re for laminar flow and the Swamee–Jain approximation of the Colebrook formula for turbulent flow. darcy_weisbach() gives the pressure lost, head_loss() the same loss written as a height of liquid, and pump_power() the power a pump must supply. plot_friction() draws the friction factor against the Reynolds number for several roughnesses (a Moody chart). In the transitional range (Re 2300 to 4000) the flow is unpredictable and the turbulent formula is only a guide.

Formulas and arguments
CallFormulaWhat each argument means
pipe.friction_factor(reynolds, diameter, roughness)
  • reynolds () — the Reynolds number
  • diameter () — inner diameter, m
  • roughness () — height of the bumps on the pipe wall, m (default 0, smooth)
Returns the Darcy friction factor (no unit)
pipe.darcy_weisbach(f, length, diameter, rho, v)
  • f — friction factor from friction_factor()
  • length () — length of the pipe, m
  • diameter () — inner diameter, m
  • rho () — density of the fluid, kg/m³
  • v — average speed of the flow, m/s
Returns the pressure lost to friction , in Pa
pipe.head_loss(f, length, diameter, v, g)
  • f — friction factor
  • length () — length of the pipe, m
  • diameter () — inner diameter, m
  • v — average speed, m/s
  • g — gravitational field strength, m/s² (default 9.81)
Returns the friction loss as a height of liquid , in m
pipe.pump_power(flow_rate, pressure_rise, efficiency)
  • flow_rate () — flow rate, m³/s
  • pressure_rise () — pressure the pump adds, Pa
  • efficiency () — pump efficiency, 0 to 1 (default 1)
Returns the power the pump needs, in W
pipe.plot_friction(diameter, roughness, title, save)
  • diameter () — pipe diameter, m (default 0.05)
  • roughness () — list of wall roughnesses to draw, m
  • title — chart heading
  • save — file name to write the chart to
Returns draws friction factor against Reynolds number (a Moody chart)
pipe.friction_factor(1000, 0.05)              # 0.064   (laminar: 64 / Re)
pipe.friction_factor(1e5, 0.05)               # 0.0178626   (smooth pipe, turbulent)
pipe.friction_factor(1e5, 0.05, 1e-4)         # 0.0253329   (wall bumps 0.1 mm high)
pipe.darcy_weisbach(0.02, 100, 0.05, 1000, 2) # 80000   (Pa lost along 100 m of 5 cm pipe at 2 m/s)
pipe.head_loss(0.02, 100, 0.05, 2)            # 8.15494   (m of water)
pipe.pump_power(0.01, 200000, 0.7)            # 2857.14   (W for 10 litres per second against 200 kPa)
pipe.plot_friction(0.05)
pipe.plot_friction(): a Moody chart. Laminar flow follows 64/Re, and rough pipes lose more in turbulent flow
pipe.plot_friction(): a Moody chart. Laminar flow follows 64/Re, and rough pipes lose more in turbulent flow

pipe.venturi(A1, A2, pressure_drop, rho) pipe.pitot(pressure_difference, rho) pipe.orifice(Cd, area, pressure_drop, rho)

Measuring flow. These meters all work by Bernoulli's principle: faster flow means lower pressure, so a pressure difference tells you the speed. A venturi meter is a pipe that narrows, and venturi() gives the flow rate from the pressure fall between the wide and narrow parts. A Pitot tube faces into a flow and pitot() gives its speed (the default fluid is air, as on an aircraft). orifice() gives the flow through a hole or orifice plate; Cd is the discharge coefficient, about 0.6 for a sharp-edged hole.

Formulas and arguments
CallFormulaWhat each argument means
pipe.venturi(A1, A2, pressure_drop, rho)
  • A1 () — area of the wide part, m²
  • A2 () — area of the narrow part, m²
  • pressure_drop () — pressure fall from wide to narrow part, Pa
  • rho () — density of the liquid, kg/m³ (default 1000, water)
Returns the flow rate , in m³/s
pipe.pitot(pressure_difference, rho)
  • pressure_difference () — stagnation minus static pressure, Pa
  • rho () — density of the fluid, kg/m³ (default 1.2, air)
Returns the speed of the flow, in m/s
pipe.orifice(Cd, area, pressure_drop, rho)
  • Cd () — discharge coefficient, about 0.6 for a sharp hole
  • area () — area of the hole, m²
  • pressure_drop () — pressure difference across it, Pa
  • rho () — density of the liquid, kg/m³ (default 1000, water)
Returns the flow rate , in m³/s
pipe.venturi(0.01, 0.005, 20000)              # 0.0365148   (m³/s)
pipe.pitot(600)                               # 31.6228   (an aircraft's airspeed from 600 Pa, m/s)
pipe.orifice(0.62, 0.001, 50000)              # 0.0062   (m³/s)

drag.force(Cd, rho, area, v) drag.lift(Cl, rho, area, v) drag.power(...) drag.dynamic_pressure(rho, v) drag.stokes(radius, speed, viscosity) drag.stall_speed(...)

Drag and lift. At everyday speeds, the push of air or water on a moving body grows as the square of the speed: force() is ½ ρ Cd A v², and the power needed to overcome it grows as the cube, which is why high speeds cost so much fuel. lift() is the same formula with a lift coefficient, and stall_speed() is the slowest speed at which a wing can hold up an aircraft's weight. dynamic_pressure() is ½ρv². For a small sphere moving slowly through a thick fluid, stokes() gives the drag. Typical drag coefficients: 0.3 for a car, 0.47 for a sphere, about 1 for a skydiver.

Formulas and arguments
CallFormulaWhat each argument means
drag.dynamic_pressure(rho, v)
  • rho () — density of the fluid, kg/m³ (air is 1.2)
  • v — speed of the flow, m/s
Returns the dynamic pressure , in Pa
drag.force(Cd, rho, area, v)
  • Cd () — drag coefficient (about 1 for a skydiver, 0.3 for a car, 0.47 for a sphere)
  • rho () — density of the fluid, kg/m³ (air is 1.2)
  • area () — frontal area facing the flow, m²
  • v — speed relative to the fluid, m/s
Returns the drag force , in N
drag.lift(Cl, rho, area, v)
  • Cl () — lift coefficient (about 1.5 for a wing at take-off)
  • rho () — density of the fluid, kg/m³ (air is 1.2)
  • area () — wing area, m²
  • v — speed relative to the air, m/s
Returns the lift force , in N
drag.power(Cd, rho, area, v)
  • Cd () — drag coefficient (about 1 for a skydiver, 0.3 for a car, 0.47 for a sphere)
  • rho () — density of the fluid, kg/m³ (air is 1.2)
  • area () — frontal area facing the flow, m²
  • v — speed, m/s
Returns the power needed to overcome drag, in W
drag.stokes(radius, speed, viscosity)
  • radius () — radius of the sphere, m
  • speed () — its (slow) speed, m/s
  • viscosity () — dynamic viscosity, Pa s (water is about 0.001)
Returns the drag force , in N
drag.stall_speed(mass, Cl_max, rho, area, g)
  • mass () — mass of the aircraft, kg
  • Cl_max () — the greatest lift coefficient of the wing
  • rho () — density of the fluid, kg/m³ (air is 1.2)
  • area () — wing area, m²
  • g — gravitational field strength, m/s² (default 9.81)
Returns the slowest speed that keeps it flying, in m/s
drag.dynamic_pressure(1.2, 30)                # 540   (Pa)
drag.force(0.3, 1.2, 2.2, 30)                 # 356.4   (a car at 30 m/s (108 km/h), N)
drag.power(0.3, 1.2, 2.2, 30)                 # 10692   (W)
drag.power(0.3, 1.2, 2.2, 60)                 # 85536   (double the speed: 8 times the power)
drag.lift(1.2, 1.2, 16, 60)                   # 41472   (N)
drag.stall_speed(1100, 1.6, 1.225, 16)        # 26.2336   (a light aircraft, m/s)
drag.stokes(1e-3, 0.1, 1.0)                   # 0.00188496   (a 1 mm bead in syrup, N)

drag.terminal_speed(mass, Cd, rho, area) drag.fall_speed(..., t) drag.fall_distance(..., t) drag.plot_fall()

Falling through air. A falling body speeds up until the drag equals its weight, and after that it falls at a steady terminal speed. terminal_speed() gives it, and fall_speed() and fall_distance() give the speed and the distance fallen after any time, from the exact solution (which approaches the terminal speed as tanh). plot_fall() draws speed against time, beside the straight line you would get with no air. For a basic version with Stokes' law, see fluid.terminal_velocity().

Formulas and arguments
CallFormulaWhat each argument means
drag.terminal_speed(mass, Cd, rho, area, g)
  • mass () — mass of the falling body, kg
  • Cd () — drag coefficient (about 1 for a skydiver, 0.3 for a car, 0.47 for a sphere)
  • rho () — density of the fluid, kg/m³ (air is 1.2)
  • area () — frontal area facing the flow, m²
  • g — gravitational field strength, m/s² (default 9.81)
Returns the terminal speed , where drag equals weight, in m/s
drag.fall_speed(mass, Cd, rho, area, t, g)
  • mass () — mass, kg
  • Cd () — drag coefficient
  • rho () — density of the fluid
  • area () — frontal area, m²
  • t — time since release, s
  • g — gravitational field strength, m/s² (default 9.81)
Returns the speed after time , in m/s
drag.fall_distance(mass, Cd, rho, area, t, g)
  • mass () — mass, kg
  • Cd () — drag coefficient
  • rho () — density of the fluid
  • area () — frontal area, m²
  • t — time since release, s
  • g — gravitational field strength, m/s² (default 9.81)
Returns the distance fallen , in m
drag.plot_fall(mass, Cd, rho, area, duration, title, save)
  • mass () — mass, kg (default 80)
  • Cd () — drag coefficient
  • rho () — density of the air
  • area () — frontal area, m²
  • duration — time to draw, s
  • title — chart heading
  • save — file name to write the chart to
Returns draws speed against time, beside the fall with no air
drag.terminal_speed(80, 1.0, 1.2, 0.7)        # 43.227   (a skydiver, spread out, m/s)
drag.terminal_speed(80, 0.7, 1.2, 0.18)       # 101.887   (head down, m/s)
drag.fall_speed(80, 1, 1.2, 0.7, 5)           # 35.1273   (speed after 5 s, m/s)
drag.fall_distance(80, 1, 1.2, 0.7, 5)        # 102.844   (fallen in 5 s, m (122.6 with no air))
drag.fall_speed(80, 1, 1.2, 0.7, 60)          # 43.227   (long after: terminal speed)
drag.plot_fall(80, 1.0, 1.2, 0.7)
drag.plot_fall(): speed levels off at the terminal speed instead of growing without limit
drag.plot_fall(): speed levels off at the terminal speed instead of growing without limit

drag.sound_speed(T, gamma, molar_mass) drag.mach(v, T)

The speed of sound. Sound travels faster in warmer gas, and slower in heavier gas. sound_speed() gives it for any gas from its temperature, its heat-capacity ratio gamma, and its molar mass; the defaults are for air. mach() gives a speed as a multiple of the speed of sound: below 1 is subsonic, 1 is the speed of sound, above 1 is supersonic.

Formulas and arguments
CallFormulaWhat each argument means
drag.sound_speed(T, gamma, molar_mass)
  • T — temperature, kelvin (default 293.15)
  • gamma () — ratio of heat capacities (default 1.4, air)
  • molar_mass () — molar mass of the gas, kg/mol (default 0.0289647, air)
Returns the speed of sound , in m/s
drag.mach(v, T, gamma, molar_mass)
  • v — speed of the body, m/s
  • T — temperature of the air, kelvin
  • gamma () — ratio of heat capacities
  • molar_mass () — molar mass of the gas, kg/mol
Returns the Mach number (no unit): 1 is the speed of sound
drag.sound_speed()                            # 343.235   (air at 20 °C, m/s)
drag.sound_speed(273.15)                      # 331.32   (air at 0 °C)
drag.sound_speed(293.15, 5/3, 0.004)          # 1007.76   (helium)
drag.mach(300)                                # 0.874037
drag.mach(340.3)                              # 0.991449

surface.capillary_rise(surface_tension, radius) surface.laplace(surface_tension, radius, surfaces) surface.energy(...) surface.young(...) surface.weber(...)

Surface tension. The surface of a liquid acts like a stretched skin, with a tension γ in N/m (water is 0.0728). It pulls water up a thin tube, capillary_rise(), higher the thinner the tube; liquids that do not wet the wall, such as mercury, are pushed down instead, and the answer is negative. It squeezes a drop or bubble, so the pressure inside is higher by laplace(): 2γ/r for a drop, 4γ/r for a soap bubble with two surfaces. energy() is the energy stored in a surface, film_force() the pull of a soap film on a wire, and capillary_length() the size below which surface tension beats gravity. young() gives the contact angle of a drop on a surface from the three surface tensions, and weber() compares inertia with surface tension.

Formulas and arguments
CallFormulaWhat each argument means
surface.capillary_rise(surface_tension, radius, rho, contact_angle, g)
  • surface_tension () — surface tension, N/m (water is 0.0728)
  • radius () — radius of the tube, m
  • rho () — density of the liquid, kg/m³
  • contact_angle () — angle the liquid makes with the wall, degrees (default 0)
  • g — gravitational field strength, m/s² (default 9.81)
Returns the height the liquid climbs, in m (negative if it is pushed down, as mercury is)
surface.laplace(surface_tension, radius, surfaces)
  • surface_tension () — surface tension, N/m
  • radius () — radius of the drop or bubble, m
  • surfaces — 1 for a drop (default), 2 for a soap bubble
Returns the extra pressure inside, in Pa
surface.energy(surface_tension, area)
  • surface_tension () — surface tension, N/m
  • area () — area of the surface, m²
Returns the energy stored in the surface, in J
surface.film_force(surface_tension, length, surfaces)
  • surface_tension () — surface tension, N/m
  • length () — length of the wire the film pulls on, m
  • surfaces () — number of surfaces (default 2, a soap film)
Returns the pull of the film, in N
surface.capillary_length(surface_tension, rho, g)
  • surface_tension () — surface tension, N/m
  • rho () — density of the liquid, kg/m³ (default 1000, water)
  • g — gravitational field strength, m/s² (default 9.81)
Returns the length below which surface tension beats gravity, in m
surface.young(gamma_solid_vapour, gamma_solid_liquid, gamma_liquid_vapour)
  • gamma_solid_vapour () — surface tension between solid and vapour
  • gamma_solid_liquid () — between solid and liquid
  • gamma_liquid_vapour () — between liquid and vapour (all in N/m)
Returns the contact angle of a drop on the surface, in degrees
surface.weber(rho, v, length, surface_tension)
  • rho () — density of the liquid, kg/m³
  • v — speed, m/s
  • length () — size of the drop or jet, m
  • surface_tension () — surface tension, N/m
Returns the Weber number (no unit): large means inertia beats surface tension
surface.capillary_rise(0.0728, 0.0005)        # 0.029684   (water in a 1 mm wide tube, m)
surface.capillary_rise(0.4, 0.0005, 13600, 140) # -0.00918684   (mercury in glass: pushed down, m)
surface.laplace(0.0728, 1e-3)                 # 145.6   (a 1 mm water drop, Pa)
surface.laplace(0.025, 0.01, 2)               # 10   (a soap bubble of radius 1 cm, Pa)
surface.energy(0.0728, 0.01)                  # 0.000728   (J in 100 cm² of water surface)
surface.film_force(0.025, 0.1)                # 0.005   (N)
surface.capillary_length(0.0728)              # 0.00272415   (m)
surface.young(0.05, 0.02, 0.0728)             # 65.6639   (degrees)
surface.weber(1000, 1, 0.001, 0.0728)         # 13.7363

hydro.dam_force(width, depth) hydro.center_of_pressure(depth) hydro.dam_moment(width, depth) hydro.gate_force(width, height, top_depth) hydro.manometer(height)

Liquid against walls. The pressure of a liquid grows with depth, so it pushes harder at the bottom of a wall than at the top. dam_force() is the total push on a wall holding back a depth H: ½ ρ g w H², which grows as the square of the depth. The push acts two-thirds of the way down (center_of_pressure()), and dam_moment() is its turning effect about the foot of the wall. gate_force() is the push on a rectangular gate below the surface. manometer() gives the pressure difference shown by a column of liquid (the default is mercury).

Formulas and arguments
CallFormulaWhat each argument means
hydro.dam_force(width, depth, rho, g)
  • width () — width of the wall, m
  • depth () — depth of the liquid against it, m
  • rho () — density of the liquid, kg/m³ (default 1000, water)
  • g — gravitational field strength, m/s² (default 9.81)
Returns the total push of the liquid on the wall, in N
hydro.center_of_pressure(depth)
  • depth () — depth of the liquid against the wall, m
Returns how far below the surface the push acts, in m
hydro.dam_moment(width, depth, rho, g)
  • width () — width of the wall, m
  • depth () — depth of the liquid, m
  • rho () — density of the liquid, kg/m³ (default 1000, water)
  • g — gravitational field strength, m/s² (default 9.81)
Returns the turning effect about the foot of the wall, in N m
hydro.gate_force(width, height, top_depth, rho, g)
  • width () — width of the gate, m
  • height () — height of the gate, m
  • top_depth () — depth of its top edge below the surface, m (default 0)
  • rho () — density of the liquid, kg/m³ (default 1000, water)
  • g — gravitational field strength, m/s² (default 9.81)
Returns the push on the gate, in N
hydro.manometer(height, rho, g)
  • height () — difference in height of the liquid columns, m
  • rho () — density of the manometer liquid, kg/m³ (default 13600, mercury)
  • g — gravitational field strength, m/s² (default 9.81)
Returns the pressure difference , in Pa
hydro.dam_force(100, 20)                      # 1.962e+08   (a 100 m wide dam holding 20 m of water, N)
hydro.center_of_pressure(20)                  # 13.3333   (m below the surface)
hydro.dam_moment(100, 20)                     # 1.308e+09   (N m)
hydro.dam_force(100, 40) / hydro.dam_force(100, 20) # 4   (twice as deep: 4 times the push)
hydro.gate_force(2, 1, 3)                     # 68670   (a 2 m by 1 m gate with its top 3 m down, N)
hydro.manometer(0.76)                         # 101396   (760 mm of mercury, Pa)

hydro.barometric(height) hydro.spinning_surface(omega, r) hydro.drain_time(tank_area, hole_area, height)

The atmosphere, spinning liquids, and draining tanks. Air pressure falls with height: barometric() gives it for air at a steady temperature, halving about every 5.5 km. A liquid in a spinning bucket climbs the wall, and spinning_surface() gives how much higher the surface is at a distance r from the axis than at the centre (it is a parabola). drain_time() is how long a tank takes to empty through a hole at its base, from Torricelli's speed √(2gh).

Formulas and arguments
CallFormulaWhat each argument means
hydro.barometric(height, P0, T, molar_mass, g)
  • height () — height above the ground, m
  • P0 () — pressure at the ground, Pa (default 101325)
  • T — air temperature, kelvin (default 288.15)
  • molar_mass () — molar mass of air, kg/mol (default 0.0289647)
  • g — gravitational field strength, m/s² (default 9.81)
Returns the air pressure at that height, in Pa
hydro.spinning_surface(omega, r, g)
  • omega () — angular speed of the bucket, rad/s
  • r — distance from the axis, m
  • g — gravitational field strength, m/s² (default 9.81)
Returns how much higher the liquid surface is at than at the centre, in m
hydro.drain_time(tank_area, hole_area, height, Cd, g)
  • tank_area () — cross-section of the tank, m²
  • hole_area () — area of the hole at the base, m²
  • height () — starting depth of the liquid, m
  • Cd () — discharge coefficient (default 1)
  • g — gravitational field strength, m/s² (default 9.81)
Returns the time for the tank to empty, in s
hydro.barometric(0)                           # 101325   (sea level, Pa)
hydro.barometric(5000)                        # 55998.9   (at 5 km)
hydro.barometric(8849)                        # 35475.4   (the top of Everest)
hydro.spinning_surface(10, 0.1)               # 0.0509684   (10 rad/s, 10 cm from the axis, m)
hydro.drain_time(1, 0.001, 2, 0.62)           # 1029.92   (a 1 m² tank, 2 m deep, 10 cm² hole, s)

channel.manning(n, hydraulic_radius, slope) channel.froude(v, depth) channel.wave_speed(depth) channel.regime(froude) channel.critical_depth(q) channel.weir(Cd, width, head)

Rivers, canals and weirs. Water in an open channel runs downhill, and manning() gives its speed from the roughness n of the bed, the hydraulic radius (cross-section area over wetted perimeter) and the slope. The Froude number compares the flow speed with the speed of a wave on the water, wave_speed(): below 1 the flow is subcritical (waves can travel upstream), above 1 it is supercritical (they cannot). critical_depth() is the depth at which the flow is exactly critical, and weir() gives the flow over a rectangular weir from the height of the water above its crest.

Formulas and arguments
CallFormulaWhat each argument means
channel.manning(n, hydraulic_radius, slope)
  • n — Manning roughness of the channel (0.012 smooth concrete, 0.03 earth, 0.05 weedy)
  • hydraulic_radius () — cross-section area divided by wetted perimeter, m
  • slope () — drop in height per length of channel (no unit)
Returns the speed of the flow, in m/s
channel.froude(v, depth, g)
  • v — speed of the flow, m/s
  • depth () — depth of the water, m
  • g — gravitational field strength, m/s² (default 9.81)
Returns the Froude number (no unit)
channel.wave_speed(depth, g)
  • depth () — depth of the water, m
  • g — gravitational field strength, m/s² (default 9.81)
Returns the speed of a shallow-water wave, in m/s
channel.regime(froude)
  • froude () — the Froude number
Returns the word "subcritical", "critical" or "supercritical"
channel.critical_depth(flow_per_width, g)
  • flow_per_width () — flow in m³/s for each metre of channel width
  • g — gravitational field strength, m/s² (default 9.81)
Returns the depth at which the flow is critical, in m
channel.weir(Cd, width, head, g)
  • Cd () — discharge coefficient (about 0.6)
  • width () — width of the weir, m
  • head () — height of the water surface above the weir crest, m
  • g — gravitational field strength, m/s² (default 9.81)
Returns the flow rate over the weir, in m³/s
channel.manning(0.03, 1.5, 0.001)             # 1.38125   (an earth channel, m/s)
channel.froude(2, 1)                          # 0.638551   (2 m/s in 1 m of water)
channel.wave_speed(2)                         # 4.42945   (m/s)
channel.regime(0.64)                          # 'subcritical'
channel.regime(2)                             # 'supercritical'
channel.critical_depth(3)                     # 0.971683   (3 m³/s per metre of width, m)
channel.weir(0.62, 2, 0.3)                    # 0.601675   (m³/s)

Electricity & Magnetism

Charges and fields, capacitors, circuits, magnetic forces, induction, alternating current, and electromagnetic waves. Every function is shown with its formula, and each table says which argument is which letter of it.

charge.coulomb(q1, q2, r) circuit.ohm(V=, I=, R=) magnet.cyclotron_radius(m, v, q, B) ac.resonance(L, C)

How the electricity and magnetism functions work. Everything is in SI units: coulombs, volts, amperes, ohms, farads, henries, teslas, hertz, metres and seconds, and angles in degrees (add mode="rad" for radians). There are seven namespaces: charge (forces, fields and potentials), capacitor, circuit (DC circuits, RC and RL), magnet (magnetic forces and fields), induction (Faraday's law, generators, transformers, inductors), ac (alternating current and resonance), and emwave (light and other electromagnetic waves). The constants are in constants: k_e, eps0, mu0, e, m_e, m_p and c.

Solve for any missing variable. circuit.ohm(V=, I=, R=), circuit.power_law(...) and capacitor.law(...) take their variables by name: give all but one and the law finds the last. Add explain=True to print the equation, what you gave, and the result. Every function below is shown with its formula, and each table says which argument is which letter of the formula, and what comes back.

charge.coulomb(1e-6, -2e-6, 0.1)              # -1.79751   (N: negative means they attract)
circuit.ohm(V=12, R=4)                        # 3   (amperes)
circuit.ohm(V=12, R=4, explain=True)          # 3
magnet.cyclotron_radius(constants.m_e, 1e7, constants.e, 0.01) # 0.00568563   (metres)
constants.k_e                                 # 8.98755e+09
constants.mu0                                 # 1.25664e-06

charge.coulomb(q1, q2, r) charge.field(q, r) charge.potential(q, r) charge.potential_energy(...) charge.force_on(q, E) charge.acceleration(q, m, E)

Charges, forces and fields. Two charges push or pull each other with a force that grows with the product of the charges and falls with the square of the distance (Coulomb's law). coulomb() returns a positive number when the charges repel and a negative one when they attract. A charge also fills the space round it with an electric field (field(), force per unit charge) and a potential (potential(), energy per unit charge). potential_energy() is the energy stored by a pair. force_on() and acceleration() give what a field does to a charge. Charges are in coulombs: a proton has constants.e = 1.6×10⁻¹⁹ C.

Formulas and arguments
CallFormulaWhat each argument means
charge.coulomb(q1, q2, r)
  • q1 () — first charge, C
  • q2 () — second charge, C
  • r — distance between them, m
Returns the force , in N: positive pushes them apart, negative pulls them together ( N m²/C²)
charge.field(q, r)
  • q — the charge making the field, C
  • r — distance from it, m
Returns the electric field , in N/C (pointing away if is positive)
charge.potential(q, r)
  • q — the charge, C
  • r — distance from it, m
Returns the electric potential , in volts
charge.potential_energy(q1, q2, r)
  • q1 () — first charge, C
  • q2 () — second charge, C
  • r — distance between them, m
Returns the energy stored by the pair, in J (negative for opposite charges)
charge.force_on(q, E)
  • q — the charge, C
  • E — the electric field it is in, N/C
Returns the force on the charge, in N
charge.acceleration(q, m, E)
  • q — the charge, C
  • m — its mass, kg
  • E — the electric field, N/C
Returns the acceleration , in m/s²
charge.coulomb(1e-6, -2e-6, 0.1)              # -1.79751   (N, attraction)
charge.coulomb(1e-6, 2e-6, 0.2) / charge.coulomb(1e-6, 2e-6, 0.1) # 0.25   (twice as far: a quarter)
charge.field(1e-9, 0.1)                       # 898.755   (N/C)
charge.potential(1e-9, 0.1)                   # 89.8755   (V)
charge.potential_energy(1e-6, -2e-6, 0.1)     # -0.179751   (J)
charge.force_on(2e-6, 500)                    # 0.001   (N)
charge.acceleration(constants.e, constants.m_e, 1000) # 1.75882e+14   (an electron in 1000 N/C, m/s²)

charge.flux(q_enclosed) charge.field_line(linear_density, r) charge.field_sheet(surface_density) charge.field_conductor(surface_density)

Gauss's law and symmetric charges. The electric flux out of any closed surface depends only on the charge inside it, flux(). For charges with a simple shape, Gauss's law gives the field directly: field_line() for a long line of charge, field_sheet() for a large flat sheet (the field does not depend on the distance), and field_conductor() just outside a conductor's surface, where the field is twice as strong for the same charge density.

Formulas and arguments
CallFormulaWhat each argument means
charge.flux(q_enclosed)
  • q_enclosed () — the total charge inside the closed surface, C
Returns the electric flux out of the surface, in N m²/C (Gauss's law)
charge.field_line(linear_density, r)
  • linear_density () — charge per metre of the line, C/m
  • r — distance from the line, m
Returns the electric field , in N/C
charge.field_sheet(surface_density)
  • surface_density () — charge per square metre of the sheet, C/m²
Returns the electric field near a large flat sheet, in N/C
charge.field_conductor(surface_density)
  • surface_density () — charge per square metre on the conductor's surface, C/m²
Returns the electric field just outside it, in N/C
charge.flux(1e-9)                             # 112.941   (N m²/C)
charge.field_line(1e-9, 0.1)                  # 179.751   (a line of 1 nC/m, 10 cm away, N/C)
charge.field_sheet(1e-6)                      # 56470.5   (a sheet of 1 µC/m², N/C)
charge.field_conductor(1e-6)                  # 112941   (twice the sheet)

charge.dipole_moment(q, separation) charge.dipole_torque(p, E, angle) charge.field_at(point, charges) charge.potential_at(point, charges) charge.plot_field(charges)

Dipoles and many charges. Two opposite charges close together form a dipole. dipole_moment() measures it and dipole_torque() is the turning effect of a field on it. For any collection of charges, give a list of (q, x, y) and field_at() adds up the field at a point as a pair (Ex, Ey), while potential_at() adds the potentials as ordinary numbers. plot_field() draws the field lines: they leave positive charges and end on negative ones.

Formulas and arguments
CallFormulaWhat each argument means
charge.dipole_moment(q, separation)
  • q — size of each charge, C
  • separation () — distance between the charges, m
Returns the dipole moment , in C m
charge.dipole_torque(p, E, angle)
  • p — dipole moment, C m
  • E — electric field, N/C
  • angle () — angle between the dipole and the field, degrees (default 90); use mode="rad" for radians
Returns the torque , in N m
charge.field_at(point, charges)
  • point () — where you want the field, m
  • charges () — a list of charges, each with its position
Returns the field as a pair , in N/C
charge.potential_at(point, charges)
  • point () — where you want the potential, m
  • charges () — a list of charges, each with its position
Returns the potential , in volts (potentials add as plain numbers)
charge.plot_field(charges, size, title, save)
  • charges () — list of charges (default: a dipole)
  • size — half-width of the picture, m
  • title — chart heading
  • save — file name to write the chart to
Returns draws the electric field lines
charge.dipole_moment(1e-9, 0.01)              # 1e-11   (C m)
charge.dipole_torque(1e-11, 1000)             # 1e-08   (N m)
charge.field_at((0, 0.1), [(1e-9, -0.1, 0), (-1e-9, 0.1, 0)]) # (635.516, 0)   (N/C, on the bisector)
charge.potential_at((0, 0.1), [(1e-9, -0.1, 0), (-1e-9, 0.1, 0)]) # 0   (0 on the bisector)
charge.plot_field([(1e-9, -0.1, 0), (-1e-9, 0.1, 0)])
charge.plot_field(): the field lines of a dipole run from the positive charge to the negative one
charge.plot_field(): the field lines of a dipole run from the positive charge to the negative one

capacitor.law(Q=, C=, V=) capacitor.parallel_plate(area, gap, k) capacitor.sphere(radius) capacitor.energy(C, V) capacitor.series(*caps) capacitor.parallel(*caps)

Capacitors. A capacitor stores charge on two plates, in proportion to the voltage across it: Q = C V, and law() solves that for whichever of the three you leave out. parallel_plate() gives the capacitance of two flat plates from their area and gap; k is the dielectric constant of the material between them, which multiplies it (water is about 80). sphere() is the capacitance of an isolated sphere. energy() and energy_density() give the energy stored. In parallel capacitors add, and in series their reciprocals add, the opposite of resistors. field_between() is the field between plates.

Formulas and arguments
CallFormulaWhat each argument means
capacitor.law(Q=, C=, V=)
  • Q — charge stored, C
  • C — capacitance, F
  • V — voltage across it, V (give two, it finds the third)
Returns the missing one of , ,
capacitor.parallel_plate(area, gap, k)
  • area () — area of each plate, m²
  • gap () — distance between the plates, m
  • k — dielectric constant of the material between them (default 1, vacuum)
Returns the capacitance , in F
capacitor.sphere(radius)
  • radius () — radius of the sphere, m
Returns the capacitance of the isolated sphere, in F
capacitor.energy(C, V)
  • C — capacitance, F
  • V — voltage across it, V
Returns the energy stored, in J
capacitor.series(*caps)
  • caps () — any number of capacitances, F
Returns the combined capacitance, in F
capacitor.parallel(*caps)
  • caps () — any number of capacitances, F
Returns the combined capacitance, in F
capacitor.field_between(V, gap)
  • V — voltage between the plates, V
  • gap () — distance between them, m
Returns the electric field between them, in V/m
capacitor.energy_density(E, k)
  • E — electric field, V/m
  • k — dielectric constant (default 1)
Returns the energy per cubic metre , in J/m³
capacitor.law(C=1e-6, V=12)                   # 1.2e-05   (charge, C)
capacitor.parallel_plate(0.01, 1e-3)          # 8.85419e-11   (two 10 cm × 10 cm plates 1 mm apart, F)
capacitor.parallel_plate(0.01, 1e-3, 80)      # 7.08335e-09   (with water between)
capacitor.sphere(6.371e6)                     # 0.000708869   (the Earth, F)
capacitor.energy(1e-6, 100)                   # 0.005   (J)
capacitor.series(1e-6, 2e-6, 3e-6)            # 5.45455e-07
capacitor.parallel(1e-6, 2e-6)                # 3e-06
capacitor.field_between(12, 0.001)            # 12000   (V/m)

circuit.ohm(V=, I=, R=) circuit.power_law(P=, V=, I=) circuit.resistance(rho, length, area) circuit.series(*R) circuit.parallel(*R) circuit.divider(...) circuit.terminal_voltage(...)

DC circuits. ohm() is Ohm's law, V = I R, and power_law() is P = V I; both solve for whichever variable you leave out. resistance() gives the resistance of a wire from its material, length and thickness. Resistors in series add (series()), and in parallel their reciprocals add (parallel()), so the combination is always smaller than the smallest. divider() and current_divider() split a voltage or a current between two resistors. power_resistor() and energy() give the heat produced and the energy used. terminal_voltage() is what a real battery delivers once its internal resistance is counted. drift_velocity() is the very slow average speed of the electrons in a wire.

Formulas and arguments
CallFormulaWhat each argument means
circuit.ohm(V=, I=, R=)
  • V — voltage, V
  • I — current, A
  • R — resistance, ohm (give two, it finds the third)
Returns the missing one of , ,
circuit.power_law(P=, V=, I=)
  • P — power, W
  • V — voltage, V
  • I — current, A (give two, it finds the third)
Returns the missing one of , ,
circuit.resistance(rho, length, area)
  • rho () — resistivity of the material, ohm m (copper is 1.68×10⁻⁸)
  • length () — length of the wire, m
  • area () — cross-section area of the wire, m²
Returns the resistance , in ohms
circuit.series(*resistors)
  • resistors () — any number of resistances, ohm
Returns the total resistance, in ohms
circuit.parallel(*resistors)
  • resistors () — any number of resistances, ohm
Returns the total resistance, in ohms
circuit.divider(V, R1, R2)
  • V — voltage across the pair, V
  • R1 () — first resistor, ohm
  • R2 () — second resistor, ohm
Returns the voltage across , in volts
circuit.current_divider(I, R1, R2)
  • I — total current, A
  • R1 () — first resistor, ohm
  • R2 () — second resistor in parallel with it, ohm
Returns the current through , in A
circuit.power_resistor(I, R)
  • I — current, A
  • R — resistance, ohm
Returns the heat power , in W
circuit.energy(P, t)
  • P — power, W
  • t — time, s
Returns the energy used, in J (divide by 3.6×10⁶ for kilowatt-hours)
circuit.terminal_voltage(emf, internal_resistance, R)
  • emf () — the battery's emf, V
  • internal_resistance () — resistance inside the battery, ohm
  • R — resistance of the load, ohm
Returns the voltage across the load, in volts
circuit.temperature_resistance(R0, alpha, dT)
  • R0 () — resistance at the starting temperature, ohm
  • alpha () — temperature coefficient, 1/K (copper: 0.0039)
  • dT () — temperature change, K
Returns the new resistance , in ohms
circuit.drift_velocity(I, n, area, q)
  • I — current, A
  • n — charge carriers per m³ (copper: 8.5×10²⁸)
  • area () — cross-section area of the wire, m²
  • q — charge of each carrier, C (default: an electron)
Returns the drift speed , in m/s
circuit.battery_life(capacity_Ah, I)
  • capacity_Ah — battery capacity, ampere-hours
  • I — current drawn, A
Returns the time the battery lasts, in hours
circuit.ohm(V=12, R=4)                        # 3   (A)
circuit.power_law(V=12, I=3)                  # 36   (W)
circuit.resistance(1.68e-8, 100, 1e-6)        # 1.68   (100 m of copper wire, 1 mm², ohm)
circuit.series(100, 200, 300)                 # 600
circuit.parallel(100, 200, 300)               # 54.5455
circuit.divider(12, 100, 200)                 # 8   (V across the 200 ohm resistor)
circuit.power_resistor(2, 5)                  # 20   (W)
circuit.energy(1000, 3600)                    # 3.6e+06   (J in an hour at 1 kW)
circuit.terminal_voltage(12, 0.5, 10)         # 11.4286   (V)
circuit.temperature_resistance(100, 0.0039, 50) # 119.5
circuit.drift_velocity(10, 8.5e28, 1e-6)      # 0.000734295   (m/s: less than a millimetre a second)
circuit.battery_life(2.5, 0.5)                # 5   (hours)

circuit.time_constant(R, C) circuit.rc_charge(V, R, C, t) circuit.rc_discharge(V0, R, C, t) circuit.rl_current(V, R, L, t) circuit.plot_rc()

Charging and discharging. A capacitor charged through a resistor does not fill instantly: its voltage follows V (1 - e^(-t/RC)), reaching 63% of the supply after one time constant R C. rc_charge() and rc_discharge() give the voltage at any time, and time_constant() the time constant. An inductor switched on through a resistor behaves in the same way, with L/R in place of R C: rl_current() gives the current. plot_rc() draws a capacitor charging and discharging.

Formulas and arguments
CallFormulaWhat each argument means
circuit.time_constant(R, C)
  • R — resistance, ohm
  • C — capacitance, F
Returns the time constant , in s
circuit.rc_charge(V, R, C, t)
  • V — supply voltage, V
  • R — resistance, ohm
  • C — capacitance, F
  • t — time since switching on, s
Returns the voltage across the capacitor, in volts
circuit.rc_discharge(V0, R, C, t)
  • V0 () — starting voltage across the capacitor, V
  • R — resistance, ohm
  • C — capacitance, F
  • t — time since discharging began, s
Returns the voltage across the capacitor, in volts
circuit.rl_current(V, R, L, t)
  • V — supply voltage, V
  • R — resistance, ohm
  • L — inductance, H
  • t — time since switching on, s
Returns the current in the circuit, in A
circuit.plot_rc(V, R, C, title, save)
  • V — supply voltage, V
  • R — resistance, ohm
  • C — capacitance, F
  • title — chart heading
  • save — file name to write the chart to
Returns draws the charging and discharging curves with the time constant marked
circuit.time_constant(1000, 1e-3)             # 1   (s)
circuit.rc_charge(10, 1000, 1e-3, 1)          # 6.32121   (V after one time constant)
circuit.rc_discharge(10, 1000, 1e-3, 1)       # 3.67879   (V)
circuit.rc_charge(10, 1000, 1e-3, 5)          # 9.93262   (nearly full after five)
circuit.rl_current(12, 6, 0.5, 0.1)           # 1.39761   (A)
circuit.plot_rc(10, 1000, 1e-3)
circuit.plot_rc(): after one time constant a charging capacitor is at 63% and a discharging one at 37%
circuit.plot_rc(): after one time constant a charging capacitor is at 63% and a discharging one at 37%

magnet.force_on_charge(q, v, B, angle) magnet.cyclotron_radius(m, v, q, B) magnet.cyclotron_frequency(q, B, m) magnet.helix_pitch(...) magnet.velocity_selector(E, B) magnet.force_on_wire(B, I, L, angle) magnet.plot_helix()

Magnetic forces. A magnetic field pushes a moving charge sideways, with a force q v B sin θ at right angles to the motion, so it never speeds the charge up, it only turns it. A charge moving across the field goes round in a circle of radius cyclotron_radius(), and the number of turns a second, cyclotron_frequency(), does not depend on the speed (that is how a cyclotron works). A charge moving at an angle to the field spirals along it, helix_pitch() being the distance gone in one turn. velocity_selector() is the speed that passes straight through crossed electric and magnetic fields. force_on_wire() is the force on a current-carrying wire. plot_helix() draws the spiral.

Formulas and arguments
CallFormulaWhat each argument means
magnet.force_on_charge(q, v, B, angle)
  • q — the charge, C (a proton has 1.6×10⁻¹⁹)
  • v — speed of the charge, m/s
  • B — magnetic field, T
  • angle () — angle between the velocity and the field, degrees (default 90); use mode="rad" for radians
Returns the force , in N (always at right angles to the motion)
magnet.cyclotron_radius(m, v, q, B)
  • m — mass of the particle, kg
  • v — its speed, m/s
  • q — its charge, C
  • B — magnetic field, T
Returns the radius of its circle, in m
magnet.cyclotron_frequency(q, B, m)
  • q — charge, C
  • B — magnetic field, T
  • m — mass, kg
Returns the number of turns per second , in Hz (it does not depend on the speed)
magnet.helix_pitch(m, v_parallel, q, B)
  • m — mass, kg
  • v_parallel () — speed along the field, m/s
  • q — charge, C
  • B — magnetic field, T
Returns the distance moved along the field in one turn, in m
magnet.velocity_selector(E, B)
  • E — electric field, V/m
  • B — magnetic field at right angles to it, T
Returns the speed that passes straight through, in m/s
magnet.force_on_wire(B, I, length, angle)
  • B — magnetic field, T
  • I — current, A
  • length () — length of wire in the field, m
  • angle () — angle between the wire and the field, degrees (default 90); use mode="rad" for radians
Returns the force on the wire, in N
magnet.plot_helix(m, v, q, B, angle, turns, title, save)
  • m — mass, kg (default: an electron)
  • v — speed, m/s
  • q — charge, C
  • B — magnetic field, T
  • angle () — launch angle to the field, degrees
  • turns — number of turns to draw
  • title — chart heading
  • save — file name to write the chart to
Returns draws the spiral path of the particle
magnet.force_on_charge(constants.e, 1e6, 0.5) # 8.01088e-14   (N)
magnet.cyclotron_radius(constants.m_e, 1e7, constants.e, 0.01) # 0.00568563   (m)
magnet.cyclotron_radius(constants.m_p, 1e7, constants.e, 0.01) # 10.4397   (a proton: 1836 times bigger, m)
magnet.cyclotron_frequency(constants.e, 0.01, constants.m_e) # 2.79925e+08   (Hz)
magnet.velocity_selector(1000, 0.01)          # 100000   (m/s)
magnet.force_on_wire(0.5, 10, 0.2)            # 1   (N)
magnet.plot_helix(constants.m_e, 1e7, constants.e, 0.01, angle=60)
magnet.plot_helix(): a charge launched at an angle to the field spirals along it
magnet.plot_helix(): a charge launched at an angle to the field spirals along it

magnet.field_wire(I, r) magnet.field_loop(I, radius, turns) magnet.field_solenoid(n, I) magnet.field_toroid(N, I, r) magnet.force_between_wires(...) magnet.torque_on_loop(...) magnet.hall_voltage(...)

Magnetic fields from currents. A current makes a magnetic field. field_wire() is the field round a long straight wire, field_loop() the field at the centre of a coil, field_solenoid() the nearly uniform field inside a long coil (turns per metre), and field_toroid() the field inside a ring-shaped coil. Two parallel wires carrying currents in the same direction attract: force_between_wires() gives the force, and is the basis of the definition of the ampere. magnetic_moment() and torque_on_loop() describe a coil in a field, which is how an electric motor turns. hall_voltage() is the sideways voltage across a current-carrying strip in a field, used to measure fields and to find the density of charge carriers.

Formulas and arguments
CallFormulaWhat each argument means
magnet.field_wire(I, r)
  • I — current in the wire, A
  • r — distance from the wire, m
Returns the magnetic field , in T ( T m/A)
magnet.field_loop(I, radius, turns)
  • I — current, A
  • radius () — radius of the coil, m
  • turns () — number of turns (default 1)
Returns the field at the centre of the coil, in T
magnet.field_solenoid(n, I)
  • n — turns per metre of the solenoid
  • I — current, A
Returns the field inside a long solenoid, in T
magnet.field_toroid(N, I, r)
  • N — total number of turns
  • I — current, A
  • r — distance from the centre of the ring, m
Returns the field inside the toroid, in T
magnet.force_between_wires(I1, I2, length, d)
  • I1 () — current in the first wire, A
  • I2 () — current in the second wire, A (negative if opposite)
  • length () — length of wire, m
  • d — distance between the wires, m
Returns the force between them, in N: positive attracts, negative repels
magnet.magnetic_moment(N, I, area)
  • N — number of turns
  • I — current, A
  • area () — area of the coil, m²
Returns the magnetic moment , in A m²
magnet.torque_on_loop(N, I, area, B, angle)
  • N — number of turns
  • I — current, A
  • area () — area of the coil, m²
  • B — magnetic field, T
  • angle () — angle between the field and the coil's axis, degrees (default 90)
Returns the torque , in N m
magnet.hall_voltage(I, B, n, thickness, q)
  • I — current along the strip, A
  • B — magnetic field through it, T
  • n — charge carriers per m³
  • thickness () — thickness of the strip in the direction of the field, m
  • q — charge of each carrier, C (default: an electron)
Returns the Hall voltage across the strip, in volts
magnet.field_wire(10, 0.05)                   # 4e-05   (T, 10 A at 5 cm)
magnet.field_loop(2, 0.1, 50)                 # 0.000628319   (T, 50 turns of radius 10 cm)
magnet.field_solenoid(1000, 2)                # 0.00251327   (T, 1000 turns per metre)
magnet.field_toroid(1000, 2, 0.1)             # 0.004   (T)
magnet.force_between_wires(10, 10, 1, 0.01)   # 0.002   (N, per metre of wire)
magnet.magnetic_moment(100, 2, 0.01)          # 2   (A m²)
magnet.torque_on_loop(100, 2, 0.01, 0.5)      # 1   (N m)
magnet.hall_voltage(1, 1, 8.5e28, 1e-4)       # 7.34295e-07   (V)

induction.flux(B, area, angle) induction.emf(N, dPhi, dt) induction.motional_emf(B, L, v) induction.generator_peak(N, B, area, omega) induction.generator(..., t)

Induction. A changing magnetic flux induces a voltage (Faraday's law), and the minus sign in emf() is Lenz's law: the voltage drives a current that opposes the change. flux() is the field through a loop (the angle is from the loop's normal). motional_emf() is the voltage across a rod moved through a field. A coil spun in a field makes an alternating voltage, the principle of every generator: generator_peak() gives its peak, generator() its value at any moment.

Formulas and arguments
CallFormulaWhat each argument means
induction.flux(B, area, angle)
  • B — magnetic field, T
  • area () — area of the loop, m²
  • angle () — angle between the field and the loop's normal, degrees (default 0); use mode="rad" for radians
Returns the magnetic flux , in Wb
induction.emf(N, dPhi, dt)
  • N — number of turns in the coil
  • dPhi () — change in magnetic flux, Wb
  • dt () — time over which it changes, s
Returns the induced emf , in volts (the minus sign is Lenz's law)
induction.motional_emf(B, length, v)
  • B — magnetic field, T
  • length () — length of the moving rod, m
  • v — its speed across the field, m/s
Returns the emf across the rod, in volts
induction.generator_peak(N, B, area, omega)
  • N — number of turns
  • B — magnetic field, T
  • area () — area of the coil, m²
  • omega () — angular speed, rad/s
Returns the greatest emf , in volts
induction.generator(N, B, area, omega, t)
  • N — number of turns
  • B — magnetic field, T
  • area () — area of the coil, m²
  • omega () — angular speed, rad/s
  • t — time, s
Returns the emf at time , in volts
induction.flux(0.5, 0.02)                     # 0.01   (Wb)
induction.emf(100, 0.002, 0.1)                # -2   (V: the flux rises by 2 mWb in 0.1 s)
induction.motional_emf(0.5, 0.2, 3)           # 0.3   (V)
induction.generator_peak(200, 0.5, 0.02, 314.159) # 628.318   (V: 200 turns at 50 Hz)
induction.generator(200, 0.5, 0.02, 314.159, 0.005) # 628.318   (V a quarter-turn in (peak))

induction.inductance_solenoid(N, area, length, mu_r) induction.self_emf(L, dI, dt) induction.inductor_energy(L, I) induction.field_energy_density(B) induction.transformer(V, Np, Ns) induction.transformer_current(...)

Inductors and transformers. A coil resists changes in its own current, by an amount called its inductance. inductance_solenoid() gives it from the turns, area, length and core material (mu_r, about 1000 for iron). self_emf() is the voltage across an inductor when its current changes, and inductor_energy() and field_energy_density() give the energy stored in the magnetic field. A transformer changes an alternating voltage in the ratio of the turns on its two coils: transformer() gives the voltage out, and transformer_current() the current, which changes the opposite way so that power is conserved.

Formulas and arguments
CallFormulaWhat each argument means
induction.inductance_solenoid(N, area, length, mu_r)
  • N — number of turns
  • area () — cross-section area of the coil, m²
  • length () — length of the coil, m
  • mu_r () — relative permeability of the core (default 1, air)
Returns the inductance , in H
induction.self_emf(L, dI, dt)
  • L — inductance, H
  • dI () — change in current, A
  • dt () — time over which it changes, s
Returns the emf across the inductor, in volts
induction.inductor_energy(L, I)
  • L — inductance, H
  • I — current, A
Returns the energy stored, in J
induction.field_energy_density(B, mu_r)
  • B — magnetic field, T
  • mu_r () — relative permeability (default 1)
Returns the energy per cubic metre , in J/m³
induction.transformer(V_primary, N_primary, N_secondary)
  • V_primary () — voltage into the primary coil, V
  • N_primary () — turns on the primary coil
  • N_secondary () — turns on the secondary coil
Returns the voltage out of the secondary, in volts
induction.transformer_current(I_primary, N_primary, N_secondary)
  • I_primary () — current in the primary coil, A
  • N_primary () — turns on the primary
  • N_secondary () — turns on the secondary
Returns the current in the secondary, in A (power in equals power out)
induction.inductance_solenoid(500, 1e-4, 0.2) # 0.00015708   (H)
induction.inductance_solenoid(500, 1e-4, 0.2, 1000) # 0.15708   (with an iron core)
induction.self_emf(0.5, 2, 0.1)               # -10   (V)
induction.inductor_energy(0.5, 2)             # 1   (J)
induction.field_energy_density(1)             # 397887   (J/m³ in a 1 T field)
induction.transformer(230, 1000, 50)          # 11.5   (V out)
induction.transformer_current(1, 1000, 50)    # 20   (A out)

ac.rms(peak) ac.peak(rms) ac.reactance_c(C, f) ac.reactance_l(L, f) ac.impedance(R, L, C, f) ac.phase(R, L, C, f) ac.current_rms(V, R, L, C, f)

Alternating current. Mains electricity swings back and forth, and its "230 V" is the rms value, the steady voltage that would heat the same: rms() divides the peak by √2, and peak() goes the other way. A capacitor and an inductor oppose alternating current by an amount called their reactance: a capacitor's, reactance_c(), falls as the frequency rises, and an inductor's, reactance_l(), rises. A series R, L, C circuit has an impedance from all three, impedance(), and the current lags the voltage by phase() (negative means it leads, as for a capacitor). Use L=0 or C=None to leave a part out.

Formulas and arguments
CallFormulaWhat each argument means
ac.rms(peak)
  • peak () — peak value of the voltage or current
Returns the rms value, the steady value with the same heating effect
ac.peak(rms)
  • rms () — rms value of the voltage or current
Returns the peak value
ac.reactance_c(C, f)
  • C — capacitance, F
  • f — frequency, Hz
Returns the capacitive reactance , in ohms (falls as the frequency rises)
ac.reactance_l(L, f)
  • L — inductance, H
  • f — frequency, Hz
Returns the inductive reactance , in ohms (rises with the frequency)
ac.impedance(R, L, C, f)
  • R — resistance, ohm
  • L — inductance, H (use 0 to leave it out)
  • C — capacitance, F (use None to leave it out)
  • f — frequency, Hz
Returns the impedance , in ohms
ac.phase(R, L, C, f)
  • R — resistance, ohm
  • L — inductance, H (use 0 to leave it out)
  • C — capacitance, F (use None to leave it out)
  • f — frequency, Hz
  • mode — "deg" (default) or "rad"
Returns the angle by which the current lags the voltage (negative: it leads)
ac.current_rms(V_rms, R, L, C, f)
  • V_rms () — rms voltage of the supply, V
  • R — resistance, ohm
  • L — inductance, H (use 0 to leave it out)
  • C — capacitance, F (use None to leave it out)
  • f — frequency, Hz
Returns the rms current , in A
ac.rms(325)                                   # 229.81   (V: the peak of a 230 V mains supply)
ac.peak(230)                                  # 325.269
ac.reactance_c(1e-6, 50)                      # 3183.1   (ohm)
ac.reactance_l(0.1, 50)                       # 31.4159   (ohm)
ac.impedance(100, 0.1, 1e-6, 50)              # 3153.27   (ohm)
ac.impedance(100, 0, None, 50)                # 100   (just the resistor)
ac.phase(100, 0.1, 1e-6, 50)                  # -88.1827   (degrees: mostly capacitor, so the current leads)
ac.current_rms(230, 100, 0.1, 1e-6, 50)       # 0.0729402   (A)

ac.resonance(L, C) ac.q_factor(R, L, C) ac.power(V_rms, I_rms, phase) ac.plot_resonance(L, C, resistances)

Resonance and power. At one frequency the reactances of the inductor and capacitor cancel exactly, the impedance falls to just R, and the current is at its greatest. This resonant frequency is resonance(), and it is how a radio is tuned to a station. The Q factor, q_factor(), says how sharp the peak is: low resistance gives a tall, narrow peak. power() gives the average power, which includes a cos φ "power factor": a circuit whose current is out of step with its voltage does less useful work for the same current. plot_resonance() draws the current against frequency for several resistances.

Formulas and arguments
CallFormulaWhat each argument means
ac.resonance(L, C)
  • L — inductance, H
  • C — capacitance, F
Returns the resonant frequency , in Hz
ac.q_factor(R, L, C)
  • R — resistance, ohm
  • L — inductance, H
  • C — capacitance, F
Returns the quality factor (no unit): a larger is a sharper resonance
ac.power(V_rms, I_rms, phase)
  • V_rms () — rms voltage, V
  • I_rms () — rms current, A
  • phase () — phase angle between them, degrees (default 0)
Returns the average power , in W
ac.plot_resonance(L, C, resistances, V_rms, title, save)
  • L — inductance, H
  • C — capacitance, F
  • resistances () — list of resistances to draw, ohm
  • V_rms () — rms supply voltage, V
  • title — chart heading
  • save — file name to write the chart to
Returns draws current against frequency for each resistance
ac.resonance(0.1, 1e-6)                       # 503.292   (Hz)
ac.q_factor(10, 0.1, 1e-6)                    # 31.6228
ac.impedance(100, 0.1, 1e-6, 503.292)         # 100   (at resonance, just R)
ac.power(230, 5, 30)                          # 995.929   (W)
ac.plot_resonance(0.1, 1e-6, resistances=(10, 30, 100))
ac.plot_resonance(): less resistance gives a taller, sharper peak at the resonant frequency
ac.plot_resonance(): less resistance gives a taller, sharper peak at the resonant frequency

emwave.wavelength(f) emwave.frequency(wavelength) emwave.band(...) emwave.intensity(E0) emwave.peak_field(I) emwave.radiation_pressure(I)

Electromagnetic waves. Light, radio, microwaves, X-rays and the rest are all the same thing: electric and magnetic fields travelling together at the speed of light. wavelength() and frequency() convert between the two (their product is always c), and band() names the kind of wave from either one. In a wave the electric and magnetic fields are in a fixed ratio, E = c B (b_from_e(), e_from_b()). intensity() gives the power carried per square metre from the peak electric field, and peak_field() goes back. intensity_spread() shows how it thins out with distance from a source. Light pushes on what it hits, radiation_pressure(), twice as hard if it is reflected. speed_in_medium() is the slower speed in glass or water.

Formulas and arguments
CallFormulaWhat each argument means
emwave.wavelength(frequency)
  • frequency () — frequency of the wave, Hz
Returns the wavelength , in m ( m/s)
emwave.frequency(wavelength)
  • wavelength () — wavelength, m
Returns the frequency , in Hz
emwave.band(wavelength=, frequency=)
  • wavelength () — wavelength, m (give exactly one of the two)
  • frequency () — frequency, Hz
Returns the name of the kind of wave: "radio wave", "microwave", "infrared", "visible light", "ultraviolet", "X-ray" or "gamma ray"
emwave.b_from_e(E)
  • E — electric field of the wave, V/m
Returns the magnetic field of the wave, in T
emwave.e_from_b(B)
  • B — magnetic field of the wave, T
Returns the electric field of the wave, in V/m
emwave.intensity(E0)
  • E0 () — peak electric field of the wave, V/m
Returns the intensity (power per square metre), in W/m²
emwave.peak_field(intensity)
  • intensity () — power per square metre, W/m²
Returns the peak electric field , in V/m
emwave.intensity_spread(power, r)
  • power () — power of the source, W
  • r — distance from it, m
Returns the intensity at that distance, in W/m²
emwave.radiation_pressure(intensity, reflect)
  • intensity () — power per square metre, W/m²
  • reflect — False (default) if the surface absorbs the light, True if it reflects it
Returns the pressure of the light, in Pa
emwave.speed_in_medium(n)
  • n — refractive index of the material (water 1.33, glass about 1.5)
Returns the speed of light in it, in m/s
emwave.wavelength(100e6)                      # 2.99792   (an FM radio station at 100 MHz, m)
emwave.frequency(500e-9)                      # 5.99585e+14   (green light, Hz)
emwave.band(wavelength=500e-9)                # 'visible light'
emwave.band(frequency=2.4e9)                  # 'microwave'   (WiFi)
emwave.band(wavelength=1e-10)                 # 'X-ray'
emwave.b_from_e(300)                          # 1.00069e-06   (T)
emwave.intensity(100)                         # 13.2721   (W/m²)
emwave.peak_field(1000)                       # 868.021   (V/m in sunlight-strength light)
emwave.intensity_spread(100, 2)               # 1.98944   (W/m²)
emwave.radiation_pressure(1000)               # 3.33564e-06   (Pa)
emwave.speed_in_medium(1.5)                   # 1.99862e+08   (m/s in glass)

Relativity

Time dilation and length contraction, spacetime, relativistic momentum and energy, the Doppler effect for light, black holes, and gravity's effect on clocks and light. Every function is shown with its formula, and each table says which argument is which letter of it.

lorentz.gamma(v) lorentz.time_dilation(proper_time, v) relmech.energy_total(m, v) doppler.longitudinal(f, v) gr.schwarzschild_radius(M)

How the relativity functions work. Everything is in SI units: metres, seconds, kilograms, joules, and speeds in metres per second. There are four namespaces: lorentz (time, length, velocities and spacetime in special relativity), relmech (relativistic momentum, energy, decays and collisions), doppler (the Doppler effect and redshift of light) and gr (black holes, gravity's effect on clocks and light).

Speeds as a fraction of light. Relativity is about speeds close to the speed of light, constants.c = 299,792,458 m/s, so lorentz.speed(0.8) turns the fraction 0.8 into metres per second, and lorentz.beta(v) turns it back. A speed at or above c is refused with a plain message. In gr, wherever a mass is needed you can pass a body name such as 'sun' or 'earth'. Every function below is shown with its formula, and each table says which argument is which letter of it.

lorentz.speed(0.8)                            # 2.39834e+08   (m/s)
lorentz.gamma(lorentz.speed(0.8))             # 1.66667
constants.c                                   # 2.99792e+08
gr.schwarzschild_radius('sun')                # 2954.13   (m)

lorentz.gamma(v) lorentz.time_dilation(proper_time, v) lorentz.length_contraction(proper_length, v) lorentz.velocity_addition(u, v) lorentz.rapidity(v)

Time, length and speed at high speed. Everything starts with the Lorentz factor gamma(v), which is 1 at rest and grows without limit as the speed approaches that of light. A moving clock ticks slowly: time_dilation() multiplies the time on the clock itself (its proper time) by gamma. A moving object is shorter along its direction of motion: length_contraction() divides its length at rest by gamma. Speeds do not simply add: velocity_addition() combines two speeds and never gives more than the speed of light. rapidity() is the quantity that does add. speed_from_gamma() goes the other way.

Formulas and arguments
CallFormulaWhat each argument means
lorentz.speed(beta)
  • beta () — the speed as a fraction of the speed of light (between −1 and 1)
Returns the speed , in m/s ( m/s)
lorentz.beta(v)
  • v — speed, m/s (use lorentz.speed(0.8) for 0.8 c)
Returns the speed as a fraction of the speed of light
lorentz.gamma(v)
  • v — speed, m/s (use lorentz.speed(0.8) for 0.8 c)
Returns the Lorentz factor (no unit): 1 at rest, and growing without limit as approaches
lorentz.speed_from_gamma(gamma)
  • gamma () — the Lorentz factor (at least 1)
Returns the speed that has that Lorentz factor, in m/s
lorentz.time_dilation(proper_time, v)
  • proper_time () — time between two ticks on the moving clock itself, s
  • v — speed, m/s (use lorentz.speed(0.8) for 0.8 c)
Returns the time between those ticks as measured by you, in s (longer than )
lorentz.length_contraction(proper_length, v)
  • proper_length () — length of the object at rest, m
  • v — speed, m/s (use lorentz.speed(0.8) for 0.8 c) of the object along its length
Returns the length you measure, in m (shorter than )
lorentz.velocity_addition(u, v)
  • u — speed of an object in a frame, m/s
  • v — speed of that frame relative to you, m/s
Returns the speed of the object relative to you, in m/s (never more than )
lorentz.rapidity(v)
  • v — speed, m/s (use lorentz.speed(0.8) for 0.8 c)
Returns the rapidity (no unit): rapidities of successive boosts simply add
lorentz.gamma(lorentz.speed(0.6))             # 1.25
lorentz.gamma(lorentz.speed(0.99))            # 7.08881
lorentz.time_dilation(1, lorentz.speed(0.8))  # 1.66667   (s on your clock for each 1 s on theirs)
lorentz.length_contraction(100, lorentz.speed(0.8)) # 60   (a 100 m spaceship, m)
lorentz.velocity_addition(lorentz.speed(0.8), lorentz.speed(0.8)) / constants.c # 0.97561   (not 1.6)
lorentz.rapidity(lorentz.speed(0.8))          # 1.09861
lorentz.speed_from_gamma(2) / constants.c     # 0.866025

lorentz.transform(x, t, v) lorentz.interval(dt, dx) lorentz.proper_time(dt, dx) lorentz.simultaneity(dx, v)

Events in spacetime. An event is a place and a time. transform() is the Lorentz transformation: it gives where and when the same event happens for an observer moving along x at speed v. Observers disagree about distances and times, but all agree on the spacetime interval, interval(): when it is positive one event can cause the other, and proper_time() gives the time on a clock that travels between them. simultaneity() shows that events that happen at the same time in one frame do not in another.

Formulas and arguments
CallFormulaWhat each argument means
lorentz.transform(x, t, v)
  • x — position of the event, m
  • t — time of the event, s
  • v — speed of the new frame along x, m/s
Returns the pair : where and when the same event happens in the moving frame
lorentz.interval(dt, dx)
  • dt () — time between two events, s
  • dx () — distance between them, m
Returns the spacetime interval , in m²: the same for every observer (positive: one event can cause the other)
lorentz.proper_time(dt, dx)
  • dt () — time between two events, s
  • dx () — distance between them, m
Returns the time on a clock that travels from one event to the other, in s
lorentz.simultaneity(dx, v)
  • dx () — distance between two events that are simultaneous in the first frame, m
  • v — speed, m/s (use lorentz.speed(0.8) for 0.8 c) of the new frame
Returns how far apart in time the two events seem in the moving frame, in s
lorentz.transform(1000, 0, lorentz.speed(0.6)) # (1250, -2.50173e-06)   ((x', t'))
lorentz.interval(1, 1e8)                      # 7.98755e+16   (m²)
lorentz.proper_time(10, 1e9)                  # 9.42727   (s)
lorentz.simultaneity(1000, lorentz.speed(0.6)) # -2.50173e-06   (s)

lorentz.twin_trip(distance, v) lorentz.muon_fraction(distance, v) lorentz.plot_gamma() lorentz.plot_minkowski(v)

The twins and the muons. twin_trip() is the twin paradox: a traveller goes out to a star and back at speed v and returns having aged less than the twin who stayed. It returns both times. muon_fraction() is the classic test of time dilation: muons made high in the atmosphere live only 2.2 µs, which would not get them to the ground even at nearly the speed of light, yet many arrive, because their clocks run slowly. plot_gamma() draws the Lorentz factor, and plot_minkowski() a spacetime diagram, in which a moving frame's axes tilt towards the line of light.

Formulas and arguments
CallFormulaWhat each argument means
lorentz.twin_trip(distance, v)
  • distance () — distance to the turning point, m
  • v — speed, m/s (use lorentz.speed(0.8) for 0.8 c)
Returns the pair for the round trip, in s: the traveller ages less
lorentz.muon_fraction(distance, v, lifetime)
  • distance () — distance travelled, m
  • v — speed of the muons, m/s
  • lifetime () — lifetime of a muon at rest, s (default 2.197×10⁻⁶)
Returns the fraction still alive when they arrive
lorentz.plot_gamma(title, save)
  • title — chart heading
  • save — file name to write the chart to
Returns draws the Lorentz factor against speed
lorentz.plot_minkowski(v, title, save)
  • v — speed of the moving frame, m/s (default 0.5 c)
  • title — chart heading
  • save — file name to write the chart to
Returns draws a spacetime diagram with the tilted axes of the moving frame and the light lines
lorentz.twin_trip(4.24 * 9.4607e15, lorentz.speed(0.9)) # (2.97342e+08, 1.29608e+08)   (a round trip to Proxima Centauri: (Earth seconds, traveller seconds))
lorentz.muon_fraction(10000, lorentz.speed(0.98)) # 0.0458213   (10 km down)
lorentz.muon_fraction(10000, 1000)            # 0   (at ordinary speed: none survive)
lorentz.plot_gamma()
lorentz.plot_minkowski()
lorentz.plot_gamma(): the Lorentz factor is close to 1 until the speed is a large fraction of the speed of light, then climbs steeply
lorentz.plot_gamma(): the Lorentz factor is close to 1 until the speed is a large fraction of the speed of light, then climbs steeply
lorentz.plot_minkowski(): a frame moving at half the speed of light has axes tilted towards the light line
lorentz.plot_minkowski(): a frame moving at half the speed of light has axes tilted towards the light line

relmech.momentum(m, v) relmech.energy_total(m, v) relmech.rest_energy(m) relmech.kinetic(m, v) relmech.energy_momentum(p, m) relmech.accelerate(m, q, V) relmech.plot_kinetic()

Momentum and energy near the speed of light. At high speed momentum is γ m v and the total energy is γ m c². At rest this leaves the rest energy m c², and what is left over is the kinetic energy, kinetic(), which for slow speeds is the familiar ½ m v². energy_momentum() links energy, momentum and mass without needing the speed, and speed_from_momentum() and speed_from_kinetic() find the speed back. accelerate() gives the speed a charged particle reaches after falling through a voltage. plot_kinetic() shows why the speed of light cannot be reached: the energy needed grows without limit.

Formulas and arguments
CallFormulaWhat each argument means
relmech.momentum(m, v)
  • m — mass, kg
  • v — speed, m/s (use lorentz.speed(0.8) for 0.8 c)
Returns the relativistic momentum , in kg m/s
relmech.rest_energy(m)
  • m — mass, kg
Returns the energy of the mass at rest, in J
relmech.energy_total(m, v)
  • m — mass, kg
  • v — speed, m/s (use lorentz.speed(0.8) for 0.8 c)
Returns the total energy , in J
relmech.kinetic(m, v)
  • m — mass, kg
  • v — speed, m/s (use lorentz.speed(0.8) for 0.8 c)
Returns the kinetic energy , in J (close to when )
relmech.energy_momentum(p, m)
  • p — momentum, kg m/s
  • m — mass, kg
Returns the total energy , in J
relmech.speed_from_momentum(p, m)
  • p — momentum, kg m/s
  • m — mass, kg (used to find the total energy E)
Returns the speed , in m/s
relmech.speed_from_kinetic(m, K)
  • m — mass, kg
  • K — kinetic energy, J
Returns the speed , in m/s
relmech.accelerate(m, q, V)
  • m — mass of the particle, kg
  • q — its charge, C
  • V — voltage it is accelerated through, V
Returns the speed it reaches from rest, in m/s
relmech.plot_kinetic(m, title, save)
  • m — mass, kg (default: an electron)
  • title — chart heading
  • save — file name to write the chart to
Returns draws relativistic and classical kinetic energy against speed
relmech.rest_energy(constants.m_e) / constants.eV # 510999   (eV: an electron)
relmech.momentum(constants.m_e, lorentz.speed(0.9)) # 5.63865e-22   (kg m/s)
relmech.energy_total(constants.m_e, lorentz.speed(0.9)) / constants.eV # 1.17231e+06   (eV)
relmech.kinetic(constants.m_e, lorentz.speed(0.9)) / constants.eV # 661313   (eV)
relmech.kinetic(constants.m_e, 100)           # 4.56294e-27   (slow: the same as ½ m v²)
relmech.speed_from_kinetic(constants.m_e, 1e6 * constants.eV) / constants.c # 0.941079   (a 1 MeV electron)
relmech.accelerate(constants.m_e, constants.e, 1e6) / constants.c # 0.941079   (through a million volts)
relmech.accelerate(constants.m_e, constants.e, 100) / constants.c # 0.0197807   (through 100 V)
relmech.energy_momentum(1e-22, constants.m_e) # 8.71873e-14   (J)
relmech.speed_from_momentum(1e-22, constants.m_e) / constants.c # 0.343849
relmech.plot_kinetic()
relmech.plot_kinetic(): the classical ½ m v² gives up at high speed; the relativistic kinetic energy grows without limit as v approaches c
relmech.plot_kinetic(): the classical ½ m v² gives up at high speed; the relativistic kinetic energy grows without limit as v approaches c

relmech.invariant_mass(E, p) relmech.mass_from_energy(E) relmech.decay_energy(M, m1, m2) relmech.threshold_energy(m_beam, m_target, m_products) relmech.cherenkov_angle(n, v)

Mass, decays and collisions. Mass and energy are the same thing counted in different units: mass_from_energy() converts one to the other. invariant_mass() is the mass of a system worked out from its total energy and momentum, which every observer agrees on. decay_energy() is the energy released when a particle breaks into two lighter ones, and threshold_energy() is the least kinetic energy a beam particle needs to make new ones when it hits a target at rest. A charged particle that moves faster than light does in a medium gives off a cone of blue light, the Cherenkov effect: cherenkov_angle() is the angle of the cone and cherenkov_threshold() the slowest speed that makes it.

Formulas and arguments
CallFormulaWhat each argument means
relmech.invariant_mass(E, p)
  • E — total energy of the system, J
  • p — total momentum, kg m/s
Returns the mass of the system, in kg: the same for every observer
relmech.mass_from_energy(E)
  • E — energy, J
Returns the mass that corresponds to it, in kg
relmech.decay_energy(M, m1, m2)
  • M — mass of the particle that decays, kg
  • m1 () — mass of the first product, kg
  • m2 () — mass of the second product, kg
Returns the energy released, in J (shared as motion of the products)
relmech.threshold_energy(m_beam, m_target, m_products)
  • m_beam () — mass of the particle fired, kg
  • m_target () — mass of the target at rest, kg
  • m_products () — total mass of everything made, kg
Returns the least kinetic energy the beam needs for the reaction to happen, in J
relmech.cherenkov_angle(n, v)
  • n — refractive index of the medium
  • v — speed, m/s (use lorentz.speed(0.8) for 0.8 c) of the charged particle
  • mode — "deg" (default) or "rad"
Returns the half-angle of the cone of light, in degrees
relmech.cherenkov_threshold(n)
  • n — refractive index of the medium
Returns the slowest speed that makes Cherenkov light, in m/s
relmech.mass_from_energy(1)                   # 1.11265e-17   (kg for 1 J)
relmech.invariant_mass(relmech.energy_total(constants.m_e, lorentz.speed(0.9)), relmech.momentum(constants.m_e, lorentz.speed(0.9))) # 9.10938e-31   (the electron mass again)
relmech.decay_energy(constants.m_n, constants.m_p, constants.m_e) / constants.eV # 782333   (a neutron decay, eV (the neutrino is ignored))
relmech.threshold_energy(constants.m_p, constants.m_p, 4 * constants.m_p) / (1e9 * constants.eV) # 5.62963   (making a proton–antiproton pair, GeV)
relmech.cherenkov_angle(1.33, lorentz.speed(0.99)) # 40.5821   (degrees, in water)
relmech.cherenkov_threshold(1.33) / constants.c # 0.75188

doppler.longitudinal(f, v, approaching) doppler.wavelength(wavelength, v) doppler.transverse(f, v) doppler.redshift(v) doppler.speed_from_redshift(z) doppler.aberration(angle, v)

The Doppler effect for light. A source moving towards you has its light shifted to a higher frequency, and one moving away to a lower one; for light this uses the relativistic formulas in longitudinal() and wavelength(), so it works at any speed and does not need a medium. transverse() is the effect with no classical counterpart: a source moving across your line of sight still has its frequency lowered, by time dilation. redshift() is the number astronomers quote, z, and speed_from_redshift() reverses it. aberration() shows how your own motion changes the direction from which light seems to come.

Formulas and arguments
CallFormulaWhat each argument means
doppler.longitudinal(f, v, approaching)
  • f — frequency the source sends out, Hz
  • v — speed, m/s (use lorentz.speed(0.8) for 0.8 c) of the source
  • approaching — True (default) if it moves towards you, False if away
Returns the frequency you receive, in Hz
doppler.wavelength(wavelength, v, approaching)
  • wavelength () — wavelength the source sends out, m
  • v — speed, m/s (use lorentz.speed(0.8) for 0.8 c) of the source
  • approaching — True (default) towards you, False away
Returns the wavelength you receive, in m
doppler.transverse(f, v)
  • f — frequency the source sends out, Hz
  • v — speed, m/s (use lorentz.speed(0.8) for 0.8 c) of the source across your line of sight
Returns the frequency you receive, in Hz (always lower)
doppler.redshift(v)
  • v — speed, m/s (use lorentz.speed(0.8) for 0.8 c) at which the source moves away
Returns the redshift (no unit)
doppler.speed_from_redshift(z)
  • z — the redshift
Returns the speed of recession, in m/s
doppler.aberration(angle, v)
  • angle () — direction the light comes from, measured from your direction of travel in the rest frame, degrees
  • v — speed, m/s (use lorentz.speed(0.8) for 0.8 c)
  • mode — "deg" (default) or "rad"
Returns the direction it seems to come from when you move, in the same unit
doppler.longitudinal(5e14, lorentz.speed(0.5)) # 8.66025e+14   (Hz, towards us)
doppler.longitudinal(5e14, lorentz.speed(0.5), approaching=False) # 2.88675e+14   (Hz, away)
doppler.wavelength(500e-9, lorentz.speed(0.5), approaching=False) # 8.66025e-07   (m)
doppler.transverse(5e14, lorentz.speed(0.6))  # 4e+14   (Hz)
doppler.redshift(lorentz.speed(0.5))          # 0.732051
doppler.speed_from_redshift(1) / constants.c  # 0.6
doppler.aberration(90, lorentz.speed(0.5))    # 60   (degrees)

gr.schwarzschild_radius(M) gr.photon_sphere(M) gr.isco(M) gr.hawking_temperature(M) gr.hawking_lifetime(M) gr.plot_clock_rate(M)

Black holes. Squeeze a mass inside its Schwarzschild radius, schwarzschild_radius(), and not even light can escape. Outside it, light can still orbit at photon_sphere(), and matter can orbit stably down to isco(). A black hole also glows very faintly with Hawking radiation: hawking_temperature() is how hot, and hawking_lifetime() how long it takes to evaporate. In every function you can give a mass in kilograms or the name of a body. plot_clock_rate() shows how clocks slow near the horizon.

Formulas and arguments
CallFormulaWhat each argument means
gr.schwarzschild_radius(M)
  • M — mass in kg, or a body name such as 'sun', 'earth', 'moon'
Returns the radius of the event horizon, in m ( N m²/kg²)
gr.photon_sphere(M)
  • M — mass in kg, or a body name such as 'sun', 'earth', 'moon'
Returns the radius where light can orbit, in m
gr.isco(M)
  • M — mass in kg, or a body name such as 'sun', 'earth', 'moon'
Returns the radius of the innermost stable circular orbit, in m
gr.hawking_temperature(M)
  • M — mass in kg, or a body name such as 'sun', 'earth', 'moon'
Returns the temperature of the black hole's radiation, in K
gr.hawking_lifetime(M)
  • M — mass in kg, or a body name such as 'sun', 'earth', 'moon'
Returns the time to evaporate, in s
gr.plot_clock_rate(M, title, save)
  • M — mass in kg or a body name (the shape of the curve is the same for every mass)
  • title — chart heading
  • save — file name to write the chart to
Returns draws clock rate against distance from the horizon
gr.schwarzschild_radius('sun')                # 2954.13   (m)
gr.schwarzschild_radius('earth')              # 0.00886981   (m: Earth squeezed into a marble)
gr.photon_sphere('sun')                       # 4431.19   (m)
gr.isco('sun')                                # 8862.38   (m)
gr.hawking_temperature('sun')                 # 6.16843e-08   (K)
gr.hawking_lifetime('sun') / 3.156e7          # 2.0972e+67   (years)
gr.plot_clock_rate()
gr.plot_clock_rate(): a clock ticks more slowly the closer it is to a black hole, and stops at the event horizon
gr.plot_clock_rate(): a clock ticks more slowly the closer it is to a black hole, and stops at the event horizon

gr.gravity_at(M, r) gr.time_dilation(t_far, M, r) gr.redshift(M, r) gr.height_shift(height) gr.gps_drift(altitude)

Gravity slows clocks. The closer a clock is to a mass, the more slowly it ticks. gravity_at() gives the rate compared with a clock far away, time_dilation() the time that passes, and redshift() how much the light climbing out loses. Near the ground the effect is tiny, and height_shift() is the simple version for a small height: it has been measured in a tower of 22 m. gps_drift() adds the two effects that matter for a GPS satellite: its clock gains from weaker gravity and loses from its speed, and the net gain of about 38 millionths of a second a day would put positions out by about 10 km a day if it were not corrected.

Formulas and arguments
CallFormulaWhat each argument means
gr.gravity_at(M, r)
  • M — mass in kg, or a body name such as 'sun', 'earth', 'moon'
  • r — distance from its centre, m (outside the horizon)
Returns how fast a clock there ticks compared with one far away (1 far away, 0 at the horizon)
gr.time_dilation(t_far, M, r)
  • t_far () — time that passes far away, s
  • M — mass in kg, or a body name such as 'sun', 'earth', 'moon'
  • r — distance from its centre, m
Returns the time that passes on the clock at , in s
gr.redshift(M, r)
  • M — mass in kg, or a body name such as 'sun', 'earth', 'moon'
  • r — distance from its centre where the light starts, m
Returns the redshift of the light when it reaches far away
gr.height_shift(height, g)
  • height () — height the light climbs, m
  • g — gravitational field strength, m/s² (default 9.80665)
Returns the fractional shift in frequency (no unit)
gr.gps_drift(altitude)
  • altitude () — height of the satellite above Earth, m (default 2.02×10⁷, a GPS satellite)
Returns the seconds a satellite clock gains on the ground each day (the first term is gravity, the second its speed)
gr.gravity_at('sun', 6.957e8)                 # 0.999998   (on the surface of the Sun)
gr.time_dilation(86400, 'earth', 6.371e6)     # 86400   (seconds on the ground for a day far away)
gr.redshift('sun', 6.957e8)                   # 2.12314e-06   (light leaving the Sun)
gr.height_shift(22.5)                         # 2.45506e-15   (light climbing a 22.5 m tower)
gr.gps_drift() * 1e6                          # 38.5125   (microseconds gained per day)

gr.light_deflection(M, b) gr.perihelion_shift(M, a, e)

Bending light and turning orbits. Mass curves space, so light that passes a star is bent: light_deflection() gives the angle in radians, and for the Sun it is 1.75 seconds of arc, the figure measured in the eclipse of 1919 that made Einstein famous. perihelion_shift() is the small extra turn that a planet's orbit makes each time round, in radians; for Mercury it explains the 43 seconds of arc a century that Newton's gravity could not.

Formulas and arguments
CallFormulaWhat each argument means
gr.light_deflection(M, b)
  • M — mass in kg, or a body name such as 'sun', 'earth', 'moon'
  • b — closest distance of the ray to its centre, m
Returns the angle the ray is bent through, in radians
gr.perihelion_shift(M, a, e)
  • M — mass in kg, or a body name such as 'sun', 'earth', 'moon' of the central body
  • a — semi-major axis of the orbit, m
  • e — eccentricity of the orbit (0 to 1)
Returns the extra turn of the orbit each time round, in radians
gr.light_deflection('sun', 6.957e8)           # 8.49253e-06   (radians, grazing the Sun)
math.degrees(gr.light_deflection('sun', 6.957e8)) * 3600 # 1.75171   (arcseconds)
gr.perihelion_shift('sun', 5.79e10, 0.2056) * 206264.806 * 415.2 # 42.9993   (Mercury, arcseconds per century)

Nuclear & Particle Physics

The structure of nuclei, nuclear reactions, radiation and dose, fission reactors, fusion, and the particles of the Standard Model. Every function is shown with its formula, and each table says which argument is which letter of it.

nucleus.semf_binding(Z, A) radiation.attenuation(I0, mu, x) reactor.fuel_use(power, time) fusion.reaction_energy(name) particle.check(reactants, products)

How the nuclear and particle functions work. These sit beside the nuclear namespace (binding energy and radioactive decay) and cover the rest of the subject in five namespaces: nucleus (the structure of nuclei, nuclear reactions, scattering), radiation (shielding and dose), reactor (fission and chain reactions), fusion (energy released and the barrier to overcome) and particle (the particles of the Standard Model, conservation laws, and collisions).

Units. Nuclear energies are in MeV and nuclear masses in atomic mass units (u), unless a function says otherwise. radiation is in SI units: grays, sieverts, metres. Particle masses are in MeV/c², and lifetimes in seconds. Every function below is shown with its formula, and each table says which argument is which letter of it.

nucleus.semf_binding(26, 56)                  # 495.384   (MeV, iron-56)
radiation.half_value_layer(0.5)               # 1.38629   (m)
fusion.reaction_energy('DT')                  # 17.5893   (MeV)
particle.mass('proton')                       # 938.272   (MeV/c²)

nucleus.semf_binding(Z, A) nucleus.semf_terms(Z, A) nucleus.semf_per_nucleon(Z, A) nucleus.stable_z(A) nucleus.is_magic(n) nucleus.plot_binding()

Binding energy from the liquid-drop model. Treat the nucleus as a drop of liquid, and its binding energy is the sum of five terms: a volume term (each nucleon is held by its neighbours), a surface term (those at the surface have fewer neighbours), the electric Coulomb repulsion of the protons, an asymmetry term (nuclei prefer equal numbers of protons and neutrons), and a pairing term. semf_binding() adds them up, and semf_terms() shows each one. It works well for heavy nuclei and badly for the lightest. stable_z() finds the number of protons that makes a nucleus with A nucleons most tightly bound, and is_magic() says whether a number fills a nuclear shell. plot_binding() draws the famous curve whose peak, near iron, explains why fusion of light nuclei and fission of heavy ones both release energy.

Formulas and arguments
CallFormulaWhat each argument means
nucleus.semf_binding(Z, A)
  • Z — number of protons
  • A — number of nucleons (protons + neutrons)
Returns the binding energy , in MeV (, , , MeV; is the pairing term, below)
nucleus.semf_terms(Z, A)
  • Z — number of protons
  • A — number of nucleons (protons + neutrons)
Returns a dict of the five terms (volume, surface, coulomb, asymmetry, pairing) in MeV, with MeV
nucleus.semf_per_nucleon(Z, A)
  • Z — number of protons
  • A — number of nucleons (protons + neutrons)
Returns the binding energy for each nucleon, in MeV
nucleus.stable_z(A)
  • A — number of nucleons
Returns the number of protons that makes a nucleus with nucleons most tightly bound
nucleus.is_magic(n)
  • n — number of protons or of neutrons
Returns True if is 2, 8, 20, 28, 50, 82 or 126, the numbers that fill a nuclear shell
nucleus.plot_binding(title, save)
  • title — chart heading
  • save — file name to write the chart to
Returns draws the binding energy for each nucleon against
nucleus.semf_binding(26, 56)                  # 495.384   (MeV, iron-56)
nucleus.semf_per_nucleon(26, 56)              # 8.84614   (MeV (the measured value is 8.79))
nucleus.semf_binding(92, 238)                 # 1814.73   (MeV, uranium-238)
nucleus.semf_terms(26, 56)                    # {'volume': 882, 'surface': -260.543, 'coulomb': -120.796, 'asymmetry': -6.77143, 'pairing': 1.49399}
nucleus.stable_z(56)                          # 26
nucleus.stable_z(238)                         # 92
nucleus.is_magic(82)                          # True   (lead)
nucleus.is_magic(83)                          # False
nucleus.plot_binding()
nucleus.plot_binding(): the binding energy for each nucleon peaks near iron, so fusing light nuclei or splitting heavy ones both release energy
nucleus.plot_binding(): the binding energy for each nucleon peaks near iron, so fusing light nuclei or splitting heavy ones both release energy

nucleus.q_value(reactants, products) nucleus.alpha_decay(parent_mass, daughter_mass) nucleus.beta_q(parent_mass, daughter_mass, kind) nucleus.daughter_number(N0, half_life_parent, half_life_daughter, t) nucleus.specific_activity(half_life, molar_mass)

Reactions and decays. Any nuclear reaction releases (or needs) the energy of the mass that goes missing: q_value() takes the atomic masses before and after, in u, and gives the energy in MeV. alpha_decay() gives both the energy released and the share that the alpha particle carries off (the lighter particle gets most of it). beta_q() does the same for beta decay; for beta-plus decay two electron masses must also be found. daughter_number() follows a decay chain: how many nuclei of the daughter there are at time t, as the parent decays into it and it decays in turn. specific_activity() is how radioactive a gram of a pure element is.

Formulas and arguments
CallFormulaWhat each argument means
nucleus.q_value(reactants, products)
  • reactants () — atomic masses before the reaction, u (one number or a list)
  • products () — atomic masses after, u
Returns the energy released, in MeV (negative: the reaction needs energy)
nucleus.alpha_decay(parent_mass, daughter_mass)
  • parent_mass () — atomic mass of the nucleus that decays, u
  • daughter_mass () — atomic mass of the nucleus left, u
Returns the pair in MeV: the energy released, and the part of it carried off by the alpha particle
nucleus.beta_q(parent_mass, daughter_mass, kind)
  • parent_mass () — atomic mass of the parent, u
  • daughter_mass () — atomic mass of the daughter, u
  • kind — "minus" (default) or "plus"
Returns the energy released by the beta decay, in MeV
nucleus.daughter_number(N0, half_life_parent, half_life_daughter, t)
  • N0 () — number of parent nuclei at the start
  • half_life_parent () — half-life of the parent, s (so )
  • half_life_daughter () — half-life of the daughter, s
  • t — time passed, s
Returns the number of daughter nuclei at time
nucleus.specific_activity(half_life, molar_mass)
  • half_life () — half-life, s
  • molar_mass () — molar mass, g/mol
Returns the decays each second in a gram, in Bq/g ( per mol)
nucleus.q_value([2.01410177812, 3.0160492779], [4.00260325413, 1.00866491595]) # 17.5893   (deuterium + tritium)
nucleus.alpha_decay(238.0507882, 234.0436012) # (4.26973, 4.19794)   (uranium-238 (Q, alpha energy), MeV)
nucleus.beta_q(14.003241989, 14.003074004)    # 0.156477   (carbon-14, MeV)
nucleus.beta_q(11.011433613, 11.009305406, 'plus') # 0.960414   (carbon-11, MeV)
nucleus.daughter_number(1e6, 10, 2, 5)        # 132583   (parent half-life 10 s, daughter 2 s, after 5 s)
nucleus.specific_activity(5.05e10, 226)       # 3.65743e+10   (radium-226, Bq/g (one curie))

nucleus.closest_approach(K, Z1, Z2) nucleus.scattering_angle(K, Z1, Z2, b) nucleus.coulomb_barrier(Z1, A1, Z2, A2)

Charged particles meet nuclei. A charged particle that is fired at a nucleus is pushed back by the electric force. closest_approach() is how near a particle of a given energy gets in a head-on approach, which sets an upper limit on the size of the nucleus: this is how Rutherford found it. scattering_angle() gives the angle a particle is turned through if it was aimed off-centre by the impact parameter b. coulomb_barrier() is the height of the electric barrier two nuclei must climb (or tunnel through) to touch.

Formulas and arguments
CallFormulaWhat each argument means
nucleus.closest_approach(K, Z1, Z2)
  • K — kinetic energy of the particle, MeV
  • Z1 () — charge of the particle, in units of
  • Z2 () — charge of the nucleus, in units of
Returns the distance it reaches in a head-on approach, in m
nucleus.scattering_angle(K, Z1, Z2, b)
  • K — kinetic energy of the particle, MeV
  • Z1 () — charge of the particle, in units of
  • Z2 () — charge of the nucleus, in units of
  • b — impact parameter (how far off-centre it was aimed), m
  • mode — "deg" (default) or "rad"
Returns the angle it is turned through (Rutherford scattering)
nucleus.coulomb_barrier(Z1, A1, Z2, A2)
  • Z1 () — protons in the first nucleus
  • A1 () — nucleons in the first nucleus
  • Z2 () — protons in the second nucleus
  • A2 () — nucleons in the second nucleus
Returns the height of the electric barrier between them, in MeV ( fm)
nucleus.closest_approach(5, 2, 79)            # 4.55029e-14   (5 MeV alpha at gold, m)
nucleus.scattering_angle(5, 2, 79, 1e-13)     # 25.6349   (degrees)
nucleus.scattering_angle(5, 2, 79, 1e-14)     # 132.546   (aimed closer: turned further)
nucleus.coulomb_barrier(1, 2, 1, 3)           # 0.444076   (deuterium and tritium, MeV)
nucleus.coulomb_barrier(1, 1, 1, 1)           # 0.599985   (two protons, MeV)

radiation.attenuation(I0, mu, x) radiation.half_value_layer(mu) radiation.thickness(fraction, mu) radiation.mean_free_path(mu) radiation.inverse_square(I1, r1, r2) radiation.plot_attenuation()

Shielding. A beam of gamma rays or X-rays that crosses a material loses intensity exponentially: attenuation() gives what is left after a thickness x, using the material's attenuation coefficient mu in 1/m. Each half_value_layer() of material halves the beam, and thickness() gives the thickness that cuts it down to any fraction. mean_free_path() is the average distance a photon goes before it interacts. Away from a small source the intensity also falls as inverse_square(). plot_attenuation() shows several materials side by side.

Formulas and arguments
CallFormulaWhat each argument means
radiation.attenuation(I0, mu, x)
  • I0 () — intensity of the beam going in
  • mu () — attenuation coefficient of the material, 1/m
  • x — thickness crossed, m
Returns the intensity left, in the same unit as
radiation.half_value_layer(mu)
  • mu () — attenuation coefficient, 1/m
Returns the thickness that halves the beam, in m
radiation.thickness(fraction, mu)
  • fraction () — fraction of the beam allowed through (0 to 1)
  • mu () — attenuation coefficient, 1/m
Returns the thickness needed, in m
radiation.mean_free_path(mu)
  • mu () — attenuation coefficient, 1/m
Returns the average distance before an interaction, in m
radiation.inverse_square(I1, r1, r2)
  • I1 () — intensity at distance
  • r1 () — first distance from a point source, m
  • r2 () — second distance, m
Returns the intensity at
radiation.plot_attenuation(mus, title, save)
  • mus () — attenuation coefficients to draw, 1/m
  • title — chart heading
  • save — file name to write the chart to
Returns draws the fraction of the beam that gets through against thickness, with half-value layers marked
radiation.attenuation(100, 0.5, 2)            # 36.7879   (a beam of 100 through 2 m)
radiation.half_value_layer(0.5)               # 1.38629   (m)
radiation.thickness(0.01, 0.5)                # 9.21034   (to cut the beam to 1%, m)
radiation.mean_free_path(0.5)                 # 2   (m)
radiation.inverse_square(100, 1, 3)           # 11.1111   (three times as far: a ninth)
radiation.plot_attenuation()
radiation.plot_attenuation(): the beam halves each half-value layer, and a denser material (larger mu) needs far less thickness
radiation.plot_attenuation(): the beam halves each half-value layer, and a denser material (larger mu) needs far less thickness

radiation.absorbed_dose(energy, mass) radiation.equivalent_dose(dose, radiation_type) radiation.effective_dose(*parts) radiation.gamma_dose_rate(activity, photon_energy, r, mu_en_rho, yield_)

Dose. The absorbed dose is the energy absorbed in each kilogram of tissue, in grays (absorbed_dose()). Different radiation does different harm for the same dose, so the equivalent dose in sieverts weights it: 1 for photons and electrons, 2 for protons, 10 for neutrons and 20 for alpha particles (equivalent_dose() takes a name or a number). Different tissues are also not equally sensitive, so the effective dose adds up the equivalent doses with a weight for each tissue. gamma_dose_rate() gives the dose rate in air at a distance from a point source of gamma rays, from its activity.

Formulas and arguments
CallFormulaWhat each argument means
radiation.absorbed_dose(energy, mass)
  • energy () — energy absorbed, J
  • mass () — mass that absorbs it, kg
Returns the absorbed dose , in grays (Gy)
radiation.equivalent_dose(dose, radiation_type)
  • dose () — absorbed dose, Gy
  • radiation_type () — the radiation weighting: a number, or 'photon' or 'electron' (1), 'proton' (2), 'neutron' (10), 'alpha' (20)
Returns the equivalent dose , in sieverts (Sv)
radiation.effective_dose(*parts)
  • parts () — pairs of the equivalent dose to a tissue (Sv) and that tissue's weighting factor
Returns the whole-body effective dose , in Sv
radiation.gamma_dose_rate(activity, photon_energy, r, mu_en_rho, yield_)
  • activity () — activity of the source, Bq
  • photon_energy () — energy of each photon, MeV
  • r — distance from the source, m
  • mu_en_rho () — energy-absorption coefficient of the absorber, m²/kg
  • yield_ () — photons per decay (default 1)
Returns the dose rate in Gy/s
radiation.absorbed_dose(0.35, 70)             # 0.005   (0.35 J in a 70 kg person, Gy)
radiation.equivalent_dose(0.01, 'alpha')      # 0.2   (Sv)
radiation.equivalent_dose(0.01, 'neutron')    # 0.1   (Sv)
radiation.effective_dose((0.01, 0.12), (0.02, 0.08)) # 0.0028   (two tissues, Sv)
radiation.gamma_dose_rate(3.7e10, 1.17, 1, 0.0027, 2) # 2.98045e-06   (one curie of cobalt-60 at 1 m, Gy/s)

reactor.fissions_per_second(power, energy) reactor.fuel_use(power, time, molar_mass, energy) reactor.energy_per_kg(molar_mass, energy)

Fission power. Each fission of uranium-235 releases about 200 MeV. fissions_per_second() says how many fissions give a thermal power, fuel_use() how much uranium that burns in a given time, and energy_per_kg() how much energy a kilogram yields: about 20 million times as much as a kilogram of coal.

Formulas and arguments
CallFormulaWhat each argument means
reactor.fissions_per_second(power, energy)
  • power () — thermal power, W
  • energy () — energy released by each fission, MeV (default 200)
Returns the number of fissions each second
reactor.fuel_use(power, time, molar_mass, energy)
  • power () — thermal power, W
  • time () — time, s
  • molar_mass () — molar mass of the fuel, g/mol (default 235)
  • energy () — energy for each fission, MeV (default 200)
Returns the mass of fuel fissioned, in kg
reactor.energy_per_kg(molar_mass, energy)
  • molar_mass () — molar mass of the fuel, g/mol (default 235)
  • energy () — energy for each fission, MeV (default 200)
Returns the energy from a kilogram of fuel, in J/kg
reactor.fissions_per_second(3e9)              # 9.36226e+19   (a 3 GW (thermal) station)
reactor.fuel_use(3e9, 86400)                  # 3.15654   (kg of uranium-235 in a day)
reactor.energy_per_kg()                       # 8.21152e+13   (J/kg)

reactor.eta(nu, sigma_f, sigma_a) reactor.four_factor(eta, f, p, epsilon) reactor.six_factor(...) reactor.multiplication(N0, k, generations) reactor.reactor_period(k, generation_time) reactor.moderator_loss(A) reactor.collisions_to_slow(A, E0, E)

The chain reaction. A reactor works when each fission leads on average to exactly one more: the multiplication factor k is 1. eta() is the neutrons made for each one absorbed by the fuel. four_factor() and six_factor() multiply the factors that decide k in an infinite and a finite reactor. multiplication() shows how the neutron count grows (or dies) over generations, and reactor_period() how quickly. The neutrons from fission are fast but fission is most likely with slow ones, so a moderator slows them: moderator_loss() is the average energy loss in a collision (largest for light nuclei) and collisions_to_slow() the number of collisions needed.

Formulas and arguments
CallFormulaWhat each argument means
reactor.eta(nu, sigma_f, sigma_a)
  • nu () — neutrons released by each fission
  • sigma_f () — fission cross-section
  • sigma_a () — absorption cross-section, in the same unit
Returns the neutrons made for each neutron absorbed by the fuel
reactor.four_factor(eta, f, p, epsilon)
  • eta () — neutrons made per neutron absorbed in the fuel
  • f — fraction of thermal neutrons absorbed by the fuel
  • p — chance a neutron escapes being captured while slowing
  • epsilon () — fast-fission factor
Returns the multiplication factor of an infinitely large reactor
reactor.six_factor(eta, f, p, epsilon, p_fnl, p_tnl)
  • eta () — neutrons made per neutron absorbed
  • f — thermal utilisation
  • p — resonance escape probability
  • epsilon () — fast-fission factor
  • p_fnl () — chance a fast neutron does not leak out (default 1)
  • p_tnl () — chance a thermal neutron does not leak out (default 1)
Returns the multiplication factor : above 1 the chain reaction grows, at 1 it holds steady, below 1 it dies away
reactor.multiplication(N0, k, generations)
  • N0 () — neutrons in the first generation
  • k — multiplication factor
  • generations () — number of generations
Returns the number of neutrons after generations
reactor.reactor_period(k, generation_time)
  • k — multiplication factor (above 1)
  • generation_time () — time for one generation of neutrons, s
Returns the time for the neutron population to grow by a factor , in s
reactor.moderator_loss(A)
  • A — mass number of the moderator nucleus
Returns the average logarithmic energy loss for each collision (1 for hydrogen)
reactor.collisions_to_slow(A, E0, E)
  • A — mass number of the moderator nucleus
  • E0 () — starting energy of the neutron
  • E — final energy, in the same unit
Returns the number of collisions needed
reactor.eta(2.43, 585, 681)                   # 2.08744   (uranium-235)
reactor.four_factor(1.65, 0.71, 0.87, 1.02)   # 1.03959
reactor.six_factor(1.65, 0.71, 0.87, 1.02, 0.97, 0.99) # 0.998317
reactor.multiplication(1000, 1.01, 100)       # 2704.81   (k slightly above 1)
reactor.multiplication(1000, 0.99, 100)       # 366.032   (slightly below)
reactor.reactor_period(1.001, 1e-4)           # 0.1   (s)
reactor.moderator_loss(1)                     # 1   (hydrogen)
reactor.collisions_to_slow(1, 2e6, 0.025)     # 18.1975   (from 2 MeV to thermal in water)
reactor.collisions_to_slow(12, 2e6, 0.025)    # 115.343   (in graphite)

fusion.reaction_energy(name) fusion.gamow_energy(Z1, Z2, reduced_mass) fusion.tunnel_probability(E, Z1, Z2, reduced_mass) fusion.power_density(n, sigma_v, energy) fusion.triple_product(n, T, tau) fusion.mass_loss_rate(power)

Fusion. reaction_energy() gives the energy released by the main fusion reactions: 'DT' (the easiest on Earth), 'DD-n', 'DD-p', 'DHe3', and 'pp', the Sun's chain. Two nuclei must get close against their electric repulsion, and gamow_energy() and tunnel_probability() estimate how hard that is: even at the Sun's core, nuclei get through by quantum tunnelling, not by brute force. power_density() is the fusion power in a hot gas, and triple_product() is the Lawson number that decides whether a plasma can sustain itself. mass_loss_rate() is the mass a star turns into light each second.

Formulas and arguments
CallFormulaWhat each argument means
fusion.reaction_energy(name)
  • name — the reaction: 'DT', 'DD-n', 'DD-p', 'DHe3', or 'pp' (four hydrogen into one helium)
Returns the energy released, in MeV
fusion.gamow_energy(Z1, Z2, reduced_mass)
  • Z1 () — protons in the first nucleus
  • Z2 () — protons in the second nucleus
  • reduced_mass () — reduced mass of the pair, , in u
Returns the Gamow energy , in MeV ()
fusion.tunnel_probability(E, Z1, Z2, reduced_mass)
  • E — energy of the collision, MeV
  • Z1 () — protons in the first nucleus
  • Z2 () — protons in the second nucleus
  • reduced_mass () — reduced mass of the pair, u
Returns the chance of tunnelling through the electric barrier
fusion.power_density(n, sigma_v, energy)
  • n — total ion density, 1/m³
  • sigma_v () — reaction rate parameter, m³/s
  • energy () — energy released by each reaction, MeV
Returns the fusion power in each cubic metre, in W/m³
fusion.triple_product(n, T, tau)
  • n — ion density, 1/m³
  • T — temperature, keV
  • tau () — energy confinement time, s
Returns the Lawson triple product, in keV s/m³ (deuterium–tritium needs about to ignite)
fusion.mass_loss_rate(power)
  • power () — power radiated, W
Returns the mass turned into energy each second, in kg/s
fusion.kev_to_kelvin(T)
  • T () — temperature in keV
Returns the temperature in kelvin
fusion.reaction_energy('DT')                  # 17.5893   (MeV)
fusion.reaction_energy('pp')                  # 26.731   (MeV)
fusion.reaction_energy('DHe3')                # 18.3531   (MeV)
fusion.gamow_energy(1, 1, 0.5)                # 0.489565   (two protons, MeV)
fusion.tunnel_probability(0.001, 1, 1, 0.5)   # 2.45894e-10   (at 1 keV)
fusion.power_density(1e20, 1e-22, 17.59)      # 704557   (W/m³)
fusion.triple_product(1e20, 15, 3)            # 4.5e+21   (keV s/m³)
fusion.mass_loss_rate(3.828e26)               # 4.25922e+09   (the Sun, kg/s)
fusion.kev_to_kelvin(10)                      # 1.16045e+08   (K)

particle.properties(name) particle.mass(name) particle.hadron(quarks) particle.check(reactants, products) particle.plot_masses()

Particles and conservation laws. particle.names lists the particles in the table, from the electron to the Higgs boson, and properties() and mass() give their mass, charge, baryon number, lepton number, and spin. hadron() builds a particle out of quarks: write them as letters, with a capital letter for an antiquark, so 'uud' is a proton and 'uD' a positive pion. check() tests a reaction against the rules every reaction must obey: charge, baryon number and lepton number must be the same before and after, and the energy must be there. It is how you find out why a decay does not happen. plot_masses() shows how widely the masses vary.

Formulas and arguments
CallFormulaWhat each argument means
particle.properties(name)
  • name — name of a particle, such as 'proton' (see particle.names)
Returns a dict: mass (MeV/c²), charge, baryon_number, lepton_number, spin
particle.mass(name)
  • name — name of a particle, such as 'proton' (see particle.names)
Returns the mass in MeV/c²
particle.hadron(quarks)
  • quarks — letters for the quarks (u, d, s, c, b, t); a capital letter is an antiquark, so 'uud' is a proton and 'uD' is a pion+
Returns a dict: charge (units of ), baryon_number, kind (baryon, meson, antibaryon), and name (such as 'proton', or None)
particle.check(reactants, products)
  • reactants — list of particle names before
  • products — list of particle names after
Returns a dict saying whether charge, baryon_number, lepton_number and energy (enough mass) are conserved, and allowed if all four are
particle.plot_masses(title, save)
  • title — chart heading
  • save — file name to write the chart to
Returns draws a bar chart of particle masses on a log scale
particle.names                                # ('electron', 'positron', 'muon', 'antimuon', 'tau', 'antitau', 'neutrino_e', 'antineutrino_e', 'neutrino_mu', 'antineutrino_mu', 'proton', 'antiproton', 'neutron', 'antineutron', 'pion+', 'pion-', 'pion0', 'kaon+', 'kaon-', 'kaon0', 'lambda', 'photon', 'W+', 'W-', 'Z', 'higgs', 'gluon')
particle.mass('proton')                       # 938.272   (MeV/c²)
particle.properties('neutron')                # {'mass': 939.565, 'charge': 0, 'baryon_number': 1, 'lepton_number': 0, 'spin': 0.5}
particle.hadron('uud')                        # {'charge': 1, 'baryon_number': 1, 'kind': 'baryon', 'name': 'proton'}
particle.hadron('uD')                         # {'charge': 1, 'baryon_number': 0, 'kind': 'meson', 'name': 'pion+'}
particle.hadron('uds')                        # {'charge': 0, 'baryon_number': 1, 'kind': 'baryon', 'name': 'lambda'}
particle.check(['neutron'], ['proton', 'electron', 'antineutrino_e']) # {'charge': True, 'baryon_number': True, 'lepton_number': True, 'energy': True, 'allowed': True}
particle.check(['neutron'], ['proton', 'electron']) # {'charge': True, 'baryon_number': True, 'lepton_number': False, 'energy': True, 'allowed': False}   (missing the neutrino)
particle.check(['proton'], ['positron', 'pion0']) # {'charge': True, 'baryon_number': False, 'lepton_number': False, 'energy': True, 'allowed': False}   (the proton does not decay)
particle.plot_masses()
particle.plot_masses(): the masses of the particles span more than five orders of magnitude
particle.plot_masses(): the masses of the particles span more than five orders of magnitude

particle.decay_distance(lifetime, energy, mass) particle.width(lifetime) particle.lifetime(width) particle.force_range(mass) particle.cm_energy_fixed_target(...) particle.cm_energy_collider(...) particle.bending_radius(p, B)

Decays and accelerators. decay_distance() is how far a fast unstable particle travels before it decays (time dilation makes it much further than c times its lifetime at rest). A particle that lives for a very short time has an uncertain energy, its width(), and lifetime() goes the other way. A force carried by a particle has a range of about ħ/(m c): force_range() shows why the electromagnetic force (massless photon) reaches far, the nuclear force (pion) reaches about a femtometre, and the weak force (W boson) is far shorter. cm_energy_fixed_target() and cm_energy_collider() give the energy actually available to make new particles, and show why colliders beat fixed targets. bending_radius() and momentum_from_radius() link a particle's momentum to the radius of its path in a magnetic field, which is how accelerators are sized.

Formulas and arguments
CallFormulaWhat each argument means
particle.decay_distance(lifetime, energy, mass)
  • lifetime () — lifetime of the particle at rest, s
  • energy () — total energy of the particle, MeV
  • mass () — its mass, MeV/c²
Returns the average distance it travels before decaying, in m
particle.width(lifetime)
  • lifetime () — lifetime, s
Returns the energy width , in MeV
particle.lifetime(width)
  • width () — energy width, MeV
Returns the lifetime , in s
particle.force_range(mass)
  • mass () — mass of the particle that carries the force, MeV/c²
Returns the range of the force, in m
particle.cm_energy_fixed_target(E_beam, m_beam, m_target)
  • E_beam () — total energy of the beam particle, MeV
  • m_beam () — its mass, MeV/c²
  • m_target () — mass of the target at rest, MeV/c²
Returns the energy available in the centre-of-mass frame, in MeV
particle.cm_energy_collider(E1, E2, m1, m2)
  • E1 () — total energy of the first beam, MeV
  • E2 () — total energy of the second beam, MeV
  • m1 () — mass of the first particle, MeV/c² (default 0)
  • m2 () — mass of the second, MeV/c² (default 0)
Returns the energy available when they collide head-on, in MeV
particle.bending_radius(p, B, charge)
  • p — momentum, GeV/c
  • B — magnetic field, T
  • charge () — charge in units of (default 1)
Returns the radius of its circle, in m
particle.momentum_from_radius(B, r, charge)
  • B — magnetic field, T
  • r — radius of the circle, m
  • charge () — charge in units of (default 1)
Returns the momentum , in GeV/c
particle.decay_distance(2.197e-6, 1e4, 105.658) # 62333.9   (a 10 GeV muon, m)
particle.decay_distance(2.6e-8, 1000, 139.57) # 55.3007   (a 1 GeV charged pion, m)
particle.width(1e-23)                         # 65.8212   (MeV)
particle.lifetime(2.0)                        # 3.29106e-22   (s for a width of 2 MeV)
particle.force_range(139.57)                  # 1.41382e-15   (the pion, m)
particle.force_range(80379)                   # 2.45496e-18   (the W boson, m)
particle.cm_energy_fixed_target(1e6 + 938.272, 938.272, 938.272) # 43359.7   (a 1 TeV proton on a proton at rest, MeV)
particle.cm_energy_collider(6.5e6, 6.5e6)     # 1.3e+07   (6.5 TeV beams colliding, MeV)
particle.bending_radius(7000, 8.33)           # 2803.06   (a 7 TeV proton in 8.33 T, m)
particle.momentum_from_radius(1, 1)           # 0.299792   (GeV/c)

particle.barns(x) particle.mean_free_path(n, sigma) particle.interaction_probability(n, sigma, thickness) particle.reaction_rate(luminosity, sigma)

Cross-sections and rates. How likely a particle is to interact is measured by its cross-section, the effective target area of a nucleus, in square metres, or in barns (1 barn = 10−28 m²; barns() converts). With n targets in each cubic metre, mean_free_path() is how far a particle goes before it interacts and interaction_probability() the chance it interacts in a given thickness. At a collider, reaction_rate() gives the number of events each second from the luminosity and the cross-section.

Formulas and arguments
CallFormulaWhat each argument means
particle.barns(x)
  • x — cross-section in barns
Returns the cross-section in m²
particle.mean_free_path(n, sigma)
  • n — targets per cubic metre
  • sigma () — cross-section of each target, m²
Returns the average distance before an interaction, in m
particle.interaction_probability(n, sigma, thickness)
  • n — targets per cubic metre
  • sigma () — cross-section of each target, m²
  • thickness () — thickness of the material, m
Returns the chance that a particle interacts
particle.reaction_rate(luminosity, sigma)
  • luminosity () — collider luminosity, 1/(m² s)
  • sigma () — cross-section of the process, m²
Returns the number of events each second
particle.barns(1)                             # 1e-28   (m²)
particle.mean_free_path(8.5e28, particle.barns(1)) # 0.117647   (m)
particle.interaction_probability(1e28, particle.barns(1), 0.1) # 0.0951626
particle.reaction_rate(1e30, 1e-37)           # 1e-07   (events/s)

Units

One function covers ten categories — length, mass, volume, area, time, speed, data, energy, pressure, and temperature. convert() figures out the category from the units themselves, so you never specify it; it just refuses to convert across categories (length to mass, say) with a clear error.

convert(value, from_unit, to_unit)

Converts between units within the same category. Case-insensitive, with generous aliases (km / kilometer / kilometers). Raises a clear error for incompatible categories.

convert(10, "km", "miles")
# 6.213711922373339

convert(98.6, "F", "C")
# 37.0

convert(1, "gb", "mb")
# 1024.0

Casting & Rounding

Type conversion and rounding, fixing two classic Python gotchas.

cast(value, target)

Casts a value — or an entire DataFrame column — to another type. Fixes Python's bool("False") gotcha (which is actually True!) and parses string literals like "[1, 2, 3]" into real lists instead of character-splitting.

cast("123", "int")
# 123

cast("False", "bool")
# False   (Python's own bool("False") is True!)

cast("[1, 2, 3]", "list")
# [1, 2, 3]

c["age"] = cast(c["age"], "int")   # casts a whole column

round(value, places=0)

Rounds using standard round-half-up — the way Excel works — instead of Python's "banker's rounding," where round(2.5) is actually 2.

round(2.5)
# 3   (Python's own round(2.5) gives 2)

round(1250, -2)
# 1300   (Python's own round(1250, -2) gives 1200)

Data Operations

The typical flow: load() a file into a variable, inspect it with loadh()/loadc()/loads(), optionally clean() it, then query it with sql() — real SQL, powered by DuckDB, run directly against your DataFrame by its variable name.

load(path)

Loads a CSV, Excel file, or Google Sheet into a DataFrame and previews it.

c = load("data.csv")
c = load("data.xlsx")
c = load("https://docs.google.com/spreadsheets/d/.../edit")

loadh(df) · loadc(df) · loads(df)

Preview just the column headers, the record count, or the schema (columns + datatypes).

loadh(c)
loadc(c)
loads(c)

sql(query)

Runs any SQL query against DataFrames already loaded in your script, by variable name. Powered by DuckDB — full SQL, including JOINs and window functions.

d = sql("select * from c where age > 30 order by age desc")

xtract(df, filename=None)

Exports a DataFrame to a CSV file in the current directory.

xtract(d)
xtract(sql("select * from c"), filename="filtered.csv")

clean(df)

Standardizes column names, strips whitespace from text columns, drops exact duplicates, and drops fully-empty rows. Returns a new DataFrame.

c = clean(c)

Lookups & Excel Functions

If you already know these from spreadsheets, they behave the same way here. Every function auto-detects your loaded dataset, so once you've called load(), you don't need to pass it in again — unless you've loaded more than one dataset, in which case pass it explicitly as df=.

vlook(value, return_col, df=None, default=...)

Searches the FIRST column (like Excel's VLOOKUP) and returns the matching row's value. Auto-detects your loaded dataset.

c = load("employees.csv")
vlook("Alice", "salary")
# 85000

vlook("Zoe", "salary", default=0)
# 0  (no error, since a fallback was given)

xlook(value, lookup_col, return_col, df=None, default=...)

Searches ANY column you specify (like Excel's XLOOKUP) — not limited to the first.

xlook("Engineering", "department", "name")
# "Alice"

countif(col, value, df=None)

Counts rows where col equals value.

countif("department", "Engineering")
# 3

sumif(col, value, sum_col, df=None)

Sums sum_col for rows where col equals value.

sumif("department", "Engineering", "salary")
# 250000

avgif(col, value, avg_col, df=None)

Averages avg_col for rows where col equals value.

avgif("department", "Engineering", "salary")
# 83333.33

AI

Bring your own API key from Anthropic, OpenAI, Gemini, or any OpenAI-compatible provider (Groq, Mistral, etc.) — the provider is detected automatically from the key's format, so you rarely need to specify it. Nothing here works without a key; everything else in Fyrefly does.

ask(data, question, api_key=None, provider=None, model=None, base_url=None)

Converts a plain-English question into SQL and runs it. Provider auto-detected from your key's format (Anthropic, OpenAI, Gemini, or Groq/any OpenAI-compatible API). Pass a tuple of DataFrames for joins across any number of datasets. Transient provider errors retry automatically.

c = load("sales.csv")

ask(c, "total sales by region last quarter?", api_key="AIza...")

ask((c, d), "top customers by total sales?", api_key="...")

insights(df, api_key=None, provider=None, model=None, base_url=None)

Generates a short, plain-English summary of a dataset — trends, ranges, and data-quality issues.

insights(c, api_key="AIza...")

Visualization

No API key, no separate charting library to learn. Call viz() directly with a kind, or use the shorter viz.bar(), viz.scatter(), etc. once you've loaded a single dataset — both produce the same chart.

viz(df, x=None, y=None, kind="auto", hue=None, title=None, figsize=(8,5), save=None)

Visualizes a DataFrame. kind="auto" picks a sensible chart based on your columns' types — heatmap, bar, histogram, or scatter.

viz(c)                               # auto -> correlation heatmap
viz(c, "age", "salary")              # auto -> scatter
viz(c, "department", "salary", kind="box")

viz.bar / .scatter / .line / .box / .violin / .hist / .kde / .pie / .heatmap / .pairplot

Namespace shortcuts for each chart type. Auto-detects your loaded DataFrame when exactly one is in scope.

c = load("data.csv")

viz.bar("department", "salary")
viz.scatter("age", "salary")
viz.heatmap()
viz.pairplot(hue="department")

API Testing

One function, api.gui(), opens a free Postman-style window in your browser, where you build a request, send it, and read a plain-English answer. Nothing here needs an account or a key from Fyrefly. Treat API keys like passwords: keep them in environment variables or in a secret variable, not in files you share.

api.gui(url=None, method="get", send=None, headers=None, token=None, ..., port=0, open_browser=True, block=None, folder=None)

A free Postman-style tool for everyday API work. It opens a window in your browser. Pick the method, fill in the address, use the Params, Authorization, Headers, Body, Variables, Checks, Save to variable and Settings tabs, then press Send (or Ctrl+Enter). The body can be JSON, form-data with file upload, a plain web form, or text/XML, and every row has an on/off checkbox. Keep several requests open in tabs, and save any of them (Ctrl+S). Every argument is optional: pass an address, a method, send=, headers= or token= and the window opens already filled in.

Collections and environments. Save requests into collections and folders shown in the sidebar, reopen, copy, rename or delete them, and run a whole collection in order with a pass/fail table. An environment (UAT, SIT, Production...) is a named set of variables, switched from a dropdown and used as {{HOST}} anywhere. Built in: {{$guid}}, {{$timestamp}}, {{$isoTimestamp}}, {{$randomInt}}.

Secrets stay secret. Mark a variable as secret and it is hidden on screen and kept in a separate secrets.json (owner-only, listed in .gitignore). A token typed straight into a request is not written to the saved file, so collections are safe to share. Save to variable: log in once, and the token is copied from the answer into {{TOKEN}} for every later request. Click "to var" next to any value in the Preview tab, or set it up in the Save to variable tab.

Bring your work over. Import a Postman collection (v2.x), a Postman environment, or a pasted curl command, and export a collection as Postman JSON.

The answer. Status, time and size, then these tabs: Report (a plain-English report; when something fails, a HOW TO FIX IT list says what is missing and which tab to fix it in), Body (colored JSON, a find box with match count, Copy, and Download), Preview (a collapsible JSON tree, HTML in a locked box, or an image), Table, Headers (with any redirects followed), and Request (what was sent, and the same call as curl, secrets hidden). Under the address bar, a line shows what {{variables}} will resolve to, and flags any that have no value.

Settings. Check SSL certificates on or off (for internal servers), follow redirects, a proxy, a default timeout, cookies that persist between calls, and a history of what you sent (tokens are never kept).

Where your work lives. In one plain folder, ~/fyrefly_api. It is created for you the first time, and created again if you delete it (the saved work in it is gone). folder="..." picks another place and folder=False saves nothing.

~/fyrefly_api/
  collections/Shop API.json      requests, safe to share or commit
  environments/UAT.json          variables; secret values are blank here
  secrets.json                   the secret values (not shared, in .gitignore)
  history.json  settings.json

Safe by design. It runs only on your computer: the page is served from 127.0.0.1, needs a one-time secret that is part of the address it prints, rejects other websites and other host names, and loads nothing from the internet. Tokens and passwords go into hidden boxes, are sent only to the API you are testing, and are never shown in a report. In a script it waits until you click Quit or press Ctrl+C. In a notebook it returns at once, and window.stop() closes it. Pass open_browser=False to only print the address, or port=8765 to choose the port.

from fyrefly import api

api.gui()                                                     # opens the window
api.gui("https://api.example.com/users", method="post",
        send="name=Asha, age=30", token="YOUR_TOKEN")         # opens it already filled in

# Fyrefly API window is running. It is only reachable from this computer.
# Saved requests and environments live in: /home/you/fyrefly_api
# If it does not open by itself, paste this into your browser:
#   http://127.0.0.1:52817/?t=...
# Click Quit in the window, or press Ctrl+C here, to close it.