add(*nos)
Adds any count of numbers.
add(5, 4, 6, 10)
# 25
add(1, 2, 3, 4, 5)
# 15
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.
Fyrefly is a single package — every function in this reference comes from one install.
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.
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
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.
| Call | Formula | What each argument means |
|---|---|---|
si(p, n, r) |
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".
| Call | Formula | What each argument means |
|---|---|---|
ci(p, n, r) |
Returns the compound interest earned (the principal is not included) | |
ci(p, n, r, "semi-annual") |
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.
| Call | Formula | What each argument means |
|---|---|---|
si_ci_diff(p, r, n=2) |
Returns the difference between compound and simple interest | |
si_ci_diff(p, r, n=3) |
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.
| Call | Formula | What each argument means |
|---|---|---|
profit(sp, cp) |
Returns the profit | |
profit(sp, cp, "%") |
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.
| Call | Formula | What each argument means |
|---|---|---|
loss(sp, cp) |
Returns the loss | |
loss(sp, cp, "%") |
Returns the pair (loss, loss percentage) |
loss(80, 100)
# 20
loss(80, 100, "%")
# (20, 20.0)
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.
| Call | Formula | What each argument means |
|---|---|---|
ap.tn(a, d, n) |
Returns the th term | |
ap.sn(a, d, n) |
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.
| Call | Formula | What each argument means |
|---|---|---|
gp.tn(a, r, n) |
Returns the th term | |
gp.sn(a, r, n) |
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.
| Call | Formula | What each argument means |
|---|---|---|
hp.tn(a, d, n) |
Returns the th term | |
hp.sn(a, d, n) |
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
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.
| Call | Formula | What each argument means |
|---|---|---|
avg(*nos) |
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."
| Call | Formula | What each argument means |
|---|---|---|
avgc(*pairs) |
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.
| Call | Formula | What each argument means |
|---|---|---|
mode(*nos) |
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.
| Call | Formula | What each argument means |
|---|---|---|
median(*nos) |
Returns the median |
median(1, 3, 2)
# 2
median(1, 2, 3, 4)
# 2.5
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.
| Call | Formula | What each argument means |
|---|---|---|
hyp(a, b) |
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.
| Call | Formula | What each argument means |
|---|---|---|
slope(point1, point2) |
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.
| Call | Formula | What each argument means |
|---|---|---|
centroid(point1, point2, point3) |
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)²).
| Call | Formula | What each argument means |
|---|---|---|
dist(point1, point2) |
Returns the distance |
dist((6, 2), (10, 5))
# 5.0
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.
| Call | Formula | What each argument means |
|---|---|---|
circle.area(r) |
Returns the area | |
circle.perimeter(r) |
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.
| Call | Formula | What each argument means |
|---|---|---|
semicircle.area(r) |
Returns the area | |
semicircle.perimeter(r) |
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.
| Call | Formula | What each argument means |
|---|---|---|
square.area(s) |
Returns the area | |
square.perimeter(s) |
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.
| Call | Formula | What each argument means |
|---|---|---|
rectangle.area(l, w) |
Returns the area | |
rectangle.perimeter(l, w) |
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.
| Call | Formula | What each argument means |
|---|---|---|
triangle.area(base, height) |
Returns the area | |
triangle.area_sides(a, b, c) |
Returns the area (Heron's formula); an error if they cannot form a triangle | |
triangle.perimeter(a, b, c) |
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.
| Call | Formula | What each argument means |
|---|---|---|
equilateral_triangle.area(s) |
Returns the area | |
equilateral_triangle.perimeter(s) |
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.
| Call | Formula | What each argument means |
|---|---|---|
parallelogram.area(base, height) |
Returns the area | |
parallelogram.perimeter(a, b) |
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.
| Call | Formula | What each argument means |
|---|---|---|
rhombus.area(d1, d2) |
Returns the area | |
rhombus.perimeter(s) |
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.
| Call | Formula | What each argument means |
|---|---|---|
kite.area(d1, d2) |
Returns the area | |
kite.perimeter(a, b) |
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.
| Call | Formula | What each argument means |
|---|---|---|
trapezoid.area(a, b, height) |
Returns the area | |
trapezoid.perimeter(a, b, c, d) |
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.
| Call | Formula | What each argument means |
|---|---|---|
ellipse.area(a, b) |
Returns the area | |
ellipse.perimeter(a, b) |
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.
| Call | Formula | What each argument means |
|---|---|---|
polygon.area(n_sides, length) |
Returns the area of a regular polygon | |
polygon.perimeter(n_sides, length) |
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.
| Call | Formula | What each argument means |
|---|---|---|
sector.area(r, angle) |
Returns the area of the sector | |
sector.arc_length(r, angle) |
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.
| Call | Formula | What each argument means |
|---|---|---|
annulus.area(outer_r, inner_r) |
Returns the area of the ring |
annulus.area(10, 6)
# 201.06192982974676
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.
| Call | Formula | What each argument means |
|---|---|---|
cube.volume(s) |
Returns the volume | |
cube.lsa(s) |
Returns the four side faces | |
cube.tsa(s) |
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.
| Call | Formula | What each argument means |
|---|---|---|
cuboid.volume(l, w, h) |
Returns the volume | |
cuboid.lsa(l, w, h) |
Returns the four side faces | |
cuboid.tsa(l, w, h) |
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.
| Call | Formula | What each argument means |
|---|---|---|
sphere.volume(r) |
Returns the volume | |
sphere.surface_area(r) |
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.
| Call | Formula | What each argument means |
|---|---|---|
hemisphere.volume(r) |
Returns the volume | |
hemisphere.csa(r) |
Returns the curved dome only | |
hemisphere.tsa(r) |
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.
| Call | Formula | What each argument means |
|---|---|---|
cylinder.volume(r, h) |
Returns the volume | |
cylinder.csa(r, h) |
Returns the side only | |
cylinder.tsa(r, h) |
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.
| Call | Formula | What each argument means |
|---|---|---|
cone.volume(r, h) |
Returns the volume | |
cone.csa(r, h) |
Returns the curved surface | |
cone.tsa(r, h) |
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.
| Call | Formula | What each argument means |
|---|---|---|
frustum.volume(r1, r2, h) |
Returns the volume | |
frustum.csa(r1, r2, h) |
Returns the curved surface | |
frustum.tsa(r1, r2, h) |
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.
| Call | Formula | What each argument means |
|---|---|---|
pyramid.volume(s, h) |
Returns the volume | |
pyramid.lsa(s, h) |
Returns the four triangular faces | |
pyramid.tsa(s, h) |
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.
| Call | Formula | What each argument means |
|---|---|---|
prism.volume(base_area, h) |
Returns the volume | |
prism.lsa(base_perimeter, h) |
Returns the side faces | |
prism.tsa(base_area, base_perimeter, h) |
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.
| Call | Formula | What each argument means |
|---|---|---|
torus.volume(R, r) |
Returns the volume | |
torus.surface_area(R, r) |
Returns the surface area |
torus.volume(5, 2)
# 394.78417604357435
torus.surface_area(5, 2)
# 394.78417604357435
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.
| Call | Formula | What each argument means |
|---|---|---|
eqn.l(*equations) |
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).
| Call | Formula | What each argument means |
|---|---|---|
eqn.q(*coeffs) |
Returns the roots | |
eqn.q.sum(*coeffs) |
Returns the sum of the roots (Vieta's formulas) | |
eqn.q.mul(*coeffs) |
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
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.
| Call | Formula | What each argument means |
|---|---|---|
log(x, base) |
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.
| Call | Formula | What each argument means |
|---|---|---|
exp(a, b) |
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.
| Call | Formula | What each argument means |
|---|---|---|
sqrt(no) |
Returns the square root | |
curt(no) |
Returns the cube root (negative numbers allowed) | |
nroot(no, n) |
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
Degrees by default — pass mode="rad" for radians.
sin(angle, mode="deg") · cos(...) · tan(...)The three primary trigonometric ratios.
| Call | Formula | What each argument means |
|---|---|---|
sin(angle) |
Returns the sine | |
cos(angle) |
Returns the cosine | |
tan(angle) |
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.
| Call | Formula | What each argument means |
|---|---|---|
cosec(angle) |
Returns undefined where | |
sec(angle) |
Returns undefined where | |
cot(angle) |
Returns undefined where |
cosec(30)
# 2.0
sec(60)
# 2.0
cot(45)
# 1.0
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.
| Call | Formula | What each argument means |
|---|---|---|
newton.force(m, a) |
Returns the force , in N | |
newton.mass(F, a) |
Returns the mass , in kg | |
newton.accel(F, m) |
Returns the acceleration , in m/s² | |
newton.second(F=, m=, a=) |
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) |
Returns the weight , in N | |
newton.net(*forces) |
Returns the net force, a number or an pair | |
newton.equilibrium(*forces, tol) |
Returns True if the forces cancel (Newton's first law) | |
newton.third(force) |
Returns the equal and opposite force | |
newton.components(force, angle, mode) |
Returns the pair | |
newton.resultant(f1, f2, angle=90) |
Returns the size of the combined force | |
newton.lift(m, a, g=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.
| Call | Formula | What each argument means |
|---|---|---|
motion.suvat(s=, u=, v=, a=, t=) |
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) |
Returns the final speed | |
motion.distance(u, a, t) |
Returns the distance | |
motion.time(u, v, a) |
Returns the time | |
motion.accel(u, v, t) |
Returns the acceleration | |
motion.speed(distance, time) |
Returns the average speed | |
motion.fall(height, g=9.81) |
Returns the pair (time, speed at the bottom) | |
motion.stopping_distance(speed, decel) |
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.
| Call | Formula | What each argument means |
|---|---|---|
projectile.time(v, angle, h0, g) |
Returns the time of flight , in s (from the ground, ) | |
projectile.range(v, angle, h0, g) |
Returns the horizontal distance to landing, in m (from the ground, ) | |
projectile.height(v, angle, h0, g) |
Returns the greatest height , in m | |
projectile.at(t, v, angle, h0, g) |
Returns the position | |
projectile.velocity_at(t, v, angle, g) |
Returns the velocity | |
projectile.best_angle(v, h0, g) |
Returns the angle in degrees that gives the longest range | |
projectile.path(v, angle, h0, g, steps) |
Returns a table with columns t, x, y | |
projectile.plot(v, angle, h0, g, title, save, compare) |
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])

