Fresnel reflectance at an angle, for a conductor or a dielectric¶
The reflectance of unpolarised light arriving from vacuum at a surface of complex refractive index n + ik (lusine.materials.refractive_index at one wavelength), at an angle theta from the normal: the mean of the s and p reflectances of the full Fresnel equations, not Schlick's approximation, written in real arithmetic.
With a^2 + b^2 = sqrt((n^2 - k^2 - sin^2)^2 + 4 n^2 k^2) and a^2 = (a^2 + b^2 + n^2 - k^2 - sin^2) / 2, R_s = (a^2 + b^2 - 2 a cos + cos^2) / (a^2 + b^2 + 2 a cos + cos^2), R_p = R_s (a^2 + b^2 - 2 a sin tan + sin^2 tan^2) / (a^2 + b^2 + 2 a sin tan + sin^2 tan^2), R_p's ratio here multiplied through by cos^2, so it holds at grazing: R = 1 at cos_theta = 0, and lusine.properties. fresnel_f0 at cos_theta = 1. A conductor is k > 0; a dielectric is k = 0, where a = sqrt(n^2 - sin^2) = n cos of the refracted angle. a^2, a sum that cannot be negative, is held at zero where rounding takes it below. Source: Born & Wolf, Principles of Optics, 7th ed., 14.2 (eqs. 14.2.12 to 14.2.15); Pharr, Jakob & Humphreys, Physically Based Rendering, 3rd ed., 8.2.1 (FrConductor). Checked against a reference in numpy, which has its own tests.
- Source
- Born & Wolf, Principles of Optics, 7th ed., 14.2 (eqs. 14.2.12 to 14.2.15); Hecht, Optics, 5th ed., 4.6 and 4.8.
- Proof
- reference, limit, invariant.
fresnel.trait
namespace lusine.properties;
float fresnel(float n, float k, float cos_theta)
{
let c2 = cos_theta * cos_theta;
let s2 = 1.0 - c2;
let n2 = n * n;
let k2 = k * k;
let t = n2 - k2 - s2;
let ab = sqrt(t * t + 4.0 * n2 * k2);
let half = 0.5 * (ab + t);
let a = sqrt(half < 0.0 ? 0.0 : half);
let two_ac = 2.0 * a * cos_theta;
let rs = (ab - two_ac + c2) / (ab + two_ac + c2);
let p_cross = two_ac * s2;
let p_square = ab * c2 + s2 * s2;
let rp = rs * (p_square - p_cross) / (p_square + p_cross);
return 0.5 * (rs + rp);
}
The reflectance of unpolarised light arriving from vacuum at a surface of complex index n + ik (lusine.materials.refractive_index), at an angle θ from the normal, given as cos θ. It is the full Fresnel equations, the mean of the s and p reflectances, written in real arithmetic (Pharr, Jakob & Humphreys's FrConductor, Physically Based Rendering, 3rd ed., 8.2.1), not Schlick's approximation: its only difference from the complex form is binary32's rounding, which this page measures. A dielectric is k = 0.
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.
print(trait_run.source('properties/fresnel.trait'))
// lusine.properties.fresnel: the reflectance of unpolarised light arriving from vacuum at a surface of complex
// refractive index n + ik (lusine.materials.refractive_index at one wavelength), at an angle theta from the normal:
// the mean of the s and p reflectances of the full Fresnel equations, not Schlick's approximation, written in real
// arithmetic. With a^2 + b^2 = sqrt((n^2 - k^2 - sin^2)^2 + 4 n^2 k^2) and a^2 = (a^2 + b^2 + n^2 - k^2 - sin^2) / 2,
// R_s = (a^2 + b^2 - 2 a cos + cos^2) / (a^2 + b^2 + 2 a cos + cos^2),
// R_p = R_s (a^2 + b^2 - 2 a sin tan + sin^2 tan^2) / (a^2 + b^2 + 2 a sin tan + sin^2 tan^2),
// R_p's ratio here multiplied through by cos^2, so it holds at grazing: R = 1 at cos_theta = 0, and lusine.properties.
// fresnel_f0 at cos_theta = 1. A conductor is k > 0; a dielectric is k = 0, where a = sqrt(n^2 - sin^2) = n cos of
// the refracted angle. a^2, a sum that cannot be negative, is held at zero where rounding takes it below.
// Source: Born & Wolf, Principles of Optics, 7th ed., 14.2 (eqs. 14.2.12 to 14.2.15); Pharr, Jakob & Humphreys,
// Physically Based Rendering, 3rd ed., 8.2.1 (FrConductor).
// Checked against a reference in numpy, which has its own tests.
namespace lusine.properties;
float fresnel(float n, float k, float cos_theta)
{
let c2 = cos_theta * cos_theta;
let s2 = 1.0 - c2;
let n2 = n * n;
let k2 = k * k;
let t = n2 - k2 - s2;
let ab = sqrt(t * t + 4.0 * n2 * k2);
let half = 0.5 * (ab + t);
let a = sqrt(half < 0.0 ? 0.0 : half);
let two_ac = 2.0 * a * cos_theta;
let rs = (ab - two_ac + c2) / (ab + two_ac + c2);
let p_cross = two_ac * s2;
let p_square = ab * c2 + s2 * s2;
let rp = rs * (p_square - p_cross) / (p_square + p_cross);
return 0.5 * (rs + rp);
}
The program¶
Each section runs the kernel over a lattice of items, one a case.
print(trait_run.law_program('properties/fresnel.trait', 'fresnel'))
namespace probe;
lusine.version version = { { number = 1 } };
struct item {
host float n;
host float k;
host float cos_theta;
float result_0;
};
lattice<item, 1> items(spacing = 1 m, boundary = wall);
void main() { items.result_0 = lusine.properties.fresnel(items.n, items.k, items.cos_theta); }
Reference: the complex Fresnel equations¶
For every phase in refractive_index at every wavelength, at 91 angles from 0° to 90°, the kernel against the textbook's complex form (Born & Wolf, Principles of Optics, 14.2), computed here in complex float64 from the same binary32 inputs: r_s = (cos θ − N cos θ_t)/(cos θ + N cos θ_t), r_p = (N cos θ − cos θ_t)/(N cos θ + cos θ_t), R = (|r_s|² + |r_p|²)/2, cos θ_t = √(1 − sin²θ/N²). They agree to a few parts in 10⁶ of R: the real form is exact, and the difference is binary32's.
rows = parse_rows(trait_run.source('materials/generated/optical.trait'), 'refractive_index')
theta = np.radians(np.arange(0.0, 91.0, 1.0))
cos = np.float32(np.cos(theta)).astype(float)
cos[-1] = 0.0
n = np.repeat([r['n'] for r in rows], len(cos))
k = np.repeat([r.get('k', 0.0) for r in rows], len(cos))
c = np.tile(cos, len(rows))
got = trait_run.law('properties/fresnel.trait', 'fresnel', n=n, k=k, cos_theta=c)
def exact(n, k, c):
N = np.float32(n).astype(float) + 1j * np.float32(k).astype(float)
ct = np.sqrt(1 - (1 - c * c) / (N * N) + 0j)
rs = (c - N * ct) / (c + N * ct)
rp = (N * c - ct) / (N * c + ct)
return 0.5 * (np.abs(rs) ** 2 + np.abs(rp) ** 2)
want = exact(n, k, c)
worst = float(np.max(np.abs(got - want)))
assert worst < 2e-5, worst
print(f'{len(got)} cases: the largest difference from the complex equations is {worst:.1e} of R')
4095 cases: the largest difference from the complex equations is 3.2e-07 of R
Limit: normal incidence is F0, grazing incidence is everything¶
At θ = 0 the kernel is lusine.properties.fresnel_f0 of the same n and k, to binary32's rounding (their sums differ in order); at θ = 90° it is exactly 1 for every phase.
n = np.array([r['n'] for r in rows])
k = np.array([r.get('k', 0.0) for r in rows])
at0 = trait_run.law('properties/fresnel.trait', 'fresnel', n=n, k=k, cos_theta=1.0)
f0 = trait_run.law('properties/fresnel_f0.trait', 'fresnel_f0', n=n, k=k)
assert np.all(np.abs(at0 - f0) <= 16 * EPS * f0), np.abs(at0 / f0 - 1).max()
at90 = trait_run.law('properties/fresnel.trait', 'fresnel', n=n, k=k, cos_theta=0.0)
assert np.all(at90 == 1.0), at90
print('at 0 degrees, F0 to', np.abs(at0 / f0 - 1).max(), '; at 90 degrees, exactly 1')
at 0 degrees, F0 to 8.228092600592873e-07 ; at 90 degrees, exactly 1
Invariant: Brewster's angle, and the curves¶
A dielectric at Brewster's angle, tan θ = n, reflects no p-polarised light, so the unpolarised reflectance is half of R_s there; glass's (n = 1.5) is 0.0739 at 56.3°. Copper's three channels against glass, beside Schlick's approximation F0 + (1 − F0)(1 − cos θ)⁵ drawn dashed for comparison, which the kernel does not use: for copper Schlick's is off by up to 0.04 of R.
nd = np.array([1.333, 1.5, 2.4])
cb = np.float32(1 / np.sqrt(1 + nd * nd)).astype(float)
got = trait_run.law('properties/fresnel.trait', 'fresnel', n=nd, k=0.0, cos_theta=cb)
ct = np.sqrt(1 - (1 - cb * cb) / nd ** 2)
half_rs = 0.5 * ((cb - nd * ct) / (cb + nd * ct)) ** 2
assert np.all(np.abs(got - half_rs) < 1e-6), got - half_rs
print('at Brewster: ', got, 'half R_s:', half_rs)
cos = np.linspace(0.0, 1.0, 201)
index = {(r['material'], r['wavelength']): (r['n'], r.get('k', 0.0)) for r in rows if r['state'] == 'solid'}
fig, ax = plt.subplots(figsize=(7, 3.6))
worst = 0.0
for w, colour in ((650.0, '#c33'), (550.0, '#3a3'), (450.0, '#33c')):
n_, k_ = index[('copper', w)]
r = trait_run.law('properties/fresnel.trait', 'fresnel', n=n_, k=k_, cos_theta=cos)
f0 = r[-1]
schlick = f0 + (1 - f0) * (1 - cos) ** 5
worst = max(worst, float(np.max(np.abs(schlick - r))))
ax.plot(np.degrees(np.arccos(cos)), r, color=colour, label=f'copper {w:g} nm')
ax.plot(np.degrees(np.arccos(cos)), schlick, '--', color=colour, lw=0.8)
n_, k_ = index[('soda_lime_glass', 589.3)]
glass = trait_run.law('properties/fresnel.trait', 'fresnel', n=n_, k=k_, cos_theta=cos)
ax.plot(np.degrees(np.arccos(cos)), glass, color='k', label='soda-lime glass 589.3 nm')
ax.set_xlabel('angle from the normal (deg)')
ax.set_ylabel('R')
ax.legend()
plt.show()
assert 0.03 < worst < 0.05, worst
print(f"Schlick's approximation for copper is off by up to {worst:.3f}")
at Brewster: [0.03913551 0.07396449 0.24790798] half R_s: [0.03913551 0.0739645 0.24790798]
Schlick's approximation for copper is off by up to 0.040