DEV Community: sishiya (しし屋)

Solving WiFi waves on 0.75 cm cells in WebGL2 (FDTD, and correcting for a coarse grid)

sishiya (しし屋) — Thu, 22 Oct 2026 01:31:16 +0000

How much do walls weaken WiFi, and where in a home is the signal hard to reach? I made a page that solves the electromagnetic wave equation directly on the GPU and shows the waves spreading and the received strength. Plain WebGL2, no libraries.

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

You can move the router around a two-bedroom apartment and draw extra walls. This post covers how it is solved on the GPU, and how I fixed walls weakening the signal far too much because the grid was too coarse.

The equation

A 2D wave equation for only the electric field Ez perpendicular to the floor (TM mode; in 2D, Ez alone describes it exactly):

∂²E/∂t² = (1/με) ∇²E − (σ/ε) ∂E/∂t

cell size DX = 0.75 cm, 1280×800 cells (9.6 m × 6 m)
Courant number 0.5 → one step is 12.5 ps

The wavelength at 2.4 GHz is about 12.5 cm, so one wavelength is about 17 cells.

GPU: ping-pong between two textures

State lives in two RGBA32F textures, swapping read and write each step:

Channel	Contents
r	current Ez (u(n))
g	previous Ez (u(n−1))
b	exponential moving average of u² (for signal strength)

Rendering to floating-point textures needs EXT_color_buffer_float.

Materials are stored as IDs 0–7 in an R8 texture. The shader looks up coefficients from uniform arrays by ID:

uK2[id] = KC²/εr (metal is −1, forcing Ez = 0)
uA[id] = σ·DT/2ε (loss)

The outer 40 cells are a sponge with damping growing as the cube, to remove reflections.

Converting 2D falloff to 3D

In 2D, wave strength falls off as 1/r with distance r. In real 3D free space it is 1/r².

So at startup, the average of "strength × r" in an empty field is measured as C, and received strength is:

received [dBm] = 20 dBm − FSPL_3D(d) + 10·log10(mean u² · r / C)

Distance loss uses the 3D formula (free-space path loss), and only walls, reflections and interference come from the simulation.

Stuck: glass at −14 dB at 5 GHz

Material values (relative permittivity εr and conductivity σ) come from ITU-R P.2040. Used as they are, walls weakened 5 GHz far too much:

Material (5.2 GHz)	Simulation without correction	Analytic value
Glass	14.4 dB	about 3 dB
Concrete (18 cm)	59.1 dB	23.0 dB

The cause is grid coarseness. On a 0.75 cm grid, one wavelength inside glass at 5 GHz is only about 3 cells. Wave propagation (numerical dispersion) then drifts from reality.

But a finer grid costs cells × steps (in 2D, halving the cell size costs 8×).

The fix: bisect for "values giving the same transmission loss"

So, for each material and frequency, find the εr and σ that give the analytic transmission loss on this grid, and use those:

a 1D strip test that hits a slab with a plane wave, stripDb(f, id, er, sig)
the analytic transmission through the slab with complex numbers, slabDb() (normal incidence)
bisection for the εr and σ where stripDb equals slabDb
store them in the material table as eff: [[er, σ at 2.4GHz], [er, σ at 5GHz]]

Checking, all 5 materials × 2 frequencies agree within 0.01 dB (e.g. concrete 13.37 dB at 2.437 GHz and 23.01 dB at 5.2 GHz).

This is not "bending the physical values". It cancels the coarse grid's error using a measurable quantity, transmission loss. The page says so as a "grid correction". Adding materials or changing thicknesses means regenerating them.

The gap experiment

One scene has an 8 cm-thick metal plate with gaps of 20, 8, 4 and 2 cm. When a gap is narrower than half a wavelength (6.2 cm at 2.4 GHz, 2.9 cm at 5 GHz), waves can hardly pass.

