The absorption coefficient: 4πk over the wavelength¶

The attenuation coefficient of light in a phase, alpha = 4 pi k / lambda, from its extinction coefficient k at a wavelength lambda in vacuum (lusine.materials.refractive_index): the intensity falls as exp(-alpha d) over a depth d (Beer-Lambert), 1 / alpha the depth it falls by e.

Water's red goes in metres and its blue in tens of metres; a transparent phase (k absent, none) has none. Source: Born & Wolf, Principles of Optics, 7th ed., 14.1 (eqs. 14.1.18 to 14.1.24); Hecht, Optics, 5th ed., 4.8. Checked against a reference in numpy, which has its own tests.

Source

Born & Wolf, Principles of Optics, 7th ed., 14.1 (eqs. 14.1.18 to 14.1.24); Hecht, Optics, 5th ed., 4.8.

Proof

reference, invariant.

absorption.trait

namespace lusine.properties;

float<m^-1> absorption(float k, float<m> wavelength)
{
  return 4.0 * 3.141592653589793 * k / wavelength;
}

The attenuation coefficient α = 4πk/λ of light in a phase, from its extinction coefficient k at a wavelength λ in vacuum (lusine.materials.refractive_index); the intensity falls as e^(−αd) over a depth d (Beer–Lambert). The library stores no absorption coefficient: this kernel derives it.

In [1]:

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
from lusine.library_sources import parse_rows
EPS = float(np.finfo(np.float32).eps)   # Trait's float is binary32

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 from lusine.library_sources import parse_rows EPS = float(np.finfo(np.float32).eps) # Trait's float is binary32

The kernel¶

The law as the library states it.

In [2]:

print(trait_run.source('properties/absorption.trait'))

print(trait_run.source('properties/absorption.trait'))

// lusine.properties.absorption: the attenuation coefficient of light in a phase, alpha = 4 pi k / lambda, from its
// extinction coefficient k at a wavelength lambda in vacuum (lusine.materials.refractive_index): the intensity falls
// as exp(-alpha d) over a depth d (Beer-Lambert), 1 / alpha the depth it falls by e. Water's red goes in metres and
// its blue in tens of metres; a transparent phase (k absent, none) has none.
// Source: Born & Wolf, Principles of Optics, 7th ed., 14.1 (eqs. 14.1.18 to 14.1.24); Hecht, Optics, 5th ed., 4.8.
// Checked against a reference in numpy, which has its own tests.

namespace lusine.properties;

float<m^-1> absorption(float k, float<m> wavelength)
{
  return 4.0 * 3.141592653589793 * k / wavelength;
}

The program¶

Each section runs the kernel over a lattice of items, one a case.

In [3]:

print(trait_run.law_program('properties/absorption.trait', 'absorption'))

print(trait_run.law_program('properties/absorption.trait', 'absorption'))

namespace probe;

lusine.version version = { { number = 1 } };

struct item {
  host float k;
  host float<m> wavelength;
  float<m^-1> result_0;
};
lattice<item, 1> items(spacing = 1 m, boundary = wall);

void main() { items.result_0 = lusine.properties.absorption(items.k, items.wavelength); }

Reference: every phase in the library¶