![compare=[15, 30, 45, 60, 75]: 30° and 60° land together, 45° goes farthest](images/mech-projectile-angles.png)
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.
| Call | Formula | What each argument means |
|---|---|---|
gravity.force(m1, m2, r) |
Returns the pull , in N ( N m²/kg²) | |
gravity.law(F=, m1=, m2=, r=) |
Returns the variable you leave out | |
gravity.field(M, r) |
Returns the field strength , in m/s² | |
gravity.orbit_speed(M, r) |
Returns the circular-orbit speed | |
gravity.escape_speed(M, r) |
Returns the escape speed | |
gravity.orbital_period(M, a) |
Returns the period , in s | |
gravity.kepler3(T1, a1, a2) |
Returns the new period (same time unit as ) | |
gravity.weight_on(body, m) |
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.
| Call | Formula | What each argument means |
|---|---|---|
energy.kinetic(m, v) |
Returns the kinetic energy, in J | |
energy.potential(m, h, g=9.81) |
Returns the potential energy, in J | |
energy.spring(k, x) |
Returns the stored energy, in J | |
energy.work(force, distance, angle=0) |
Returns the work , in J | |
energy.power(work, time) |
Returns the power , in W | |
energy.power_fv(force, speed) |
Returns the power , in W | |
energy.speed_after_drop(height, v0, g) |
Returns the speed at the bottom | |
energy.height_reached(speed, g) |
Returns the greatest height | |
energy.efficiency(useful, supplied) |
Returns the efficiency , as a percentage | |
energy.kinetic_law(KE=, m=, v=) |
Returns the variable you leave out | |
energy.potential_law(PE=, m=, g=, h=) |
Returns the variable you leave out | |
energy.spring_law(E=, k=, x=) |
Returns the variable you leave out | |
energy.power_law(P=, W=, t=) |
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.recoilMomentum 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.
| Call | Formula | What each argument means |
|---|---|---|
momentum.p(m, v) |
Returns the momentum , in kg m/s | |
momentum.impulse(force, time) |
Returns the impulse , in N s | |
momentum.impulse_change(m, u, v) |
Returns the impulse | |
momentum.collide(m1, u1, m2, u2, e) |
Returns the speeds afterwards | |
momentum.elastic(m1, u1, m2, u2) |
Returns the speeds afterwards () | |
momentum.inelastic(m1, u1, m2, u2) |
Returns the one common speed (, they stick) | |
momentum.energy_lost(m1, u1, m2, u2, e) |
Returns the kinetic energy lost, in J | |
momentum.recoil(m_gun, m_bullet, v_bullet) |
Returns the recoil speed of the gun (negative: backwards) | |
momentum.com(masses, positions) |
Returns the centre of mass | |
momentum.rocket(exhaust_speed, m_start, m_end) |
Returns the speed gained | |
momentum.law(p=, m=, v=) |
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.
| Call | Formula | What each argument means |
|---|---|---|
circular.accel(v, r) |
Returns the centripetal acceleration , in m/s² | |
circular.force(m, v, r) |
Returns the inward force , in N | |
circular.speed(r, period) |
Returns the speed | |
circular.period(r, v) |
Returns the time for one lap | |
circular.omega(period) |
Returns the angular speed , in rad/s | |
circular.banked_speed(r, angle, mu, g) |
Returns the fastest steady speed on the curve | |
circular.vertical_circle(r, g) |
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).
| Call | Formula | What each argument means |
|---|---|---|
friction.force(mu, normal) |
Returns the friction force , in N | |
friction.incline(angle, mu, g) |
Returns the acceleration down the slope (0 if friction holds it) | |
friction.repose(mu) |
Returns the steepest angle in degrees at which it stays put | |
friction.stopping_distance(v, mu, g) |
Returns the sliding distance | |
friction.slide_time(length, angle, mu, g) |
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.
| Call | Formula | What each argument means |
|---|---|---|
inertia.point(m, r) |
Returns the moment of inertia , in kg m² | |
inertia.ring(m, r) |
Returns about its axis | |
inertia.disc(m, r) |
Returns of a solid disc or cylinder | |
inertia.sphere(m, r) |
Returns of a solid sphere about a diameter | |
inertia.hollow_sphere(m, r) |
Returns of a thin shell | |
inertia.rod_center(m, length) |
Returns about the middle of the rod | |
inertia.rod_end(m, length) |
Returns about one end | |
inertia.plate(m, a, b) |
Returns about the centre, axis perpendicular to the plate | |
inertia.shift(I_cm, m, d) |
Returns about the new axis (parallel axis theorem) | |
rotation.torque(force, arm, angle=90) |
Returns the torque , in N m | |
rotation.angular_accel(torque, I) |
Returns the angular acceleration , in rad/s² | |
rotation.torque_law(torque=, I=, alpha=) |
Returns the variable you leave out | |
rotation.angular_momentum(I, omega) |
Returns the angular momentum | |
rotation.kinetic(I, omega) |
Returns the rotational kinetic energy, in J | |
rotation.rolling_speed(height, shape, g) |
Returns the speed at the bottom | |
rotation.lever(load, load_arm, effort_arm) |
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.
| Call | Formula | What each argument means |
|---|---|---|
shm.spring_period(m, k) |
Returns the period , in s | |
shm.pendulum_period(length, g) |
Returns the period | |
shm.hooke(k, x) |
Returns the restoring force (negative: it pulls back) | |
shm.frequency(period) |
Returns the frequency , in Hz | |
shm.omega(period) |
Returns the angular frequency , in rad/s | |
shm.position(A, period, t, phase) |
Returns the displacement | |
shm.velocity(A, period, t, phase) |
Returns the velocity | |
shm.max_speed(A, period) |
Returns the greatest speed, at the middle | |
shm.energy(k, A) |
Returns the total energy , in J | |
shm.path(A, period, cycles, steps) / shm.plot(...) |
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)

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.
| Call | Formula | What each argument means |
|---|---|---|
damping.omega0(m, k) |
Returns the natural frequency , rad/s | |
damping.gamma(m, c) |
Returns the decay rate , 1/s | |
damping.damping_ratio(m, k, c) |
Returns the damping ratio | |
damping.critical(m, k) |
Returns the critical damping constant, N s/m | |
damping.regime(m, k, c) |
Returns "underdamped", "critically damped", "overdamped", or "undamped" | |
damping.omega_damped(m, k, c) |
Returns the ringing frequency (underdamped only) | |
damping.period(m, k, c) |
Returns the time of one ring | |
damping.q_factor(m, k, c) |
Returns the quality factor | |
damping.decay_time(m, c) |
Returns the time for the amplitude to fall to | |
damping.amplitude(A, m, c, t) |
Returns the envelope size at time | |
damping.position(x0, m, k, c, t, v0) |
Returns the displacement at time (underdamped, critical, overdamped, in that order) | |
damping.driven_amplitude(F0, m, k, c, w) |
Returns the steady swing size , in m | |
damping.resonance_frequency(m, k, c) |
Returns the driving frequency of the biggest swing (0 if there is no peak) | |
damping.peak_amplitude(F0, m, k, c) |
Returns the biggest swing | |
damping.path(x0, m, k, c, v0, duration, steps) / plot(...) / plot_resonance(F0, m, k, c) |
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)


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).
| Call | Formula | What each argument means |
|---|---|---|
orbit.periapsis(a, e) |
Returns the closest distance | |
orbit.apoapsis(a, e) |
Returns the farthest distance | |
orbit.semi_major(r_peri, r_apo) |
Returns the semi-major axis | |
orbit.eccentricity(r_peri, r_apo) |
Returns the eccentricity | |
orbit.period(M, a) |
Returns the period , in s | |
orbit.speed(M, r, a) |
Returns the speed there (vis-viva) | |
orbit.total_energy(M, m, a) |
Returns the total energy , in J (negative: bound) | |
orbit.areal_velocity(M, a, e) |
Returns the area swept per second (Kepler's second law) | |
orbit.hohmann(M, r1, r2) |
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) |
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)

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.
| Call | Formula | What each argument means |
|---|---|---|
pulley.atwood(m1, m2, g) |
Returns (acceleration, tension); positive means goes down | |
pulley.on_incline(m_slope, angle, m_hang, mu, g) |
Returns (acceleration, tension); positive means the hanging mass goes down, 0 if friction holds | |
pulley.on_table(m_table, m_hang, mu, g) |
Returns (acceleration, tension) | |
pulley.lift_tension(m, a, g) |
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.
| Call | Formula | What each argument means |
|---|---|---|
fluid.pressure(F, A) |
Returns the pressure , in Pa | |
fluid.pressure_law(P=, F=, A=) |
Returns the variable you leave out | |
fluid.hydrostatic(depth, rho, g) |
Returns the extra pressure | |
fluid.absolute_pressure(depth, rho, g) |
Returns the total pressure ( Pa) | |
fluid.hydraulic(F1, A1, A2) |
Returns the force on the second piston | |
fluid.buoyancy(volume, rho, g) |
Returns the upthrust , in N | |
fluid.floats(rho_object, rho_fluid) |
Returns True if it floats | |
fluid.submerged_fraction(rho_object, rho_fluid) |
Returns the fraction under the surface (1 if it sinks) | |
fluid.continuity(A1, v1, A2) |
Returns the speed in the second pipe | |
fluid.flow_rate(A, v) |
Returns the volume per second , in m³/s | |
fluid.bernoulli(P1, v1, h1, v2, h2, rho, g) |
Returns the pressure at the second point | |
fluid.torricelli(depth, g) |
Returns the jet speed | |
fluid.reynolds(rho, v, length, viscosity) |
Returns the Reynolds number (smooth below ~2000, turbulent above ~4000) | |
fluid.terminal_velocity(radius, rho_sphere, rho_fluid, viscosity, g) |
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.
| Call | Formula | What each argument means |
|---|---|---|
elasticity.stress(F, A) |
Returns the stress , in Pa | |
elasticity.strain(extension, length) |
Returns the strain (no unit) | |
elasticity.youngs(F, A, L, dL) |
Returns Young's modulus , in Pa | |
elasticity.youngs_law(E=, F=, A=, L=, dL=) |
Returns the variable you leave out | |
elasticity.extension(F, L, A, E) |
Returns the extension , in m | |
elasticity.elastic_energy(F, extension) |
Returns the stored energy , in J | |
elasticity.series(*ks) |
Returns the combined of springs end to end | |
elasticity.parallel(*ks) |
Returns the combined of springs side by side | |
elasticity.shear_modulus(F, A, shear_strain) |
Returns the shear modulus | |
elasticity.bulk_modulus(dP, dV, V) |
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
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.hHow 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.
| Call | Formula | What each argument means |
|---|---|---|
photon.energy(wavelength=, frequency=, unit, length) |
Returns the energy of one photon ( J s) | |
photon.frequency(wavelength=, energy=, unit, length) |
Returns the frequency , in Hz | |
photon.wavelength(energy=, frequency=, unit, length) |
Returns the wavelength | |
photon.momentum(wavelength, length) |
Returns the photon's momentum , in kg m/s | |
photon.photoelectric(work_function, wavelength=, frequency=, unit, length) |
Returns the greatest kinetic energy of the electrons; 0 if no electrons come out | |
photon.threshold_frequency(work_function, unit) |
Returns the lowest frequency that frees electrons | |
photon.threshold_wavelength(work_function, unit, length) |
Returns the longest wavelength that frees electrons | |
photon.stopping_voltage(work_function, wavelength=, frequency=, unit, length) |
Returns the voltage that just stops the fastest electrons, in volts | |
photon.compton_wavelength(mass, length) |
Returns the Compton wavelength | |
photon.compton_shift(angle, mass, length) |
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.
| Call | Formula | What each argument means |
|---|---|---|
photon.wien(temperature, length) |
Returns the wavelength of the brightest glow ( m K) | |
photon.stefan(temperature, area, emissivity) |
Returns the power radiated, in W ( W/(m² K⁴)) | |
photon.planck(wavelength, temperature, length) |
Returns the brightness at that wavelength, in W/(m² m sr) ( J/K) | |
photon.plot_blackbody(temperatures, title, save) |
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](images/q-blackbody.png)
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.
| Call | Formula | What each argument means |
|---|---|---|
matter.de_broglie(mass, speed, length) |
Returns the wavelength | |
matter.de_broglie_energy(mass, energy, unit, length) |
Returns the wavelength | |
matter.electron_wavelength(volts, length) |
Returns the wavelength of an electron | |
matter.momentum_spread(position_spread, length) |
Returns the smallest possible , in kg m/s () | |
matter.position_spread(momentum_spread, length) |
Returns the smallest possible | |
matter.speed_spread(mass, position_spread, length) |
Returns the smallest possible , in m/s | |
matter.energy_spread(lifetime) |
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.
| Call | Formula | What each argument means |
|---|---|---|
bohr.energy(n, Z, unit) |
Returns the energy of level (negative: bound) | |
bohr.radius(n, Z, length) |
Returns the orbit radius ( pm) | |
bohr.speed(n, Z) |
Returns the electron's speed , in m/s () | |
bohr.ionization(n, Z, unit) |
Returns the energy to free the electron | |
bohr.transition(n_from, n_to, Z, unit) |
Returns the photon's energy (positive when given out, negative when taken in) | |
bohr.wavelength(n_from, n_to, Z, length) |
Returns the wavelength of the light (Rydberg formula, m⁻¹) | |
bohr.series_name(n_to) |
Returns "Lyman", "Balmer", "Paschen", "Brackett", "Pfund" | |
bohr.spectrum(n_to, lines, Z, length, unit) |
Returns a table of n_from, n_to, wavelength and energy | |
bohr.plot_levels(Z, n_max, jumps, title, save) |
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()

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.
| Call | Formula | What each argument means |
|---|---|---|
box.energy(n, mass, width, unit, length) |
Returns the energy of level | |
box.levels(mass, width, count, unit, length) |
Returns a list of the lowest levels | |
box.wavefunction(n, width, x, length) |
Returns the value of the wave function | |
box.probability(n, width, x1, x2, length) |
Returns the chance of finding the particle between and | |
box.transition(n_from, n_to, mass, width, unit, length) |
Returns the photon's energy | |
box.photon_wavelength(n_from, n_to, mass, width, length) |
Returns the wavelength of the light | |
box.plot(width, levels, mass, length, title, save) |
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)

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.
| Call | Formula | What each argument means |
|---|---|---|
qho.omega(mass, k) |
Returns the angular frequency , rad/s | |
qho.energy(n, omega, unit) |
Returns the energy of level | |
qho.zero_point(omega, unit) |
Returns the lowest possible energy | |
qho.spacing(omega, unit) |
Returns the gap between neighbouring levels | |
qho.turning_point(n, mass, omega, length) |
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.
| Call | Formula | What each argument means |
|---|---|---|
tunnel.decay_constant(mass, barrier, energy, unit, length) |
Returns how fast the wave dies away in the barrier, | |
tunnel.approx(mass, barrier, energy, width, unit, length) |
Returns the chance of getting through a thick barrier | |
tunnel.probability(mass, barrier, energy, width, unit, length) |
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.
| Call | Formula | What each argument means |
|---|---|---|
nuclear.radius(mass_number, length) |
Returns the nuclear radius ( fm) | |
nuclear.mass_energy(mass, unit) |
Returns the energy locked up in that mass | |
nuclear.u_to_energy(mass_u, unit) |
Returns the equivalent energy | |
nuclear.mass_defect(Z, A, atomic_mass) |
Returns the missing mass , in u | |
nuclear.binding_energy(Z, A, atomic_mass, unit) |
Returns the binding energy | |
nuclear.binding_per_nucleon(Z, A, atomic_mass, unit) |
Returns the binding energy for each particle | |
nuclear.decay_constant(half_life) |
Returns the decay constant , per that time unit | |
nuclear.half_life(decay_constant) |
Returns the half-life | |
nuclear.mean_life(half_life) |
Returns the mean life | |
nuclear.remaining(N0, t, half_life) |
Returns the number left | |
nuclear.decay_law(N=, N0=, lam=, t=) |
Returns the variable you leave out | |
nuclear.activity(N, half_life) |
Returns decays per unit time (becquerels if the half-life is in seconds) | |
nuclear.age(fraction_left, half_life) |
Returns the age (carbon dating) | |
nuclear.time_to_fraction(fraction_left, half_life) |
Returns the time until that fraction is left | |
nuclear.path(N0, half_life, half_lives, steps) / plot(...) |
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)")

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.
| Call | Formula | What each argument means |
|---|---|---|
xray.bragg_angle(spacing, wavelength, order, length, mode) |
Returns the glancing angle | |
xray.bragg_spacing(angle, wavelength, order, length, mode) |
Returns the plane spacing | |
xray.min_wavelength(volts, length) |
Returns the shortest wavelength the tube makes | |
xray.max_energy(volts, unit) |
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
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.
| Call | Formula | What each argument means |
|---|---|---|
gas.law(P=, V=, n=, T=) |
Returns the missing one of , , , ( J/(mol K)) | |
gas.pressure(n, T, V) |
Returns the pressure , in Pa | |
gas.volume(n, T, P) |
Returns the volume , in m³ | |
gas.moles(P, V, T) |
Returns the amount of gas , in mol | |
gas.temperature(P, V, n) |
Returns the temperature , in K | |
gas.combined(P1=, V1=, T1=, P2=, V2=, T2=) |
Returns the missing one of the six | |
gas.density(P, T, molar_mass) |
Returns the density , in kg/m³ | |
gas.van_der_waals(n, V, T, a, b) |
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.
| Call | Formula | What each argument means |
|---|---|---|
gas.isothermal_work(n, T, V1, V2) |
Returns the work done BY the gas, in J (negative when squeezed) | |
gas.isobaric_work(P, V1, V2) |
Returns the work done BY the gas, in J | |
gas.adiabatic_pressure(P1, V1, V2, gamma) |
Returns the final pressure , in Pa | |
gas.adiabatic_temperature(T1, V1, V2, gamma) |
Returns the final temperature , in K | |
gas.adiabatic_work(P1, V1, P2, V2, gamma) |
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.
| Call | Formula | What each argument means |
|---|---|---|
gas.cv(dof) |
Returns the molar heat capacity at constant volume, in J/(mol K) | |
gas.cp(dof) |
Returns the molar heat capacity at constant pressure, in J/(mol K) | |
gas.gamma(dof) |
Returns the ratio (no unit) | |
gas.internal_energy(n, T, dof) |
Returns the energy stored in the gas, in J | |
gas.kinetic_energy(T, unit) |
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).
| Call | Formula | What each argument means |
|---|---|---|
gas.rms_speed(T, molar_mass) |
Returns the root-mean-square speed, in m/s | |
gas.mean_speed(T, molar_mass) |
Returns the average speed, in m/s | |
gas.most_probable_speed(T, molar_mass) |
Returns the speed most molecules have, in m/s | |
gas.maxwell(v, T, molar_mass) |
Returns the share of molecules per unit of speed, , in s/m | |
gas.speed_fraction(v1, v2, T, molar_mass) |
Returns the fraction of molecules with speeds between and (0 to 1) | |
gas.mean_free_path(T, P, diameter) |
Returns the average distance between collisions , in m | |
gas.plot_maxwell(molar_mass, temperatures, title, save) |
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))

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.
| Call | Formula | What each argument means |
|---|---|---|
heat.energy(m, c, dT) |
Returns the heat needed, in J | |
heat.law(Q=, m=, c=, dT=) |
Returns the missing one of the four | |
heat.latent(m, L) |
Returns the heat needed, in J | |
heat.mix(m1, c1, T1, m2, c2, T2) |
Returns the final temperature , in the scale you used | |
heat.to_kelvin(value, scale) |
Returns the temperature in kelvin | |
heat.from_kelvin(kelvin, scale) |
Returns the temperature in that scale | |
heat.linear_expansion(L0, alpha, dT) |
Returns the change in length (same unit as ) | |
heat.volume_expansion(V0, beta, dT) |
Returns the change in volume | |
heat.newton_cooling(T0, T_env, k, t) |
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.
| Call | Formula | What each argument means |
|---|---|---|
heat.conduction(k, A, dT, L) |
Returns the rate of heat flow , in W | |
heat.resistance(L, k, A) |
Returns the thermal resistance, in K/W | |
heat.convection(h, A, dT) |
Returns the rate of heat flow , in W | |
heat.radiation(emissivity, area, T, T_surroundings) |
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.
| Call | Formula | What each argument means |
|---|---|---|
entropy.heat(Q, T) |
Returns the entropy change , in J/K | |
entropy.phase_change(m, L, T) |
Returns the entropy change , in J/K | |
entropy.heating(m, c, T1, T2) |
Returns the entropy change , in J/K | |
entropy.ideal_gas(n, T1, T2, V1, V2, dof) |
Returns the entropy change of the gas, in J/K | |
entropy.universe(Q, T_hot, T_cold) |
Returns the total entropy change, which is never negative | |
entropy.helmholtz(U, T, S) |
Returns the Helmholtz free energy , in J | |
entropy.gibbs(H, T, S) |
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.
| Call | Formula | What each argument means |
|---|---|---|
engine.first_law(dU=, Q=, W=) |
Returns the missing one of , , | |
engine.efficiency(work, heat_in) |
Returns the efficiency (0 to 1) | |
engine.carnot(T_hot, T_cold) |
Returns the greatest efficiency any engine can have between these temperatures | |
engine.refrigerator_cop(T_hot, T_cold) |
Returns the best heat removed per unit of work | |
engine.heat_pump_cop(T_hot, T_cold) |
Returns the best heat delivered per unit of work | |
engine.max_work(Q_hot, T_hot, T_cold) |
Returns the most work possible, in J | |
engine.exhaust(Q_hot, work) |
Returns the heat thrown away to the cold side, in J | |
engine.otto(compression_ratio, gamma) |
Returns the efficiency of the ideal petrol-engine cycle | |
engine.diesel(compression_ratio, cutoff_ratio, gamma) |
Returns the efficiency of the ideal diesel cycle | |
engine.plot_carnot(T_hot, T_cold, V1, V2, n, gamma, title, save) |
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)

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=.
| Call | Formula | What each argument means |
|---|---|---|
stat.boltzmann_factor(energy, T, unit) |
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) |
Returns the ratio of the populations | |
stat.partition(energies, T, degeneracies, unit) |
Returns the partition function (energies measured from the lowest level) | |
stat.probabilities(energies, T, degeneracies, unit) |
Returns a list: the chance of finding the system in each level | |
stat.average_energy(energies, T, degeneracies, unit) |
Returns the average energy over the levels | |
stat.equipartition(T, dof, unit) |
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.
| Call | Formula | What each argument means |
|---|---|---|
stat.entropy(W) |
Returns the entropy , in J/K | |
stat.gibbs_entropy(probabilities) |
Returns the entropy , in J/K | |
stat.microstates(N, n) |
Returns the number of ways | |
stat.einstein_solid(N, q) |
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.
| Call | Formula | What each argument means |
|---|---|---|
stat.fermi_dirac(energy, mu, T, unit) |
Returns the chance the state is filled, between 0 and 1 | |
stat.bose_einstein(energy, mu, T, unit) |
Returns the average number of bosons in the state | |
stat.fermi_energy(density, mass, unit) |
Returns the Fermi energy | |
stat.fermi_temperature(density, mass) |
Returns the Fermi temperature , in K | |
stat.plot_occupancy(mu, T, unit, title, save) |
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.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).
| Call | Formula | What each argument means |
|---|---|---|
stat.einstein_heat_capacity(T, theta) |
Returns the molar heat capacity, in J/(mol K) (approaches when hot) | |
stat.debye_heat_capacity(T, theta) |
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
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.
| Call | Formula | What each argument means |
|---|---|---|
pipe.poiseuille(radius, length, pressure_drop, viscosity) |
Returns the flow rate , in m³/s (valid for smooth, laminar flow) | |
pipe.pressure_drop(radius, length, flow_rate, viscosity) |
Returns the pressure difference needed, in Pa | |
pipe.resistance(radius, length, viscosity) |
Returns the hydraulic resistance , in Pa s/m³ (so that ) | |
pipe.regime(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.
| Call | Formula | What each argument means |
|---|---|---|
pipe.friction_factor(reynolds, diameter, roughness) |
Returns the Darcy friction factor (no unit) | |
pipe.darcy_weisbach(f, length, diameter, rho, v) |
Returns the pressure lost to friction , in Pa | |
pipe.head_loss(f, length, diameter, v, g) |
Returns the friction loss as a height of liquid , in m | |
pipe.pump_power(flow_rate, pressure_rise, efficiency) |
Returns the power the pump needs, in W | |
pipe.plot_friction(diameter, roughness, title, save) |
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.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.
| Call | Formula | What each argument means |
|---|---|---|
pipe.venturi(A1, A2, pressure_drop, rho) |
Returns the flow rate , in m³/s | |
pipe.pitot(pressure_difference, rho) |
Returns the speed of the flow, in m/s | |
pipe.orifice(Cd, area, pressure_drop, rho) |
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.
| Call | Formula | What each argument means |
|---|---|---|
drag.dynamic_pressure(rho, v) |
Returns the dynamic pressure , in Pa | |
drag.force(Cd, rho, area, v) |
Returns the drag force , in N | |
drag.lift(Cl, rho, area, v) |
Returns the lift force , in N | |
drag.power(Cd, rho, area, v) |
Returns the power needed to overcome drag, in W | |
drag.stokes(radius, speed, viscosity) |
Returns the drag force , in N | |
drag.stall_speed(mass, Cl_max, rho, area, g) |
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().
| Call | Formula | What each argument means |
|---|---|---|
drag.terminal_speed(mass, Cd, rho, area, g) |
Returns the terminal speed , where drag equals weight, in m/s | |
drag.fall_speed(mass, Cd, rho, area, t, g) |
Returns the speed after time , in m/s | |
drag.fall_distance(mass, Cd, rho, area, t, g) |
Returns the distance fallen , in m | |
drag.plot_fall(mass, Cd, rho, area, duration, title, save) |
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.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.
| Call | Formula | What each argument means |
|---|---|---|
drag.sound_speed(T, gamma, molar_mass) |
Returns the speed of sound , in m/s | |
drag.mach(v, T, gamma, molar_mass) |
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.
| Call | Formula | What each argument means |
|---|---|---|
surface.capillary_rise(surface_tension, radius, rho, contact_angle, g) |
Returns the height the liquid climbs, in m (negative if it is pushed down, as mercury is) | |
surface.laplace(surface_tension, radius, surfaces) |
Returns the extra pressure inside, in Pa | |
surface.energy(surface_tension, area) |
Returns the energy stored in the surface, in J | |
surface.film_force(surface_tension, length, surfaces) |
Returns the pull of the film, in N | |
surface.capillary_length(surface_tension, rho, g) |
Returns the length below which surface tension beats gravity, in m | |
surface.young(gamma_solid_vapour, gamma_solid_liquid, gamma_liquid_vapour) |
Returns the contact angle of a drop on the surface, in degrees | |
surface.weber(rho, v, length, surface_tension) |
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).
| Call | Formula | What each argument means |
|---|---|---|
hydro.dam_force(width, depth, rho, g) |
Returns the total push of the liquid on the wall, in N | |
hydro.center_of_pressure(depth) |
Returns how far below the surface the push acts, in m | |
hydro.dam_moment(width, depth, rho, g) |
Returns the turning effect about the foot of the wall, in N m | |
hydro.gate_force(width, height, top_depth, rho, g) |
Returns the push on the gate, in N | |
hydro.manometer(height, rho, g) |
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).
| Call | Formula | What each argument means |
|---|---|---|
hydro.barometric(height, P0, T, molar_mass, g) |
Returns the air pressure at that height, in Pa | |
hydro.spinning_surface(omega, r, g) |
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) |
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.
| Call | Formula | What each argument means |
|---|---|---|
channel.manning(n, hydraulic_radius, slope) |
Returns the speed of the flow, in m/s | |
channel.froude(v, depth, g) |
Returns the Froude number (no unit) | |
channel.wave_speed(depth, g) |
Returns the speed of a shallow-water wave, in m/s | |
channel.regime(froude) |
Returns the word "subcritical", "critical" or "supercritical" | |
channel.critical_depth(flow_per_width, g) |
Returns the depth at which the flow is critical, in m | |
channel.weir(Cd, width, head, g) |
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)
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.
| Call | Formula | What each argument means |
|---|---|---|
charge.coulomb(q1, q2, r) |
Returns the force , in N: positive pushes them apart, negative pulls them together ( N m²/C²) | |
charge.field(q, r) |
Returns the electric field , in N/C (pointing away if is positive) | |
charge.potential(q, r) |
Returns the electric potential , in volts | |
charge.potential_energy(q1, q2, r) |
Returns the energy stored by the pair, in J (negative for opposite charges) | |
charge.force_on(q, E) |
Returns the force on the charge, in N | |
charge.acceleration(q, m, E) |
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.
| Call | Formula | What each argument means |
|---|---|---|
charge.flux(q_enclosed) |
Returns the electric flux out of the surface, in N m²/C (Gauss's law) | |
charge.field_line(linear_density, r) |
Returns the electric field , in N/C | |
charge.field_sheet(surface_density) |
Returns the electric field near a large flat sheet, in N/C | |
charge.field_conductor(surface_density) |
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.
| Call | Formula | What each argument means |
|---|---|---|
charge.dipole_moment(q, separation) |
Returns the dipole moment , in C m | |
charge.dipole_torque(p, E, angle) |
Returns the torque , in N m | |
charge.field_at(point, charges) |
Returns the field as a pair , in N/C | |
charge.potential_at(point, charges) |
Returns the potential , in volts (potentials add as plain numbers) | |
charge.plot_field(charges, size, title, save) |
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)])

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.
| Call | Formula | What each argument means |
|---|---|---|
capacitor.law(Q=, C=, V=) |
Returns the missing one of , , | |
capacitor.parallel_plate(area, gap, k) |
Returns the capacitance , in F | |
capacitor.sphere(radius) |
Returns the capacitance of the isolated sphere, in F | |
capacitor.energy(C, V) |
Returns the energy stored, in J | |
capacitor.series(*caps) |
Returns the combined capacitance, in F | |
capacitor.parallel(*caps) |
Returns the combined capacitance, in F | |
capacitor.field_between(V, gap) |
Returns the electric field between them, in V/m | |
capacitor.energy_density(E, k) |
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.
| Call | Formula | What each argument means |
|---|---|---|
circuit.ohm(V=, I=, R=) |
Returns the missing one of , , | |
circuit.power_law(P=, V=, I=) |
Returns the missing one of , , | |
circuit.resistance(rho, length, area) |
Returns the resistance , in ohms | |
circuit.series(*resistors) |
Returns the total resistance, in ohms | |
circuit.parallel(*resistors) |
Returns the total resistance, in ohms | |
circuit.divider(V, R1, R2) |
Returns the voltage across , in volts | |
circuit.current_divider(I, R1, R2) |
Returns the current through , in A | |
circuit.power_resistor(I, R) |
Returns the heat power , in W | |
circuit.energy(P, t) |
Returns the energy used, in J (divide by 3.6×10⁶ for kilowatt-hours) | |
circuit.terminal_voltage(emf, internal_resistance, R) |
Returns the voltage across the load, in volts | |
circuit.temperature_resistance(R0, alpha, dT) |
Returns the new resistance , in ohms | |
circuit.drift_velocity(I, n, area, q) |
Returns the drift speed , in m/s | |
circuit.battery_life(capacity_Ah, I) |
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.
| Call | Formula | What each argument means |
|---|---|---|
circuit.time_constant(R, C) |
Returns the time constant , in s | |
circuit.rc_charge(V, R, C, t) |
Returns the voltage across the capacitor, in volts | |
circuit.rc_discharge(V0, R, C, t) |
Returns the voltage across the capacitor, in volts | |
circuit.rl_current(V, R, L, t) |
Returns the current in the circuit, in A | |
circuit.plot_rc(V, R, C, title, save) |
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)

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.
| Call | Formula | What each argument means |
|---|---|---|
magnet.force_on_charge(q, v, B, angle) |
Returns the force , in N (always at right angles to the motion) | |
magnet.cyclotron_radius(m, v, q, B) |
Returns the radius of its circle, in m | |
magnet.cyclotron_frequency(q, B, m) |
Returns the number of turns per second , in Hz (it does not depend on the speed) | |
magnet.helix_pitch(m, v_parallel, q, B) |
Returns the distance moved along the field in one turn, in m | |
magnet.velocity_selector(E, B) |
Returns the speed that passes straight through, in m/s | |
magnet.force_on_wire(B, I, length, angle) |
Returns the force on the wire, in N | |
magnet.plot_helix(m, v, q, B, angle, turns, title, save) |
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.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.
| Call | Formula | What each argument means |
|---|---|---|
magnet.field_wire(I, r) |
Returns the magnetic field , in T ( T m/A) | |
magnet.field_loop(I, radius, turns) |
Returns the field at the centre of the coil, in T | |
magnet.field_solenoid(n, I) |
Returns the field inside a long solenoid, in T | |
magnet.field_toroid(N, I, r) |
Returns the field inside the toroid, in T | |
magnet.force_between_wires(I1, I2, length, d) |
Returns the force between them, in N: positive attracts, negative repels | |
magnet.magnetic_moment(N, I, area) |
Returns the magnetic moment , in A m² | |
magnet.torque_on_loop(N, I, area, B, angle) |
Returns the torque , in N m | |
magnet.hall_voltage(I, B, n, thickness, q) |
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.
| Call | Formula | What each argument means |
|---|---|---|
induction.flux(B, area, angle) |
Returns the magnetic flux , in Wb | |
induction.emf(N, dPhi, dt) |
Returns the induced emf , in volts (the minus sign is Lenz's law) | |
induction.motional_emf(B, length, v) |
Returns the emf across the rod, in volts | |
induction.generator_peak(N, B, area, omega) |
Returns the greatest emf , in volts | |
induction.generator(N, B, area, omega, t) |
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.
| Call | Formula | What each argument means |
|---|---|---|
induction.inductance_solenoid(N, area, length, mu_r) |
Returns the inductance , in H | |
induction.self_emf(L, dI, dt) |
Returns the emf across the inductor, in volts | |
induction.inductor_energy(L, I) |
Returns the energy stored, in J | |
induction.field_energy_density(B, mu_r) |
Returns the energy per cubic metre , in J/m³ | |
induction.transformer(V_primary, N_primary, N_secondary) |
Returns the voltage out of the secondary, in volts | |
induction.transformer_current(I_primary, N_primary, N_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.
| Call | Formula | What each argument means |
|---|---|---|
ac.rms(peak) |
Returns the rms value, the steady value with the same heating effect | |
ac.peak(rms) |
Returns the peak value | |
ac.reactance_c(C, f) |
Returns the capacitive reactance , in ohms (falls as the frequency rises) | |
ac.reactance_l(L, f) |
Returns the inductive reactance , in ohms (rises with the frequency) | |
ac.impedance(R, L, C, f) |
Returns the impedance , in ohms | |
ac.phase(R, L, C, f) |
Returns the angle by which the current lags the voltage (negative: it leads) | |
ac.current_rms(V_rms, R, L, C, f) |
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.
| Call | Formula | What each argument means |
|---|---|---|
ac.resonance(L, C) |
Returns the resonant frequency , in Hz | |
ac.q_factor(R, L, C) |
Returns the quality factor (no unit): a larger is a sharper resonance | |
ac.power(V_rms, I_rms, phase) |
Returns the average power , in W | |
ac.plot_resonance(L, C, resistances, V_rms, title, save) |
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))

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.
| Call | Formula | What each argument means |
|---|---|---|
emwave.wavelength(frequency) |
Returns the wavelength , in m ( m/s) | |
emwave.frequency(wavelength) |
Returns the frequency , in Hz | |
emwave.band(wavelength=, frequency=) |
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) |
Returns the magnetic field of the wave, in T | |
emwave.e_from_b(B) |
Returns the electric field of the wave, in V/m | |
emwave.intensity(E0) |
Returns the intensity (power per square metre), in W/m² | |
emwave.peak_field(intensity) |
Returns the peak electric field , in V/m | |
emwave.intensity_spread(power, r) |
Returns the intensity at that distance, in W/m² | |
emwave.radiation_pressure(intensity, reflect) |
Returns the pressure of the light, in Pa | |
emwave.speed_in_medium(n) |
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)
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.
| Call | Formula | What each argument means |
|---|---|---|
lorentz.speed(beta) |
Returns the speed , in m/s ( m/s) | |
lorentz.beta(v) |
Returns the speed as a fraction of the speed of light | |
lorentz.gamma(v) |
Returns the Lorentz factor (no unit): 1 at rest, and growing without limit as approaches | |
lorentz.speed_from_gamma(gamma) |
Returns the speed that has that Lorentz factor, in m/s | |
lorentz.time_dilation(proper_time, v) |
Returns the time between those ticks as measured by you, in s (longer than ) | |
lorentz.length_contraction(proper_length, v) |
Returns the length you measure, in m (shorter than ) | |
lorentz.velocity_addition(u, v) |
Returns the speed of the object relative to you, in m/s (never more than ) | |
lorentz.rapidity(v) |
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.
| Call | Formula | What each argument means |
|---|---|---|
lorentz.transform(x, t, v) |
Returns the pair : where and when the same event happens in the moving frame | |
lorentz.interval(dt, dx) |
Returns the spacetime interval , in m²: the same for every observer (positive: one event can cause the other) | |
lorentz.proper_time(dt, dx) |
Returns the time on a clock that travels from one event to the other, in s | |
lorentz.simultaneity(dx, v) |
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.
| Call | Formula | What each argument means |
|---|---|---|
lorentz.twin_trip(distance, v) |
Returns the pair for the round trip, in s: the traveller ages less | |
lorentz.muon_fraction(distance, v, lifetime) |
Returns the fraction still alive when they arrive | |
lorentz.plot_gamma(title, save) |
Returns draws the Lorentz factor against speed | |
lorentz.plot_minkowski(v, title, save) |
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()


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.
| Call | Formula | What each argument means |
|---|---|---|
relmech.momentum(m, v) |
Returns the relativistic momentum , in kg m/s | |
relmech.rest_energy(m) |
Returns the energy of the mass at rest, in J | |
relmech.energy_total(m, v) |
Returns the total energy , in J | |
relmech.kinetic(m, v) |
Returns the kinetic energy , in J (close to when ) | |
relmech.energy_momentum(p, m) |
Returns the total energy , in J | |
relmech.speed_from_momentum(p, m) |
Returns the speed , in m/s | |
relmech.speed_from_kinetic(m, K) |
Returns the speed , in m/s | |
relmech.accelerate(m, q, V) |
Returns the speed it reaches from rest, in m/s | |
relmech.plot_kinetic(m, title, save) |
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.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.
| Call | Formula | What each argument means |
|---|---|---|
relmech.invariant_mass(E, p) |
Returns the mass of the system, in kg: the same for every observer | |
relmech.mass_from_energy(E) |
Returns the mass that corresponds to it, in kg | |
relmech.decay_energy(M, m1, m2) |
Returns the energy released, in J (shared as motion of the products) | |
relmech.threshold_energy(m_beam, m_target, m_products) |
Returns the least kinetic energy the beam needs for the reaction to happen, in J | |
relmech.cherenkov_angle(n, v) |
Returns the half-angle of the cone of light, in degrees | |
relmech.cherenkov_threshold(n) |
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.
| Call | Formula | What each argument means |
|---|---|---|
doppler.longitudinal(f, v, approaching) |
Returns the frequency you receive, in Hz | |
doppler.wavelength(wavelength, v, approaching) |
Returns the wavelength you receive, in m | |
doppler.transverse(f, v) |
Returns the frequency you receive, in Hz (always lower) | |
doppler.redshift(v) |
Returns the redshift (no unit) | |
doppler.speed_from_redshift(z) |
Returns the speed of recession, in m/s | |
doppler.aberration(angle, v) |
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.
| Call | Formula | What each argument means |
|---|---|---|
gr.schwarzschild_radius(M) |
Returns the radius of the event horizon, in m ( N m²/kg²) | |
gr.photon_sphere(M) |
Returns the radius where light can orbit, in m | |
gr.isco(M) |
Returns the radius of the innermost stable circular orbit, in m | |
gr.hawking_temperature(M) |
Returns the temperature of the black hole's radiation, in K | |
gr.hawking_lifetime(M) |
Returns the time to evaporate, in s | |
gr.plot_clock_rate(M, title, save) |
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.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.
| Call | Formula | What each argument means |
|---|---|---|
gr.gravity_at(M, r) |
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) |
Returns the time that passes on the clock at , in s | |
gr.redshift(M, r) |
Returns the redshift of the light when it reaches far away | |
gr.height_shift(height, g) |
Returns the fractional shift in frequency (no unit) | |
gr.gps_drift(altitude) |
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.
| Call | Formula | What each argument means |
|---|---|---|
gr.light_deflection(M, b) |
Returns the angle the ray is bent through, in radians | |
gr.perihelion_shift(M, a, e) |
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)
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.
| Call | Formula | What each argument means |
|---|---|---|
nucleus.semf_binding(Z, A) |
Returns the binding energy , in MeV (, , , MeV; is the pairing term, below) | |
nucleus.semf_terms(Z, A) |
Returns a dict of the five terms ( volume, surface, coulomb, asymmetry, pairing) in MeV, with MeV | |
nucleus.semf_per_nucleon(Z, A) |
Returns the binding energy for each nucleon, in MeV | |
nucleus.stable_z(A) |
Returns the number of protons that makes a nucleus with nucleons most tightly bound | |
nucleus.is_magic(n) |
Returns True if is 2, 8, 20, 28, 50, 82 or 126, the numbers that fill a nuclear shell | |
nucleus.plot_binding(title, save) |
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.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.
| Call | Formula | What each argument means |
|---|---|---|
nucleus.q_value(reactants, products) |
Returns the energy released, in MeV (negative: the reaction needs energy) | |
nucleus.alpha_decay(parent_mass, daughter_mass) |
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) |
Returns the energy released by the beta decay, in MeV | |
nucleus.daughter_number(N0, half_life_parent, half_life_daughter, t) |
Returns the number of daughter nuclei at time | |
nucleus.specific_activity(half_life, molar_mass) |
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.
| Call | Formula | What each argument means |
|---|---|---|
nucleus.closest_approach(K, Z1, Z2) |
Returns the distance it reaches in a head-on approach, in m | |
nucleus.scattering_angle(K, Z1, Z2, b) |
Returns the angle it is turned through (Rutherford scattering) | |
nucleus.coulomb_barrier(Z1, A1, Z2, A2) |
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.
| Call | Formula | What each argument means |
|---|---|---|
radiation.attenuation(I0, mu, x) |
Returns the intensity left, in the same unit as | |
radiation.half_value_layer(mu) |
Returns the thickness that halves the beam, in m | |
radiation.thickness(fraction, mu) |
Returns the thickness needed, in m | |
radiation.mean_free_path(mu) |
Returns the average distance before an interaction, in m | |
radiation.inverse_square(I1, r1, r2) |
Returns the intensity at | |
radiation.plot_attenuation(mus, title, save) |
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.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.
| Call | Formula | What each argument means |
|---|---|---|
radiation.absorbed_dose(energy, mass) |
Returns the absorbed dose , in grays (Gy) | |
radiation.equivalent_dose(dose, radiation_type) |
Returns the equivalent dose , in sieverts (Sv) | |
radiation.effective_dose(*parts) |
Returns the whole-body effective dose , in Sv | |
radiation.gamma_dose_rate(activity, photon_energy, r, mu_en_rho, yield_) |
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.
| Call | Formula | What each argument means |
|---|---|---|
reactor.fissions_per_second(power, energy) |
Returns the number of fissions each second | |
reactor.fuel_use(power, time, molar_mass, energy) |
Returns the mass of fuel fissioned, in kg | |
reactor.energy_per_kg(molar_mass, energy) |
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.
| Call | Formula | What each argument means |
|---|---|---|
reactor.eta(nu, sigma_f, sigma_a) |
Returns the neutrons made for each neutron absorbed by the fuel | |
reactor.four_factor(eta, f, p, epsilon) |
Returns the multiplication factor of an infinitely large reactor | |
reactor.six_factor(eta, f, p, epsilon, p_fnl, p_tnl) |
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) |
Returns the number of neutrons after generations | |
reactor.reactor_period(k, generation_time) |
Returns the time for the neutron population to grow by a factor , in s | |
reactor.moderator_loss(A) |
Returns the average logarithmic energy loss for each collision (1 for hydrogen) | |
reactor.collisions_to_slow(A, E0, E) |
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.
| Call | Formula | What each argument means |
|---|---|---|
fusion.reaction_energy(name) |
Returns the energy released, in MeV | |
fusion.gamow_energy(Z1, Z2, reduced_mass) |
Returns the Gamow energy , in MeV () | |
fusion.tunnel_probability(E, Z1, Z2, reduced_mass) |
Returns the chance of tunnelling through the electric barrier | |
fusion.power_density(n, sigma_v, energy) |
Returns the fusion power in each cubic metre, in W/m³ | |
fusion.triple_product(n, T, tau) |
Returns the Lawson triple product, in keV s/m³ (deuterium–tritium needs about to ignite) | |
fusion.mass_loss_rate(power) |
Returns the mass turned into energy each second, in kg/s | |
fusion.kev_to_kelvin(T) |
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.
| Call | Formula | What each argument means |
|---|---|---|
particle.properties(name) |
Returns a dict: mass (MeV/c²), charge, baryon_number, lepton_number, spin | |
particle.mass(name) |
Returns the mass in MeV/c² | |
particle.hadron(quarks) |
Returns a dict: charge (units of ), baryon_number, kind (baryon, meson, antibaryon), and name (such as 'proton', or None) | |
particle.check(reactants, products) |
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) |
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.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.
| Call | Formula | What each argument means |
|---|---|---|
particle.decay_distance(lifetime, energy, mass) |
Returns the average distance it travels before decaying, in m | |
particle.width(lifetime) |
Returns the energy width , in MeV | |
particle.lifetime(width) |
Returns the lifetime , in s | |
particle.force_range(mass) |
Returns the range of the force, in m | |
particle.cm_energy_fixed_target(E_beam, m_beam, m_target) |
Returns the energy available in the centre-of-mass frame, in MeV | |
particle.cm_energy_collider(E1, E2, m1, m2) |
Returns the energy available when they collide head-on, in MeV | |
particle.bending_radius(p, B, charge) |
Returns the radius of its circle, in m | |
particle.momentum_from_radius(B, r, charge) |
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.
| Call | Formula | What each argument means |
|---|---|---|
particle.barns(x) |
Returns the cross-section in m² | |
particle.mean_free_path(n, sigma) |
Returns the average distance before an interaction, in m | |
particle.interaction_probability(n, sigma, thickness) |
Returns the chance that a particle interacts | |
particle.reaction_rate(luminosity, sigma) |
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)
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
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)
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)
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
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...")
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 / .pairplotNamespace 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")
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() window, in picturesOne line, api.gui(), opens a free, Postman-style tool in your browser. It runs only on your computer and saves your work in a plain folder. Click any picture to enlarge it.
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.