In the simulation, the 4 cm gap gives about −30 dB at 2.4 GHz and about −14 dB at 5 GHz, and nothing gets through the 2 cm gap at either frequency.

Checking for flicker

The wave display has fine moving stripes, so a single point changes brightness 8–29 times per second. Against the photosensitive-seizure guideline (WCAG), I measured brightness changes averaged over 32 px squares. Areas changing a lot more than 3 times per second covered at most 5.9% of the screen, under the guideline (about 10% of the screen).

Main sources

ITU-R P.2040-3 (2023) (electrical properties of building materials)
ITU-R P.525 (free-space attenuation)
IEEE 802.11 (WiFi frequencies)
A. Taflove, S. C. Hagness (2005) Computational Electrodynamics: The Finite-Difference Time-Domain Method, 3rd ed.
C. Gabriel (1996) (electrical properties of the human body)

Signal thresholds, router output and body values are rough values. Nothing here indicates effects on the body.

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

Tearing the Earth with your finger: 16,000-particle SPH in a Web Worker

sishiya (しし屋) — Thu, 22 Oct 2026 01:25:55 +0000

"God's-Eye Map" shows the Earth with CesiumJS and lets you travel from the ground to space using public data that needs no API keys, such as Japan's GSI maps and the PLATEAU 3D city models.

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

Inside it there is a button that appears only while you look at Earth from space: "✋ Tear off a piece of Earth". Pinch part of the Earth and pull, and a hole opens, hot rock comes out and seawater boils.

This post covers how that calculation (SPH, a particle method for fluids) is built and where I got stuck.

At Earth's scale, gravity is more than 100 times the strength of rock, so rock can be treated as a fluid. Still, small fragments would really keep their shape better, and other things are simplified (the app labels them "estimate").

Building an Earth of particles

The Earth is about 16,000 particles. Spacing s is about 408 km, the smoothing length h = 1.1 s, and each particle has about 45 neighbors within 2h.

Layer boundaries (inner core 1221 km, core–mantle boundary 3480 km) follow PREM, the standard Earth model.

Density uses an equation of state for each layer (Murnaghan), solved so that gravity and pressure balance:

Layer	ρ0 (kg/m³)	K0 (GPa)	K′
Mantle	3500	160	4
Outer core	6100	90	4
Inner core	7800	170	4.5

The result: central pressure 340 GPa, mass 5.66×10²⁴ kg (5% lighter than the real Earth). Densities are close to PREM but not the same.

Particles on a lattice get a mass of "density × s³", and the planet is relaxed for 300 steps (multiplying speed by 0.9 each step to stop the shaking) before starting.

Inside the SPH

density is tracked with the continuity equation; summing particle masses gave too little density at the surface and wrong pressure
pressure is never below 0 (no tension), so things tear easily
artificial viscosity (α = 1, β = 2), with its work turned into heat
one step is about 10 seconds, set by the speed of sound (about 11 km/s) and particle spacing

Summing gravity over every pair is heavy, so it is computed as spherically symmetric. For each body (Earth, fragment), the space is cut into 256 shells from its center, and each particle is pulled only by the mass inside its radius.

Temperature is tracked too. Released pressure cools by adiabatic expansion, and the energy of material collapsing inward becomes heat. Compared with the melting onset temperature (solidus), it also gives a "partly molten fraction".

Running it in a Web Worker

The calculation runs in one Web Worker. One step takes about 30 ms (on a 16-core PC). Drawing happens on the main thread, with particles drawn as round point sprites in a shader.

window does not exist inside a Worker, so shared code lives on self.GV.

Stuck: how fast the "hand" moves

The hardest part was not the calculation but how to move the pinching hand.

Particles inside a sphere of radius Rg around the pinched point are pulled toward the hand with a spring (critically damped).

1. The hand was too fast
Using the finger's motion directly as the hand position stirred the Earth with the fragment, reaching 44 km/s and +350,000 °C. Now the hand approaches the pointed spot with a speed limit.