Each row of refractive_index, as committed (its wavelength in nm, so ×10⁻⁹ to metres), is held to 4πk/λ, to binary32's rounding. A row with no k (boxwood's single hand-kept index, alpha-quartz's Sellmeier index, a grass blade's PROSPECT-5 index) reads k = 0 and absorbs nothing: α = 0, exactly. Soda-lime window glass (Rubin's clear float glass) absorbs weakly, by e over 20 cm of green (550 nm) and 6.2 cm of red (650 nm), its iron's tint seen edge-on. Liquid water falls by e over 3.2 m of red (650 nm) and 35 m of blue (450 nm), why deep water is blue; silicon's blue goes in 0.4 µm.

In [4]:

rows = parse_rows(trait_run.source('materials/generated/optical.trait'), 'refractive_index')
k = np.array([r.get('k', 0.0) for r in rows])
lam = np.array([r['wavelength'] for r in rows]) * 1e-9
a = trait_run.law('properties/absorption.trait', 'absorption', k=k, wavelength=lam)
want = 4 * np.pi * np.float32(k).astype(float) / np.float32(lam).astype(float)
assert np.all(np.abs(a - want) <= 4 * EPS * want), np.abs(a / want - 1).max()
at = {(r['material'], r['state'], r['wavelength']): v for r, v in zip(rows, a)}
assert at[('boxwood', 'solid', 589.3)] == 0.0 and at[('silicon_dioxide', 'solid', 550.0)] == 0.0 and at[('grass', 'solid', 550.0)] == 0.0
assert abs(1 / at[('soda_lime_glass', 'solid', 550.0)] - 0.1989) < 0.001, 1 / at[('soda_lime_glass', 'solid', 550.0)]
red, blue = at[('water', 'liquid', 650.0)], at[('water', 'liquid', 450.0)]
assert abs(1 / red - 3.154) < 0.002 and abs(1 / blue - 35.11) < 0.02, (1 / red, 1 / blue)
si = at[('silicon', 'solid', 450.0)]
assert abs(1 / si - 0.415e-6) < 0.002e-6, 1 / si
for key in (('water', 'liquid', 650.0), ('water', 'liquid', 450.0), ('silicon', 'solid', 450.0), ('soda_lime_glass', 'solid', 550.0)):
    v = at[key]
    print(f'{key[0]:15s} {key[1]:7s} {key[2]:5g} nm: alpha {v:.4g} /m' + (f', 1/alpha {1 / v:.4g} m' if v else ''))

rows = parse_rows(trait_run.source('materials/generated/optical.trait'), 'refractive_index') k = np.array([r.get('k', 0.0) for r in rows]) lam = np.array([r['wavelength'] for r in rows]) * 1e-9 a = trait_run.law('properties/absorption.trait', 'absorption', k=k, wavelength=lam) want = 4 * np.pi * np.float32(k).astype(float) / np.float32(lam).astype(float) assert np.all(np.abs(a - want) <= 4 * EPS * want), np.abs(a / want - 1).max() at = {(r['material'], r['state'], r['wavelength']): v for r, v in zip(rows, a)} assert at[('boxwood', 'solid', 589.3)] == 0.0 and at[('silicon_dioxide', 'solid', 550.0)] == 0.0 and at[('grass', 'solid', 550.0)] == 0.0 assert abs(1 / at[('soda_lime_glass', 'solid', 550.0)] - 0.1989) < 0.001, 1 / at[('soda_lime_glass', 'solid', 550.0)] red, blue = at[('water', 'liquid', 650.0)], at[('water', 'liquid', 450.0)] assert abs(1 / red - 3.154) < 0.002 and abs(1 / blue - 35.11) < 0.02, (1 / red, 1 / blue) si = at[('silicon', 'solid', 450.0)] assert abs(1 / si - 0.415e-6) < 0.002e-6, 1 / si for key in (('water', 'liquid', 650.0), ('water', 'liquid', 450.0), ('silicon', 'solid', 450.0), ('soda_lime_glass', 'solid', 550.0)): v = at[key] print(f'{key[0]:15s} {key[1]:7s} {key[2]:5g} nm: alpha {v:.4g} /m' + (f', 1/alpha {1 / v:.4g} m' if v else ''))

water           liquid    650 nm: alpha 0.3171 /m, 1/alpha 3.154 m
water           liquid    450 nm: alpha 0.02848 /m, 1/alpha 35.11 m
silicon         solid     450 nm: alpha 2.41e+06 /m, 1/alpha 4.149e-07 m
soda_lime_glass solid     550 nm: alpha 5.027 /m, 1/alpha 0.1989 m

Invariant: water's colour through a depth¶

What reaches a depth d of liquid water, e^(−αd), at the four wavelengths: by 10 m the red is gone and the blue nearly all there. αλ/k is 4π for every phase and wavelength.

In [5]:

water = [r for r in rows if r['material'] == 'water' and r['state'] == 'liquid']
a = trait_run.law('properties/absorption.trait', 'absorption', k=[r['k'] for r in water], wavelength=[r['wavelength'] * 1e-9 for r in water])
d = np.linspace(0.0, 20.0, 201)
fig, ax = plt.subplots(figsize=(7, 3.2))
for r, v in zip(water, a):
    ax.plot(d, np.exp(-v * d), label=f"{r['wavelength']:g} nm")
ax.set_xlabel('depth (m)')
ax.set_ylabel('transmitted')
ax.legend()
plt.show()
assert np.exp(-a[-1] * 10) < 0.05 and np.exp(-a[0] * 10) > 0.7
rng = np.random.default_rng(325)
kk = np.float32(10 ** rng.uniform(-10, 1, 32)).astype(float)
lam = np.float32(rng.uniform(300e-9, 800e-9, 32)).astype(float)
a = trait_run.law('properties/absorption.trait', 'absorption', k=kk, wavelength=lam)
assert np.all(np.abs(a * lam / kk / (4 * np.pi) - 1) <= 4 * EPS)

water = [r for r in rows if r['material'] == 'water' and r['state'] == 'liquid'] a = trait_run.law('properties/absorption.trait', 'absorption', k=[r['k'] for r in water], wavelength=[r['wavelength'] * 1e-9 for r in water]) d = np.linspace(0.0, 20.0, 201) fig, ax = plt.subplots(figsize=(7, 3.2)) for r, v in zip(water, a): ax.plot(d, np.exp(-v * d), label=f"{r['wavelength']:g} nm") ax.set_xlabel('depth (m)') ax.set_ylabel('transmitted') ax.legend() plt.show() assert np.exp(-a[-1] * 10) < 0.05 and np.exp(-a[0] * 10) > 0.7 rng = np.random.default_rng(325) kk = np.float32(10 ** rng.uniform(-10, 1, 32)).astype(float) lam = np.float32(rng.uniform(300e-9, 800e-9, 32)).astype(float) a = trait_run.law('properties/absorption.trait', 'absorption', k=kk, wavelength=lam) assert np.all(np.abs(a * lam / kk / (4 * np.pi) - 1) <= 4 * EPS)