A directional light¶
A directional light, the sun's beam through the lattice from overhead.
Each cell receives what the cells before it on the beam let through: its shade S is the product of the transmittances of the cells above it, from its neighbour to the lattice's edge where the beam enters (along(-y), y down, row 0 the top, as an image's), and its light E = E0 S, on a surface facing the light. Down the beam that is Beer-Lambert's law cell by cell; an opaque cell (t = 0) shadows everything below it, and a cell is lit on its face whatever its own transmittance. The product is taken nearest first, each multiplication rounded. A beam along no axis is a ray, not yet written. Source: Pharr, Jakob & Humphreys, Physically Based Rendering, 4th ed., 11.2 (transmittance) and 12.4 (distant lights); Liou, An Introduction to Atmospheric Radiation, 2nd ed., 1.4 (the Beer-Bouguer-Lambert law). Checked against a reference in numpy, which has its own tests.
- Source
- Pharr, Jakob & Humphreys, Physically Based Rendering, 4th ed., 11.2 (transmittance) and 12.4 (distant lights); Liou, An Introduction to Atmospheric Radiation, 2nd ed., 1.4 (the Beer-Bouguer-Lambert law).
- Proof
- curve, invariant, limit.
directional.trait
namespace lusine.light;
(float<W/m2>, float) directional(float<W/m2> irradiance, float transmittance)
{
let shade = product(along(-y).transmittance);
return (irradiance * shade, shade);
}
The sun's beam through the lattice from overhead, each cell receiving what the cells above it let through. Its shade S is the product of their transmittances, product(along(-y).transmittance) (from the neighbour to the lattice's edge, nearest first, row 0 the top), the shadow map, and its light E = E0 S. A beam along no axis is a ray, which this kernel does not cover yet.
import os, sys
root = os.getcwd()
while not os.path.isfile(os.path.join(root, 'docsite.json')):
root = os.path.dirname(root)
sys.path.insert(0, os.path.join(root, 'tools'))
import numpy as np
import matplotlib.pyplot as plt
from lusine import trait_run
F32 = np.float32
U = float(np.finfo(F32).eps) / 2 # binary32's unit roundoff, 2^-24
The kernel¶
The law as the library states it, and as this page runs it.
print(trait_run.source('light/directional.trait'))
// lusine.light.directional: a directional light, the sun's beam through the lattice from overhead. Each cell receives
// what the cells before it on the beam let through: its shade S is the product of the transmittances of the cells
// above it, from its neighbour to the lattice's edge where the beam enters (along(-y), y down, row 0 the top, as an
// image's), and its light E = E0 S, on a surface facing the light. Down the beam that is
// Beer-Lambert's law cell by cell; an opaque cell (t = 0) shadows everything below it, and a cell is lit on its face
// whatever its own transmittance. The product is taken nearest first, each multiplication rounded. A beam along no
// axis is a ray, not yet written.
// Source: Pharr, Jakob & Humphreys, Physically Based Rendering, 4th ed., 11.2 (transmittance) and 12.4 (distant
// lights); Liou, An Introduction to Atmospheric Radiation, 2nd ed., 1.4 (the Beer-Bouguer-Lambert law).
// Checked against a reference in numpy, which has its own tests.
namespace lusine.light;
(float<W/m2>, float) directional(float<W/m2> irradiance, float transmittance)
{
let shade = product(along(-y).transmittance);
return (irradiance * shade, shade);
}
The program¶
A plane of cells, x across and y down, each with its transmittance, the host's; the sun's irradiance is a setting, and main() binds the kernel over the cells. sunlit(transmittance, E0) runs it once for transmittances laid out [x, y], giving the light and the shade. The planted faults below run the same program with the page's own wrong kernel in the library's place.
PROGRAM = '''namespace sun;
lusine.version version = { { number = 1 } };
struct cell { host float transmittance; float<W/m2> light; float shade; };
lattice<cell, 2> board(spacing = 1 m, boundary = wall);
host struct { float<W/m2> irradiance; } settings;
%s
void main() { (board.light, board.shade) = %s(settings.irradiance, board.transmittance); }
'''
LAW = PROGRAM % ('', 'lusine.light.directional')
print(LAW)
def sunlit(transmittance, irradiance=1000.0, program=LAW):
t = np.asarray(transmittance, dtype=float)
run = trait_run.Program(['light/directional.trait'], program, {'sun.board': list(t.shape)})
run.store('sun.board', 'transmittance', t)
run.store('sun.settings', 'irradiance', irradiance)
light, shade = run.steps('sun.main', 1, [('sun.board', 'light'), ('sun.board', 'shade')])
return light[0], shade[0]
namespace sun;
lusine.version version = { { number = 1 } };
struct cell { host float transmittance; float<W/m2> light; float shade; };
lattice<cell, 2> board(spacing = 1 m, boundary = wall);
host struct { float<W/m2> irradiance; } settings;
void main() { (board.light, board.shade) = lusine.light.directional(settings.irradiance, board.transmittance); }
Curve: Beer-Lambert down the beam¶
Under a sun overhead, through a medium letting 1 − m through each cell, the light k cells down is E0 (1 − m)^k, the exponential decay of Beer-Bouguer-Lambert cell by cell, to binary32's rounding of the k multiplications.
E0 = 1000.0
depth = np.arange(40)
def beer_holds(light, m):
want = E0 * float(F32(1.0 - m)) ** depth
return bool(np.all(np.abs(light - want[None, :]) <= (depth + 2) * U * want[None, :]))
fig, ax = plt.subplots(figsize=(7, 3.2))
for m in (0.02, 0.05, 0.1):
E, S = sunlit(np.full((5, 40), 1.0 - m), E0)
assert beer_holds(E, m), m
ax.semilogy(depth, E[2], 'o', ms=3, label=f'kernel, m = {m}')
ax.semilogy(depth, E0 * (1.0 - m) ** depth, lw=0.8)
ax.set_xlabel('depth (cells)')
ax.set_ylabel('light (W/m2)')
ax.legend()
plt.show()
Reference: bit for bit, the shade in binary32¶
Random transmittances, 9 columns of 50 cells: the kernel's shade and light are a binary32 model's bit for bit, the transmittances above a cell multiplied nearest first, each multiplication rounded, and the light E0 times the shade. The numpy reference in float64 (its rays down the second axis here) agrees within 50 × 2⁻²⁴ relative, a product of 50 roundings; where it underflows binary32, the shade is 0 or a denormal.
sys.path.insert(0, os.path.join(root, 'kernels'))
from lusine.kernels.light.directional import directional as reference
def model(transmittance, irradiance):
t = F32(transmittance)
shade = np.empty(t.shape, dtype=F32)
for x in range(t.shape[0]):
for y in range(t.shape[1]):
s = F32(1.0)
for j in range(y - 1, -1, -1): # along(-y): the neighbour first, the edge last
s = F32(s * t[x, j])
shade[x, y] = s
return F32(F32(irradiance) * shade), shade
rng = np.random.default_rng(93)
t = F32(rng.uniform(0.6, 1.0, (9, 50))).astype(float)
E, S = sunlit(t, E0)
mE, mS = model(t, E0)
assert np.array_equal(F32(E), mE) and np.array_equal(F32(S), mS)
rE, rS = reference(E0, (0, 1), t)
normal = rS > 1e-30
rel = np.max(np.abs(S - rS)[normal] / rS[normal])
assert rel <= 50 * U, rel / U
print(f'{S.size} cells bit for bit with the binary32 model; the float64 reference within {rel / U:.1f} x 2^-24')
450 cells bit for bit with the binary32 model; the float64 reference within 8.3 x 2^-24
Invariant: the shade only falls down a beam¶
The shade is a share of the open sky's, between 0 and 1, and going down a column it never grows: each cell only takes from what it lets on, a multiplication by at most 1, which binary32 cannot round up past what it was. Each column is its own: changing one leaves every other bit for bit.
t = rng.random((30, 30))
_, S = sunlit(t)
assert np.all((S >= 0) & (S <= 1)) and np.all(np.diff(S, axis=1) <= 0)
changed = t.copy()
changed[7, :10] = 0.0
_, R = sunlit(changed)
assert np.array_equal(np.delete(R, 7, axis=0), np.delete(S, 7, axis=0)) and np.all(R[7, 10:] == 0)
print('a share, never growing down a beam; the columns independent')
a share, never growing down a beam; the columns independent
Limit: an opaque cell's shadow, and an open sky¶
Under an open sky every cell has the sun's irradiance exactly. An opaque cell (t = 0) is lit on its face, whatever its own transmittance, and casts its shadow straight down, every cell below it exactly dark and nothing beside it touched.
E, S = sunlit(np.ones((20, 20)), 1000.0)
assert np.all(E == 1000.0) and np.all(S == 1.0)
def shadow_holds(S):
return bool(np.all(S[10:16, :9] == 1.0) and np.all(S[10:16, 9] == 1.0) and np.all(S[10:16, 12:] == 0.0)
and np.all(S[:10] == 1.0) and np.all(S[16:] == 1.0))
t = np.ones((40, 40))
t[10:16, 9:12] = 0.0
_, S = sunlit(t)
assert shadow_holds(S)
fig, ax = plt.subplots(1, 2, figsize=(8, 3.6))
ax[0].imshow(t.T, cmap='gray', vmin=0, vmax=1); ax[0].set_title('transmittance')
ax[1].imshow(S.T, cmap='gray', vmin=0, vmax=1); ax[1].set_title('shade, the sun overhead')
plt.show()
Planted faults¶
Two wrong kernels, written here in Trait and run through the same program in the library's place: the light taken from below (along(+y), the beam walked against the light), and the cell's own transmittance multiplied in, so an opaque cell is dark on its own face. Each fails the opaque cell's shadow and the binary32 model.
BELOW = PROGRAM % ('''
(float<W/m2>, float) beneath(float<W/m2> irradiance, float transmittance)
{
let shade = product(along(+y).transmittance);
return (irradiance * shade, shade);
}
''', 'beneath')
OWN = PROGRAM % ('''
(float<W/m2>, float) own(float<W/m2> irradiance, float transmittance)
{
let shade = product(along(-y).transmittance) * transmittance;
return (irradiance * shade, shade);
}
''', 'own')
t2 = F32(rng.uniform(0.6, 1.0, (9, 50))).astype(float)
for name, wrong in (('the beam walked from below', BELOW), ('the cell its own shadow', OWN)):
caught = [not shadow_holds(sunlit(t, program=wrong)[1]),
not np.array_equal(F32(sunlit(t2, E0, wrong)[1]), model(t2, E0)[1])]
assert all(caught), name
print(f'PLANTED, {name}: the shadow and the model each catch it')
PLANTED, the beam walked from below: the shadow and the model each catch it PLANTED, the cell its own shadow: the shadow and the model each catch it