DEV Community: sishiya (しし屋)

Solving the Schrödinger equation in the browser and matching tunneling to theory

sishiya (しし屋) — Wed, 14 Oct 2026 15:28:39 +0000

A particle hits a wall it should not be able to climb, and still gets through with some probability. That is quantum tunneling.

I made a page that shows this by solving the 1D Schrödinger equation numerically, as is. It uses Canvas 2D and no libraries.

https://sishiya.github.io/manabi-lab/quantum/?lang=en#tunnel

When a wave packet hits the wall, part of it bounces back and part leaks through to the other side. At the end, the probability of getting through is compared with the textbook formula.

This post covers the method (Visscher's method) and how I fixed the place where it disagreed with theory.

Equation and units

In units of ħ = 1 and mass m = 1, the equation is:

i ∂ψ/∂t = H ψ,   H = −(1/2) ∂²/∂x² + V(x)

The grid has N = 2000 points, dx = 0.25, time step dt = 0.02, and 80 steps per frame. The wall height is V₀ = 1.

Visscher's method: step the real and imaginary parts in turn

Splitting ψ = R + iI gives two equations:

∂R/∂t =  H I
∂I/∂t = −H R

Visscher's method updates R and I staggered by half a step, in turn: once R is new, the new R is used to advance I. It is about as simple as Euler's method, yet with a small enough time step it is stable and nearly conserves probability.

const c = 0.5 / (dx * dx);
for (let i = 1; i < N - 1; i++) R[i] += dt * (-c * (I[i+1] - 2*I[i] + I[i-1]) + V[i] * I[i]);
for (let i = 1; i < N - 1; i++) I[i] -= dt * (-c * (R[i+1] - 2*R[i] + R[i-1]) + V[i] * R[i]);

Just two loops. The second loop reads the R that the first loop just rewrote.

The initial wave packet is a Gaussian times a wave with momentum k₀ = √(2E), normalized so the total probability is 1:

const g = Math.exp(-((x - X0) ** 2) / (4 * SIG * SIG));
R[i] = g * Math.cos(k0 * x);
I[i] = g * Math.sin(k0 * x);

Absorbing the wave at the edges

If the wave bounces off the grid edges, the transmitted part comes back and spoils the count. So the 160 points at each end form an absorbing layer that shrinks the amplitude a little at a time.

The absorbed probability is summed separately for the left and right edges. What disappears on the left counts as "bounced back", and on the right as "got through", in the final fractions.

for (let j = 0; j < ABS; j++) {
  const f = 1 - 0.006 * ((ABS - j) / ABS) ** 2;
  for (const i of [j, N - 1 - j]) {
    const p = R[i] * R[i] + I[i] * I[i];
    R[i] *= f; I[i] *= f;
    const lost = p * (1 - f * f) * dx;
    if (i < N / 2) absL += lost; else absR += lost;
  }
}

Comparing with theory

For E < V₀, the textbook formula gives the probability of passing a wall of height V₀ and thickness a:

T = 1 / (1 + V₀² sinh²(κa) / (4E(V₀ − E))),   κ = √(2(V₀ − E))

This formula is for a wave with a single energy, though. The simulated packet has a small spread of energies, so it will not match exactly. The page labels it "theory (approximate)".

Where I got stuck: thin walls were way off

At first, grid points inside the wall had V = 1 and the rest 0. Then the thinner the wall, the further it drifted from theory.

The cause was that the wall thickness was being rounded to the grid spacing dx = 0.25. Tunneling depends exponentially on thickness, so even one cell at the edge makes a big difference.

So cells on the wall's edge now get a height proportional to how much of the cell overlaps the wall:

// cells at the wall's edge get partial height, so the thickness isn't rounded to the grid
const ov = Math.max(0, Math.min(x + dx/2, BC + w/2) - Math.max(x - dx/2, BC - w/2));
V[i] = V0 * ov / dx;

After that, it agrees well with theory:

Energy E	Wall thickness	Simulation	Theory
0.6	3	about 1.9%	1.8%
0.6	1.25	about 33.8%	34.0%
0.4	0.75	about 53%	53.3%
1.2 (above the wall)	1.25	about 64%	65.5%

The last row shows that even with energy above the wall, not all of it gets over: about a third bounces back. In classical mechanics it would always get over.

You can check the same values from the browser console:

const E = __qx.EXPS.tunnel; E.T.E = 0.6; E.T.w = 1.25; E.fireNew();
while (E.T.running) E.step();
E.parts();   // [bounced back, inside the wall, got through] → got through ≈ 0.338

Bonus: a Bell test fools you unless you use many pairs

The same page has a Bell test (CHSH) for entanglement. Quantum mechanics predicts S = 2√2 ≈ 2.83, while a "hidden variables" model gives S ≤ 2.

At first I used 2000 pairs each. Then even the hidden-variables model sometimes showed S = 2.06, above 2. That was just statistical scatter.

So I raised it to 20,000 pairs each and also show the standard error of S, √(Σ(1 − E²)/N). Now quantum gives S ≈ 2.83 and hidden variables S ≈ 2.0. When you say "it exceeded 2", you have to show the size of the scatter too.

Main sources

P. B. Visscher (1991) Computers in Physics 5, 596 (the time-stepping method)
D. J. Griffiths, D. F. Schroeter (2018) Introduction to Quantum Mechanics, 3rd ed., Cambridge University Press
R. P. Feynman, R. B. Leighton, M. Sands (1965) The Feynman Lectures on Physics, Vol. III
J. S. Bell (1964) Physics 1, 195
J. F. Clauser, M. A. Horne, A. Shimony, R. A. Holt (1969) Physical Review Letters 23, 880

This page is part of Manabi Lab, a series of small learning apps that run in the browser (English and Japanese). It also has double-slit, polarizer and uncertainty-principle experiments. The code and this post were written together with AI (Claude).

https://sishiya.github.io/manabi-lab/?lang=en

One garden through 14 kinds of eyes: UV, heat, echoes and compound eyes as shader passes

sishiya (しし屋) — Wed, 14 Oct 2026 15:23:27 +0000

I made a page where you walk through the same garden in first person, switching between the senses of a crow, pit viper, bat, fly, honeybee, dog, cat, mantis shrimp and more. Senses humans cannot perceive (ultraviolet, heat, echoes, polarization, near-infrared) are "translated" into colors we can see. It uses Three.js with my own shaders.

https://sishiya.github.io/manabi-lab/animal-senses/?lang=en

This post covers how one garden produces pictures for many different senses.

The views are estimates based on research and cannot be used to diagnose eyes. The app labels its explanations as "established", "estimate" or "effect".

Give objects "per-sense" properties

Besides the color the human eye sees, each vertex in the garden carries these properties (painted by paint()):

Attribute	Contents
`aCol`	color seen by the human eye (linear RGB)
`aP`	UV reflectance, daytime temperature, night temperature, whether it glows
`aA`	opacity to visible light, UV and heat
`aX`	how well echoes return, near-infrared reflectance, whether it is a smooth polarizing surface

For example, sunflower petals have a UV reflectance of 0.03 at the base and 0.58 at the tip. To us they are evenly yellow, but in UV the dark base forms a "target" pattern (values are estimates).

Opacity is split per kind of light because of garbage bags: showing that a bag hides its contents in visible light but not in UV requires different transparency for each kind of light.

Draw the same scene many times with different modes

The scene shader switches what it draws with uMode:

uMode	What it draws
0	visible light
1	ultraviolet
2	heat (temperature / 60)
3	echoes (distance / 40, and how much returns based on surface angle)
4	near-infrared (sunlight + the camera's IR light, which falls off with distance squared)
5	polarization (Fresnel equations for smooth surfaces, Rayleigh scattering for the sky)

Only the needed passes are rendered to render targets, and the compositing shader COMP_FS mixes them with a formula for each animal:

Crow: UV overlaid in purple on visible light
Pit viper: heat blurred over a 320-cell width and overlaid (the pit organ cannot see fine shapes)
Dog, cat: dichromat color formula (Viénot 1999) and blur; the cat gets brighter and softer at night
Mantis shrimp: grayscale, with color quantized to 12 hues only inside a horizontal band that sweeps up and down

You can also split the screen with a human view on the left.

Compound eyes: read a cube map on a hexagonal grid

Flies and bees see almost all the way around. So the page renders a 256 px cube map at the camera position, converts it to azimuth and elevation, and packs 336° across the screen.

That is divided into a hexagonal grid, and each hexagon reads just one direction, at its center. Mipmaps blur it by each ommatidium's acceptance angle. For the bee, a UV cube is rendered too and colors are shifted: "green → red, blue → green, UV → blue".

Bats: echoes come back at the speed of sound

The bat remembers up to six call times. For each pixel, when

age − 2d / 343 ≥ 0   (age = time since the call, d = distance to the object)

the echo has returned, so the pixel lights up and fades over 0.12 seconds. Where the call's wavefront is right now (d ≈ 343·age) is drawn as a blue line.

When a moth is ahead, the call interval shortens from 0.1 to 0.05 to 0.012 seconds (the feeding buzz).

With "hear the echoes" on, the real 80→25 kHz sweep is played stretched 10× as 8→2.5 kHz, with delay and volume set by each object's distance, direction and size.

Nearsight, farsight and presbyopia from the "distance pass"

The distance pass made for echoes also works for the focus errors of human eyes. For each pixel:

needed power = 1/distance + hyperopia (myopia is negative)

Below 0 is the myopic error; anything beyond the accommodation range (Hofstetter's minimum, 15 − 0.25×age) is the hyperopic or presbyopic error. The blur circle's diameter is "pupil diameter × error", with the pupil at 3 mm by day and 6.5 mm at night. You hold a paperback 35 cm away, so with presbyopia you can see just the book go blurry.

Color vision differences use the matrices of Machado et al. (2009).

Careful with flicker

Flies sense fast time, so they should see the flicker of a fluorescent lamp (100 brightness changes per second). On screen, though, the flicker becomes 2–25 times per second, inside the range most likely to trigger photosensitive seizures (3–30).

Measuring the first version, about 9% of the screen changed relative luminance by 0.1 or more. So now the flicker appears only when you turn on "see the lamp flicker", and the brightness swing is reduced to 85–100% (a real lamp drops to about 50%). With prefers-reduced-motion, it cannot be turned on. It is now 0%.

Small gotchas

Drawing a translucent film double-sided without depth writes put the back face over the front, so the bag's contents showed through in UV where they should be hidden → translucent objects are front-face only
Having a GLSL function and a local variable share a name (cat) is risky, so they became catV and dogV
Merging static objects into one mesh (bakeStatic()) cut draw calls from hundreds to dozens

Main sources

F. Viénot, H. Brettel, J. D. Mollon (1999) (dichromat color formula)
G. M. Machado, M. M. Oliveira, L. A. F. Fernandes (2009) (color vision matrices)
T. H. Bullock, F. P. J. Diecke (1956) (sensitivity of the pit viper's pit organ)

The other sources (about 35) are in the app.

This page is part of Manabi Lab, a series of small learning apps that run in the browser (English and Japanese). The code and this post were written together with AI (Claude).

https://sishiya.github.io/manabi-lab/?lang=en

Tracing light around a spinning black hole in the browser (WebGL2, no libraries)

sishiya (しし屋) — Wed, 14 Oct 2026 15:23:14 +0000

The black hole in the movie Interstellar has a flat disk that seems to wrap over the top and under the bottom of the hole. That is not an artistic choice. It is what gravity does to light.

To see it for myself, I made a page that traces light, ray by ray, around a spinning (Kerr) black hole and draws the result. It is plain WebGL2 with no libraries.

https://sishiya.github.io/manabi-lab/black-hole/?lang=en

Click any point on the image and the page shows where that light came from, as a side view and a top view of its path. This post covers how it works and where I got stuck.

What it does

Like ordinary ray tracing, it follows light backwards from the camera. The difference is that light does not travel in straight lines. For every pixel, the page numerically integrates the geodesic equations of Kerr spacetime.

Each ray ends in one of three ways:

It falls into the hole → black (the shadow)
It crosses the disk plane inside the disk → disk color
It escapes far away → look up the star map in its final direction

The equations

Units are G = c = M = 1 (length and time measured in black-hole masses), in Boyer–Lindquist coordinates (r, θ, φ). With the photon energy fixed at p_t = −1, there are two conserved quantities: the angular momentum b = L/E and Carter's constant q = Q/E².

Σ = r² + a²cos²θ,  Δ = r² − 2r + a²
P = r² + a² − ab,  K = (b − a)² + q,  R = P² − ΔK

dr/dλ   = Δ p_r / Σ
dθ/dλ   = p_θ / Σ
dφ/dλ   = (aP/Δ − a + b/sin²θ) / Σ
dp_r/dλ = (R′/Δ − RΔ′/Δ² − Δ′p_r²) / (2Σ)     (R′ = 4rP − Δ′K, Δ′ = 2r − 2)
dp_θ/dλ = (b²cosθ/sin³θ − a²sinθcosθ) / Σ

For light the Hamiltonian is H = 0, which removes the derivative of Σ. I derived these myself and checked them against the appendix of the paper by the team that made the movie's images (James et al. 2015). The sign of φ is flipped compared with the paper, because the paper writes the equations for tracing backwards in time.

In JavaScript (y = [r, θ, φ, p_r, p_θ]):

function deriv(y, b, q, a, out) {
  const r = y[0], th = y[1], pr = y[3], pth = y[4];
  let s = Math.sin(th); const c = Math.cos(th);
  if (Math.abs(s) < 1e-6) s = s < 0 ? -1e-6 : 1e-6;
  const S = r * r + a * a * c * c, D = r * r - 2 * r + a * a;
  const P = r * r + a * a - a * b, K = (b - a) * (b - a) + q, R = P * P - D * K;
  const Dp = 2 * r - 2, Rp = 4 * r * P - Dp * K;
  out[0] = D * pr / S;
  out[1] = pth / S;
  out[2] = (a * P / D - a + b / (s * s)) / S;
  out[3] = (Rp / D - R * Dp / (D * D) - Dp * pr * pr) / (2 * S);
  out[4] = (b * b * c / (s * s * s) - a * a * s * c) / S;
  return out;
}

This is stepped with 4th-order Runge–Kutta, with a negative step because we go backwards.

The camera: an observer swept along with space

Near a spinning black hole, space itself is dragged around in the direction of the spin. So the camera is a ZAMO (zero-angular-momentum observer, Bardeen, Press and Teukolsky 1972), which rotates along with that dragging.

From the pixel direction n, the page computes the photon energy E_F that this observer measures, and the conserved quantities b and q:

const EF = 1 / (z.alpha + z.omega * z.varpi * nph);   // photon energy seen by the camera
const pr = EF * z.rho / Math.sqrt(z.D) * nr;
const pth = EF * z.rho * nth;
const b = EF * z.varpi * nph;                         // angular momentum L/E
const q = pth * pth + c * c * (b * b / (s * s) - a * a); // Carter's constant

alpha is how fast clocks tick there, omega is the angular velocity of frame dragging, and varpi is the circumferential radius.

Step size

Paths bend more sharply close to the hole, so the step is proportional to the distance from the horizon:

const STEP_K = 0.06, STEP_MIN = 0.004;
function stepLen(r, rh) { return Math.max(STEP_K * (r - rh), STEP_MIN); }

A ray stops when:

r < r₊ + 0.01 → the hole (r₊ is the event horizon)
it crosses the equatorial plane between the disk's inner and outer edges → the disk
r > max(2 r_cam, 60) and moving outward → the sky
900 steps pass without any of these → the hole (light circling the photon orbit)

Where I got stuck: a vertical line down the middle

The first render had a flickering vertical line down the center of the image.

Rays in the central column have b ≈ 0, so they pass right over or under the spin axis. There sinθ is small, and the b/sin²θ term in dφ/dλ becomes very large. With a step based only on r, φ jumped a lot in a single step and the integration broke down.

The fix was to shrink the step so that one step never turns more than 0.04 rad:

function rk4(y, h, b, q, a) {
  deriv(y, b, q, a, _k1);
  const lim = 0.04 / (Math.abs(_k1[2]) + Math.abs(_k1[1]) + 1e-6);
  if (Math.abs(h) > lim) h = Math.sign(h) * lim;
  // ...then an ordinary RK4 step
}

The line disappeared.

A related problem: rays with b exactly 0 cross the axis, and θ goes negative (or beyond π). Computing the disk hit point with θ = π/2 then swapped the near and far sides. The shader adds π to φ when sinθ < 0, and the CPU side interpolates between path points to find the position.

Disk color: the frequency shift

The disk gas follows Keplerian orbits. The inner edge is the ISCO (innermost stable circular orbit, 3.83M for a = 0.6), and the outer edge is 18.7M, the same as the figure in the paper.

The frequency shift g of the light is:

g = E_F / (u^t (1 − Ω b))
Ω = 1/(r^{3/2} + a)
u^t = (r^{3/2} + a) / (r^{3/4} √(r^{3/2} − 3r^{1/2} + 2a))

At the movie's viewing angle (a = 0.6), g is about 1.36 on the left side, where the gas comes toward us, and about 0.76 on the band above the hole.

For clarity, the movie left out this shift in color and brightness (the paper explains why). So the page lets you switch between "movie" and "realistic". In realistic mode the color temperature is multiplied by g and the brightness by g⁴: the left side gets bright and the right side gets dim.

Speed: two levels

With up to 900 integration steps per pixel, the full screen cannot be traced every frame.

Coarse level (about 160,000 pixels): retraced in one go whenever the camera moves
Fine level (2× screen resolution, up to 5 million pixels): starts 140 ms after the camera stops, in 192 px tiles from the center outward, as many tiles per frame as the frame time allows

The tracing pass only writes its result (kind, r and φ where it hit, g) to an RGBA32F texture. The disk pattern (which rotates), color and bloom are painted every frame by separate shaders.

Also, the camera's φ (dragging left and right) is a rotation around the spin axis, so no retracing is needed: the page just rotates the sky and the disk pattern. Rays are retraced only when θ, distance or spin change.

On my PC (NVIDIA), the whole fine level takes about 20 ms.

Writing the same equations twice, on GPU and CPU

The ray diagram for a clicked point is retraced on the CPU (double precision, up to 4000 steps). Since the same equations live in both the shader and JavaScript, I keep a check that compares them at the same point: the value read back from the GPU texture against the CPU result.

__bh.renderNow();
__bh.readGeo(0.04, 0.5);       // GPU: [1, ~6.8, ~3.0, ~0.76]  (disk, r, φ, g)
__bh.pickAt(0.04, 0.5).rHit;   // CPU: ~6.8

What the bands in the movie really are

At the movie's angle (a = 0.6, r = 28, θ = 85°), tracing up the central vertical line gives:

Height on screen	What you see	How much the light bent
Upper part	Top face of the disk behind the hole (the band above the hole)	~41°
Just above the shadow edge	Bottom face of the disk in front (light loops under and comes back)	~210°
Just below the shadow edge	Top face of the disk in front	~270°
Lower part	Bottom face of the disk behind the hole (the band below the hole)	~100°

The band above the hole is the top of the disk behind the hole. Its light bends over the top of the hole to reach us. The band below is the underside of the same disk, and its light reaches us by passing under the hole.

Still missing

The thin rings at the shadow edge (second and higher images) are only 1–3 pixels wide and can look dotted
Aberration from camera motion, disk thickness, and the disk lighting itself are not included

Main sources

O. James, E. von Tunzelmann, P. Franklin, K. S. Thorne (2015) "Gravitational lensing by spinning black holes in astrophysics, and in the movie Interstellar", Classical and Quantum Gravity 32, 065001 (arXiv:1502.03808)
J. M. Bardeen, W. H. Press, S. A. Teukolsky (1972) "Rotating black holes: locally nonrotating frames, energy extraction, and scalar synchrotron radiation", The Astrophysical Journal 178
N. I. Shakura, R. A. Sunyaev (1973) "Black holes in binary systems. Observational appearance", Astronomy and Astrophysics 24
K. S. Thorne (2014) The Science of Interstellar, W. W. Norton

This page is part of Manabi Lab, a series of small learning apps that run in the browser (English and Japanese). The code and this post were written together with AI (Claude).

https://sishiya.github.io/manabi-lab/?lang=en

Diving through the digestive tract in 3D: cutting holes where tubes meet, in the shader

sishiya (しし屋) — Wed, 14 Oct 2026 15:13:18 +0000

I made a page where you enter through the mouth or nose, as a piece of food or a capsule endoscope, and travel in 3D through the throat, esophagus, stomach, small and large intestine, or down the windpipe into the lungs. Along the way, stomach acid changes the color and the gates (sphincters) open and close. The view uses Three.js; the tube drawing and the body's calculations are my own.

https://sishiya.github.io/manabi-lab/body-dive/?lang=en

This post covers how the tubes are drawn, the hardest part (where tubes meet), collision, and how time runs inside the body.

A tube = a center line + a radius

Every path in the body is a center line with a radius. Each point is [x, y, z, radius, region], in centimeters. Tube widths are roughly life-size, but the small and large intestines are heavily shortened.

prepTube() resamples the center line evenly and gives each sample a radius, region and a parallel-transport frame (N/B/T). The mesh is built from that, and each vertex also carries the type of fold for its region:

nasal turbinates, tongue, cartilage rings of the windpipe, stomach folds, circular folds of the small intestine, haustra of the colon...

We look from the inside, so tubes are drawn double-sided with normals flipped toward the viewer.

Moving things in the shader

The stomach swelling and shrinking, gates tightening, and peristaltic waves are done in the vertex shader by changing the tube radius. The fragment shader handles holes, noise bumps, lighting and fog.

The hardest part: where tubes meet

In the throat, the path from the nose and the path from the mouth merge, and then split into the windpipe and esophagus. If you simply overlap separate tubes, one tube's wall stands in the way inside the other.

So each tube's shader receives the other tube's shape as a list of capsules (segment + radius), and discards any of its own wall that lies inside the other tube, cutting a hole.

Capsules are only near the junction (with a quick sphere test first)
Folds are removed near the junction; otherwise the hole's edge shifts and you can see the darkness outside

That led to a rule: never put a gate that changes width, or the swelling stomach, at a junction. The hole is static, so if the tube moves, a gap appears.

Because the shape is approximated with capsules, a thin gap can still show at the hole's edge (a known issue).

Collision

You can move within radius × factor − 0.1 cm of each segment (every third sample of the center line), leaving 0.1 cm from the wall.

A bigger margin looks safer, but then you cannot pass through junctions. Once, I could no longer get from the mouth into the throat. After any change, I check that all three entrances (mouth, nose, anus) can be traveled through.

When a gate is closed, its segments get a negative radius and are dropped. On top of that, crossesClosed() cancels any move that crosses a gate. That double guard keeps you from jumping through a thin gate in a single frame.

Time inside the body

Food takes hours to pass the stomach, so the time speed-up differs by region:

Region	Speed-up
Mouth, esophagus	×1
Stomach, small intestine	×60
Large intestine	×240

The body calculation simStep(minutes, region) advances things like:

Stomach: acid changes color and breaks food down. The pylorus opens only when contents become a paste, or during the strong fasting contraction (MMC phase III)
Small intestine: bile turns it yellow. The ileocecal valve opens after 150 minutes in the small intestine
Large intestine: bacterial fermentation, water absorption, browning. It stops in the sigmoid colon and moves to the rectum with a mass movement
Stool form (Bristol scale) is decided from water content

If you choose sake, blood alcohol is computed too (Widmark formula: 60 kg body, 42 L distribution volume, 6.5 g/hour elimination).

The third-person camera

There is also a third-person camera that shows you (the food, capsule, etc.) from behind. At first it sat directly behind you, but near walls it had nowhere to go, moved into the body and flickered. Now the camera is placed behind you along the tube's center line.

Estimates and effects

Shortened lengths, transit times, color changes and so on are "estimates" or "effects". The start screen's "About this app" separates what is established from what is estimated.

For example, the app says the small intestine is "about 3 m in the living body". The often-quoted 6 m is the length after death (OpenStax).

Main sources

J. G. Betts et al. (2022) Anatomy and Physiology 2e, OpenStax
StatPearls (chapters on the esophagus, stomach, duodenum and rectum)
S. J. Lewis, K. W. Heaton (1997) (Bristol stool scale)
MEXT, Standard Tables of Food Composition in Japan (8th edition, 2023 supplement)

The other sources (18 in all) are in the app.

This page is part of Manabi Lab, a series of small learning apps that run in the browser (English and Japanese). The code and this post were written together with AI (Claude).

https://sishiya.github.io/manabi-lab/?lang=en

Recomputing 4.6 billion years of Earth with a box model that never reads the timeline

sishiya (しし屋) — Wed, 14 Oct 2026 15:13:13 +0000

I made a page that watches 4.6 billion years of Earth, from its birth to the appearance of humans, from one fixed spot. Move the time slider and the land, sky, sea and living things change to match that age. The 3D view uses Three.js.

https://sishiya.github.io/manabi-lab/evolution/?lang=en

The scenery is drawn from what research has found (the record). But the record alone does not show why things happened.

So I built a box model that never reads the timeline. It computes temperature, CO2, oxygen and sea animals only from outside drivers such as the Sun, volcanoes and weathering, and you can compare it with the record. This post covers how the model is built, and where it does and does not match.

Inputs are only "outside drivers"

runModel() reads nothing but the outside drivers in FORCE:

Brightness of the Sun (dimmer in the past)
Volcanoes (and giant eruptions)
How easily rock weathers (amount of continent, supercontinent breakup, mountain building)
Weathering by plants, methane-making microbes, burial of organic carbon, how reducing volcanic gas is
Earth's orbit, humans

It does not read the timeline (when snowball Earths or mass extinctions happened). The point is to see whether those come out of the calculation.

The grid has 4000 points (past and future), stepped in 0.02-million-year slices. One run takes about 0.22 seconds.

The CO2 thermostat

CO2 changes with the difference between what volcanoes add and what rock weathering removes:

dC = 700 × (volcanoes − weathering)
weathering = weath × bio × (plants) × (C/280)^0.3 × exp((background temp − 14)/13) × √ocean

Weathering is faster when it is warm, so warming lowers CO2 and cools the planet, and cooling raises CO2 and warms it. That is Earth's long-term thermostat. During a snowball Earth, weathering drops to 3%.

Temperature and snowball Earth

temp = 14 + Sun 140×(S−1) + CO2 effect + methane 18×exp(−O2/0.003) + orbit + water vapor + ice −40×(I − 0.03)

With x = ln(C/280), the CO2 effect is 4.33x + 0.25x² for x ≥ 0 and −12×(1 − exp(4.33x/12)) for x < 0.

Ice grows partially as it cools. Below 5 °C the whole planet freezes, and it thaws when the ice-free temperature would exceed 30 °C.

This includes both "ice reflects sunlight → even colder" and "weathering stops while frozen → volcanic CO2 builds up → eventually it melts".

Oxygen and the ozone layer

Oxygen is set by organic carbon burial B minus reducing volcanic gases R. R was larger in the past (0.75×(A/2.4e9)^1.2 + 0.15), and oxygen starts to rise once B − R exceeds 0.002. More oxygen removes methane's greenhouse effect, which changes the temperature too.

Results compared with the record

What emerged on its own:

Event	Model	The app's record
Huronian glaciation	snowball 2.396–2.239 billion years ago	about 2.4–2.1 billion years ago
Cryogenian snowballs (two)	734–687 and 683–649 million years ago	two, about 720–635 million years ago
End-Permian (P-T) extinction	sea animal diversity −50%	the largest mass extinction
End-Cretaceous (K-Pg) extinction	−36%	dinosaurs and others die out

The Late Devonian and end-Triassic extinctions also appear in the model. Today comes out at about 15.6 °C and about 440 ppm CO2.

What did not match:

The end-Ordovician extinction comes out at only about −10% (the record is about −50%)
Average mismatch over all time: temperature ±6 °C, CO2 ×1.7, oxygen ±1.4%, sea animals ±27%

The "how it works" panel overlays the record (dotted) and the model (solid) in four graphs and shows the mismatch as it is.

The future and "what if"

Since the model never reads the record, it can keep running into the future.

Human CO2 is removed with the Bern carbon-cycle response (the ocean takes thousands of years, rock weathering tens to hundreds of thousands)
In the middle scenario (SSP2-4.5), 2100 is +2.45 °C (about the same as IPCC AR6's +2.7 °C)
In about 600 million years CO2 drops below 150 ppm and plants can no longer grow (matches the literature)
The oceans are lost in 2–3 billion years (the literature says 1–2 billion; the model is late)

You can also make "what if" worlds. For example, with "no oxygen-making life", methane's greenhouse effect remains and none of the three snowball Earths happen.

A gotcha: bulk replace in PowerShell

I used PowerShell to rewrite code in bulk and failed three times:

With only one replacement pair, @( @("a","b") ) gets flattened and replaces one character at a time (write @(,@("a","b")))
Adding 'a'+"`n"+'b' inside an array: the comma binds first and breaks the pair
PowerShell expands ${...} inside double quotes

Now the replace function refuses to run unless each pair is a two-element array.

Main sources

The sources for the record and the model's shape are listed in the app under "About this app". The future comparison uses IPCC AR6 (2021).

This page is part of Manabi Lab, a series of small learning apps that run in the browser (English and Japanese). The code and this post were written together with AI (Claude).

https://sishiya.github.io/manabi-lab/?lang=en