2. The hand was too slow
With a limit of 3 km/s, the hole filled in before it tore off, and you could not see a deep hole. Now it is 8 km/s while pulling and 15 km/s while carrying.

3. Carrying in a straight line scraped the ground
Carrying the fragment in a straight line skimmed the surface and blasted away the ground (+100,000 °C). Now it first moves out to 14,000 km from the center, then swings around.

While you pull (until it tears), the calculation runs more than 300× faster.

Stuck: deciding what "cracked"

I wanted the cracked surface around the hole to glow as hot rock.

At first, "cracked" was decided by how far a particle moved from its starting place. Then the whole Earth wobbling made a wide area count as cracked.

So it now uses spacing to neighboring particles. If the spacing to surface particles that were neighbors at the start (up to 6) changes by 35% or more, it is "cracked". When the whole Earth wobbles, neighbor spacing barely changes.

Oceans and steam

The ocean averages 3.7 km, much thinner than a particle, so it is a separate set of light particles (about 3600). Ocean locations come from checking an Earth photo (NASA's Blue Marble) at 7000 points on the sphere and keeping only blue spots.

At first, water particles slipped into gaps between surface particles, touched the inner rock below, and half the ocean was "boiling". The fix was to treat only rock that is truly exposed at the surface as hot.

Main sources

A. M. Dziewonski, D. L. Anderson (1981) Preliminary reference Earth model, Physics of the Earth and Planetary Interiors 25
J. J. Monaghan (1992) Smoothed particle hydrodynamics, Annual Review of Astronomy and Astrophysics 30
T. Katsura et al. (2010) Physics of the Earth and Planetary Interiors 183 (adiabatic temperature of the mantle)
T. Lay, J. Hernlund, B. A. Buffett (2008) Nature Geoscience 1 (heat at the core–mantle boundary)
M. M. Hirschmann (2000) Geochemistry, Geophysics, Geosystems 1 / G. Fiquet et al. (2010) Science 329 (mantle melting onset)

Sources and terms of use for maps, photos and 3D city models are shown 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

Teaching spring-and-joint creatures to walk with a genetic algorithm (hand-written 2D physics)

sishiya (しし屋) — Thu, 22 Oct 2026 01:20:33 +0000

I made a page where creatures made only of joints (balls) and muscles (springs) evolve with a genetic algorithm (GA) to travel as far as possible. Besides flat ground, you can pick a slope, ice, the Moon or underwater. The physics and the GA are both hand-written, no libraries.

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

Creatures start with random shapes and rhythms, and over the generations they find bodies and ways of moving that go forward in their environment. This post covers how the bodies and physics are built, and how I kept the GA from exploiting gaps in the physics.

Genes = the body's blueprint

One creature's genes hold:

positions of 3–8 joints (inside a 1 m square) and the foot friction coefficient (0.2–1.2)
for each muscle: the two joints it connects, contracted and extended lengths (0.15–1.6 m), spring constant (150–600 N/m), when in the cycle it starts contracting and for how long
the body's rhythm period (0.4–2.0 seconds)

A muscle is a damped spring (damping 0.5 of critical) that heads toward the short length during its contraction window and the long length otherwise.

Physics: 15 seconds at 1/120-second steps

joints are 2 kg balls, moved for 15 seconds at 1/120-second steps (1800 steps)
the ground pushes back penetrating joints and removes downward speed, with no bounce
sideways speed is reduced by friction and a "grip on the ground" force

The score is the center of mass's x after 15 seconds minus its starting x.

Environments just change gravity or friction:

Environment	What changes
Slope	gravity tilted 10° (the picture is rotated)
Ice	friction × 0.1
Moon	gravity 1.62 m/s²
Underwater	gravity 0 (neutral buoyancy), quadratic drag on speed perpendicular to each muscle rod

The slope started at 15°, but even after 100 generations creatures only got 4 m, so it became 10°.

Preventing bodies from passing through themselves

The first version had no collision between body parts. That allowed motions impossible in a plane, like a joint passing through its own muscle rod. Where the physics has a gap, a GA will find solutions that use it.

So from v3, there is collision, collideBody():

joints are 5.5 cm-radius balls, muscles 6 cm-thick rods
penetration is pushed back in proportion to mass, and approaching speed is removed (no bounce)
each joint remembers which side of a muscle it is on, and is put back if it ends up on the other side (no tunneling even when fast)

After adding collision, I checked that momentum is conserved and that mechanical energy does not increase for a body with its muscles stopped.

Speed is about 100 ms per generation of 100 creatures (15–35 ms before collision).

A body whose physics broke scores 0

A body with strong springs might make the integration diverge and fly off. Since the score is distance traveled, that would count as "went far".

So if speed exceeds 60 m/s, the calculation is treated as broken and the score is 0.

Only bodies buildable in a plane are born

A body whose joints are too close or whose muscles cross cannot exist in a plane. So bodies are born only if they satisfy validBody():

joints at least 13 cm apart
joints at least 10.5 cm from muscles not attached to them
no muscles crossing

If a mutation breaks the shape, it retries up to 10 times; otherwise the parent is kept. If crossover's mixing of joint positions breaks it, positions come from one parent.

Inside the GA

Two selection methods:

Top half: the top half survives as is, and each makes one child
Tournament: only the top two survive; for the rest, the stronger of two random creatures becomes a parent

Mutations nudge values on a log-normal scale, and also change the body: adding or removing joints, adding or removing muscles.

Crossover starts from one parent's shape and takes same-numbered joints, muscles linking the same two points, and the period from the other parent with 50% chance. Between parents of different shapes, only the overlapping parts mix (a simplification).

A small gotcha

All scripts share names at the top level, so naming the body thickness BODY collided with the font setting BODY in the drawing code and stopped loading. Now I grep for any new uppercase name before using it.

Main sources

J. H. Holland (1975) Adaptation in Natural and Artificial Systems
K. Sims (1994) "Evolving Virtual Creatures", SIGGRAPH 94

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

Computing how 'old' the air in each room is (ventilation network + 2D flow + age of air)

sishiya (しし屋) — Thu, 22 Oct 2026 01:15:12 +0000

How much faster does air get replaced with two windows open instead of one? With only the exhaust fan running, where does air stagnate?

I made a page that computes this for a two-bedroom apartment and colors how old the air is in each spot (the age of air). No libraries.

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

The calculation has three stages:

How much air the whole home exchanges (ventilation network)
How air moves through the home (2D flow)
How long ago the air at each spot came in from outside (age of air)

This post covers the equations for each and where I got stuck in stage 3.

1. How much air: find the indoor pressure by bisection

Windows, the front door, vents and leaks use the orifice equation:

Q = αA · sign(Δp) · √(2|Δp| / ρ)      Δp = wind pressure on the outer wall − indoor pressure p

Exhaust fans use a straight-line approximation of "less flow the harder they pull" (Q = q·(1 − Δp/p_max)). Wind pressure comes from a wind-pressure-coefficient formula for low-rise buildings (Swami and Chandra), using the angle between the wall and the wind.

Once indoor pressure p is known, every opening's flow is known. So bisection finds the p where inflow equals outflow. Closed doors still connect through the gap underneath, so normally the whole home is treated as one space.

2. Flow in the home: a 2D MAC grid

A horizontal plane about 1 m above the floor is split into 10 cm cells (96×60), on a MAC grid with a 0.05-second step.

Order of one step:

Forcing from jets (circulator fan, air conditioner, wind coming in through open windows)
Advection (semi-Lagrangian)
Viscosity and drag
Projection (pressure equation, 40 SOR iterations, starting from the previous solution)

Flows from stage 1 enter as sources in cells just inside openings and sinks in fan cells.

One trick: 2D flow = "real flow ÷ ceiling height 2.4 m". That keeps the volume equal to the real one, so the exchange time V/Q comes out right.

Closed doors stay fluid cells, but the face permeability is set to 0.05. What is pushed by pressure difference (through the gap under the door) gets through; stirring flow mostly does not.

3. Age of air: "air ages one second per second"

The age of air is how long it has been since that air came in from outside. Thinking of it like this, it can be solved with an ordinary advection–diffusion equation:

air coming in from outside is 0 seconds old
in every cell, air ages one second per second (source S = 1)

Each cell gets a finite-volume, upwind, implicit (backward Euler) equation:

(DX/Δt)(a − a_old) + Σout u·a − Σin u·a_n + q_out·a + (D/DX) Σ(a − a_n) = DX · S

Smoke concentration uses the same equation with S = 0.

Stuck 1: time would not advance

At first I solved this system with a few Gauss–Seidel sweeps. With a large Δt per frame, the age of air barely moved.

Adding diffusion makes the problem elliptic, and a few sweeps do not converge. Switching to BiCGSTAB (with diagonal preconditioning) made it advance properly.

Stuck 2: subtracting twice at inflow cells

In cells where air enters through a window, I had added a term "as much air leaves this cell as enters". But entering air leaves through the faces to the neighbors, so that subtracted it twice. Now the outflow term is only for exit cells and for exchange at a single open window.

The check: the age at the exit equals V/Q

Both mistakes were found with one check.

The mean age of the air leaving (weighted by flow) equals the nominal time constant τ = V/Q, regardless of the flow pattern (Sandberg). However much stagnation there is inside, this holds at the exit.

At first it came out at 1/10 of τ, which is how I noticed. After the fixes, every scene gives 0.999–1.000 of τ.

The same ratio gives the air exchange efficiency:

air exchange efficiency = τ / (2 × mean age of air in the home)

Speed tricks

Flow is computed at most 3 × 0.05 seconds per frame (about 9× speed)
At faster settings (1 hour per second, etc.), the flow pattern is frozen and only the age of air and smoke advance
When one frame's Δt exceeds τ/6, it is split in two

A frame's computation stays under 15 ms.

Main sources

M. Sandberg (1981) "What is ventilation efficiency?", Building and Environment 16
D. Etheridge, M. Sandberg (1996) Building Ventilation: Theory and Measurement, Wiley
M. V. Swami, S. Chandra (1988) ASHRAE Transactions 94 (wind pressure coefficients for low-rise buildings)
P. R. Warren, L. M. Parkins (1985) "Single-sided ventilation through open windows"

Fan flows, leak sizes and so on are typical values for products and homes.

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

Solving an influenza-and-immunity model in the browser (ODEs, RK4, and three gotchas)

sishiya (しし屋) — Thu, 22 Oct 2026 01:09:51 +0000

When a virus or bacterium gets into the body, how does the immune system fight it? I made a page that computes the course of an infection with a model made of ordinary differential equations and shows it as graphs and a microscopic cross-section. Canvas 2D, no libraries.

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

You can pick influenza, COVID-19, norovirus, Staphylococcus aureus, pneumococcus or E. coli, and try medicines and vaccines. This post covers how the influenza model is built and where the numerics tripped me up.

The numbers in this app are rough values, tuned to literature model shapes so that well-known courses (days of fever and so on) come out. They cannot be used for diagnosis or treatment decisions.

Base: a target-cell-limited model

Influenza builds on the "target-cell-limited model with an eclipse phase" of Baccam et al. (2006). Cells of the nose and throat lining are split into five fractions (summing to 1):

Symbol	State
U	not yet infected
R	protected by interferon
E	infected but not yet releasing virus (eclipse)
I	releasing virus
D	dead or damaged

The total is N0 = 4×10⁸ cells (Baccam 2006, upper airway). Virus V is kept as a count.

Main values:

eclipse kE = 4/day (about 6 hours)
infected-cell lifespan del = 1.5/day
virus released per cell Y = 2000 (particles; far fewer are infectious)

On top of that I added interferon, NK cells, macrophages, killer T cells, B cells and antibodies.

A small dose does not infect: saturating clearance

In the original model, a single virus will always grow under some conditions. In reality, a small amount of virus does not cause infection.

So clearance by macrophages and mucus was made saturating:

clearance rate = Cm · M · V / (1 + V/Vs)     (Cm = 150, Vs = 2000)

While V is small, each virus is cleared quickly and clearance beats growth. When V gets large, clearance levels off and cannot keep up.

As a result, in a healthy adult illness starts above roughly 100,000 viruses (an estimate). Human nasal challenge studies use around 10⁵–10⁷ TCID50, so I think the order of magnitude is not off.

For a healthy adult with 100,000 viruses:

fever on day 2, up to 38.8 °C
peak on days 3–4, fever gone on day 6
virus gone on day 8

Gotcha 1: adding antibodies blew up to 10¹³

I was using ordinary 4th-order Runge–Kutta (RK4) with a fixed step. When antibody neutralization was fast (around 1000 per day), the integration diverged and counts jumped to 10¹³.

When something decays quickly, a too-large step overshoots (a stiff equation). So inside each step, the solver now subdivides automatically so that the fastest decay rate × step < 0.8:

let left = h, n = 0;
while (left > 1e-12 && n++ < 400) {
  path.f(y, P, k1);
  let rate = 0;
  for (const v of vs) if (k1[v] < 0 && y[v] > 1e-30) rate = Math.max(rate, -k1[v] / y[v]);
  const s = Math.min(left, 0.8 / Math.max(rate, 1e-9));
  // ...an ordinary RK4 step of size s (intermediate values clamped at 0)
  left -= s;
}

The decay rate -k1[v] / y[v] is "what fraction is lost per day", so its inverse is a rough timescale for that quantity.

Gotcha 2: 10⁻⁵⁰ of a cell growing back

While testing norovirus, the virus started growing again months after it should have been gone.

It turned out infected cells were left at absurdly small amounts like 10⁻⁵⁰. In the equations, even that can slowly grow.

So when almost nothing is left, it is set to 0:

// less than one particle / cell left = none (otherwise 10⁻⁵⁰ of a cell can grow back in the equations)
if (path.kind === 'virus') { if (y.V < 1e-3 && (y.I + y.E) * path.N0 < 1e-3) { y.V = 0; y.I = 0; y.E = 0; } }

A model that uses continuous quantities for counts needs "counts are whole numbers" put in somewhere.

Gotcha 3: inflammation without bacteria

In the S. aureus model, inflammation never ended even after the bacteria were gone.

The inflammation signal also came from "tissue damage", and damage also came from neutrophils' collateral harm. Signal → neutrophils → damage → signal made a loop. The fix: damage only produces a signal while bacteria are present.

Medicine: make "time" part of the state

How well a drug works depends on time since the first dose. So time t itself is one of the state variables (dt/dt = 1). Then any equation can compute the drug effect from y.t.

Oseltamivir blocks 85% of virus release from cells (5 days, half-life about 8 hours). Starting about 5 hours after fever begins shortens the fever by about a day, roughly the shortening reported in clinical trials.

How I check it

Every time, all six pathogens × four bodies × every memory option run to the end, checking for errors:

for (const pk of PATHOGEN_ORDER) for (const body of BODY_ORDER)
  for (const mem of Object.keys(PATHOGENS[pk].memOptions)) {
    const B = BODIES[body];
    __ib.set({ pk, body, mem, inn: B.inn, adp: B.adp, dose: PATHOGENS[pk].doseDefault, speed: 2, playing: true });
    __ib.frame(1200);
  }
__ib.err

Main sources

P. Baccam, C. Beauchemin, C. A. Macken, F. G. Hayden, A. S. Perelson (2006) "Kinetics of influenza A virus infection in humans", Journal of Virology 80, 7590
P. F. M. Teunis et al. (2008) Journal of Medical Virology (how many noroviruses can infect)
S. D. Elek (1957) Annals of the New York Academy of Sciences (how many S. aureus cause skin infection)
Immunology textbooks (such as Janeway's Immunobiology)

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

Inside a pupa you can rewind: deriving every frame from t alone (Canvas 2D)

sishiya (しし屋) — Fri, 16 Oct 2026 12:26:59 +0000

A caterpillar becomes a pupa and then an adult. Inside the pupa, much of the larval body breaks down and the adult body is rebuilt.

I made a page that shows that inside as a side-view diagram, which you can scrub back and forth with a single time slider. The model insect is the tobacco hornworm moth (Manduca sexta), the best-studied species. Canvas 2D, no libraries.

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

This post is about how I made sure that "dragging the slider back returns the same picture".

The rule: everything on screen comes from t

Time t is days since pupation. Wandering (when the larva walks around looking for a place to burrow) is −4, and emergence is 18 (18–20 days in the lab).

Ordinary animation advances the previous frame's state a little at a time. Then you cannot build the picture for a time you scrub back to.

So everything on screen is determined only by t and the real-time clock clock. No previous-frame state is kept. Even the fragments scattered from breaking organs are computed from t.

One timetable `TL` holds everything

What happens when lives in one table, TL. Each row is tagged with how solid its source is:

sure: known with dates for Manduca sexta
est: an estimate

The "estimate" labels in the on-screen text follow these tags. To change timing, you only edit this table.

Hormone levels are functions of t too. Ecdysone (the molting hormone), for example, has a peak in the prepupa; the adult peak rises from P2.5, is highest at P7–9 and drops sharply at P10. Only the shape (0–1) is shown, not absolute amounts.

Every stage uses "the same 15 boundaries"

Larva, prepupa, pupa and adult look completely different. Drawing each stage separately would make it impossible to connect them.

So every stage's body is described by the same 15 segment boundaries:

s = 0 is the tip of the head, 1 the head/prothorax boundary, ..., 14 the tip of the tail
each boundary has position X, dorsal height D, ventral V and backbone height Y

Organs are placed at (s, v), with v = +1 on the back surface and −1 on the belly. Morphing from larva to pupa is then just blending values boundary by boundary. Organs that have both a larval and an adult position blend by each organ's own progress.

Breaking organs: fix the randomness by place

Larval muscles and salivary glands break down early in the pupa. breakingTube() draws this:

split the tube into points
give each point a slightly shifted "progress" using randomness
go "connected → fragments → gone" by that progress

// a = how much of the organ remains (1 → 0). f = how broken each point is
const f = R.map((p, i) => clamp((1 - a) * 1.7 - hash(seed, i) * 0.7, 0, 1));

The randomness is a fixed value (a hash) made from the point's index i. Drawing new random numbers every frame would move the fragments each time you scrub back. With a hash, the same t always gives the same fragments in the same places. Blood cells gather to clean them up.

Where I got stuck: wings poking out of the back

Wings start as small pouches inside the thorax and turn inside out in the prepupa.

At first I placed the wings pointing "from the base toward the belly", and the wing body poked out of the back. Putting the base high on the side of the thorax, with the leading edge sloping down toward the belly, matched a real pupa.

Also, rotating the wings straight to the adult orientation right after emergence made them look, from outside, like "wings appearing out of nowhere". So now they expand in the same place and orientation as in the pupa (15–70 minutes after emergence), then rotate at the base and fold (70–85 minutes). These times are estimates: I could not find a study that gives when Manduca starts and finishes expanding its wings, and the app labels them "estimated" on screen.

Sharing the slider's time

Many things happen within hours around emergence. In plain days, that would be a sliver of the slider. So the slider-to-time mapping changes by section:

Share of slider	Time
76%	−6 to 18 days (until emergence)
8%	emergence to 15 min (out of the case, out of the soil)
7%	to 100 min (wings expand)
3%	to 6 hours
6%	to 1.5 days

How I check it

Every time, I run through all time in the three views (see-through, cut open, outside) and check for errors:

for (const m of ['in', 'cut', 'out']) {
  __mm.VIEW.mode = m;
  for (let i = 0; i <= 300; i++) {
    __mm.setT(uToT(i / 300));
    if (__mm.err) throw __mm.err;
  }
}

Because the picture depends only on t, calling setT at scattered times gives the same pictures as playing in order.

Main sources

The app's "About this app" lists 15. A few of them:

J. W. Truman, L. M. Riddiford (1999) Nature (hormones of complete metamorphosis)
M. Jindra, S. R. Palli, L. M. Riddiford (2013) Annual Review of Entomology (Kr-h1, broad, E93)
L. M. Schwartz, J. W. Truman (1982) Science (muscles between abdominal segments die after emergence)
K. A. Hegstrom, L. M. Riddiford, J. W. Truman (1998) Journal of Neuroscience (remodeling of larval muscles and the flight muscles)

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

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

DEV Community: sishiya (しし屋)

Solving WiFi waves on 0.75 cm cells in WebGL2 (FDTD, and correcting for a coarse grid)

The equation

GPU: ping-pong between two textures

Converting 2D falloff to 3D

Stuck: glass at −14 dB at 5 GHz

The fix: bisect for "values giving the same transmission loss"

The gap experiment

Checking for flicker

Main sources

Tearing the Earth with your finger: 16,000-particle SPH in a Web Worker

Building an Earth of particles

Inside the SPH

Running it in a Web Worker

Stuck: how fast the "hand" moves

Stuck: deciding what "cracked"

Oceans and steam

Main sources

Teaching spring-and-joint creatures to walk with a genetic algorithm (hand-written 2D physics)

Genes = the body's blueprint

Physics: 15 seconds at 1/120-second steps

Preventing bodies from passing through themselves

A body whose physics broke scores 0

Only bodies buildable in a plane are born

Inside the GA

A small gotcha

Main sources

Computing how 'old' the air in each room is (ventilation network + 2D flow + age of air)

1. How much air: find the indoor pressure by bisection

2. Flow in the home: a 2D MAC grid

3. Age of air: "air ages one second per second"

Stuck 1: time would not advance

Stuck 2: subtracting twice at inflow cells

The check: the age at the exit equals V/Q

Speed tricks

Main sources

Solving an influenza-and-immunity model in the browser (ODEs, RK4, and three gotchas)

Base: a target-cell-limited model

A small dose does not infect: saturating clearance

Gotcha 1: adding antibodies blew up to 10¹³

Gotcha 2: 10⁻⁵⁰ of a cell growing back

Gotcha 3: inflammation without bacteria

Medicine: make "time" part of the state

How I check it

Main sources

Inside a pupa you can rewind: deriving every frame from t alone (Canvas 2D)

The rule: everything on screen comes from t

One timetable TL holds everything

Every stage uses "the same 15 boundaries"

Breaking organs: fix the randomness by place

Where I got stuck: wings poking out of the back

Sharing the slider's time

How I check it

Main sources

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

Equation and units

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

Absorbing the wave at the edges

Comparing with theory

Where I got stuck: thin walls were way off

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

Main sources

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

Give objects "per-sense" properties

Draw the same scene many times with different modes

Compound eyes: read a cube map on a hexagonal grid

Bats: echoes come back at the speed of sound

Nearsight, farsight and presbyopia from the "distance pass"

Careful with flicker

Small gotchas

Main sources

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

What it does

The equations

The camera: an observer swept along with space

Step size

Where I got stuck: a vertical line down the middle

Disk color: the frequency shift

Speed: two levels

Writing the same equations twice, on GPU and CPU

One timetable `TL` holds everything