"""What a spotlight makes of the cut disc, and what lands on the wall. The disc is a mirror with grooves in it. That is the whole optical situation, and it has one governing fact. The cone -------- A groove is PRISMATIC: the tip sweeps its own shape along the direction of travel, so the groove is invariant along its own axis t. That includes the bottom — a sphere dragged along a line sweeps a CYLINDER, not a sphere — so every surface normal in a groove is perpendicular to t. Reflection is r = i - 2(i.n)n, and if n.t = 0 then r.t = i.t - 2(i.n)(n.t) = i.t The component along the groove is conserved. So every ray leaving a groove of direction t lies on a cone about t whose half-angle is the angle the incoming ray already made with t. A groove cannot scatter light anywhere it likes; it can only move it around that cone. On a flat wall a cone draws a conic section, so each groove direction contributes an arc, and the picture on the wall is a superposition of arcs weighted by how much groove runs in each direction. This is why brushed metal has a streak instead of a highlight, and it is why the projection carries the machine's handwriting: the arcs are set by the directions the machine chose to travel. Numbers for this tip, from stylus.py ------------------------------------ Groove depth 26.93 um, width 100 um (specified; the load follows). The wall of the 120 degree cone stands at 30 degrees to the surface, so its normal is 30 degrees off vertical, and a normally incident ray leaves it at 60 degrees from vertical. At exactly normal incidence that reflected ray runs PARALLEL to the opposite wall — a 120 degree V neither retroreflects (as a 90 degree V would) nor traps the ray. Anything else that matters is computed here, not assumed. Because the groove is prismatic, whether a ray escapes or hits a second facet is decided entirely by the two-dimensional cross-section, using only the components of the ray across and normal to the groove. So shadowing and second bounces are exact here, not approximated. Self-test: python3 optics.py --self-test """ from __future__ import annotations import argparse import json import math import os import numpy as np import stylus as st # --------------------------------------------------------------------------- # the metal # --------------------------------------------------------------------------- # Polished 304 in the visible, from a published dispersion fit for 304 (Castelli # et al.): n and k both rise monotonically toward the red, giving a normal # incidence reflectance of about 62% and a slightly warm tint. Note that for a # metal, unpolarised reflectance is FLAT or very slightly falling from 0 to 70 # degrees and only climbs toward unity in the last degree or two — it does not # rise steadily toward grazing. STEEL_N = {"r": 1.920, "g": 1.640, "b": 1.360} # 650 / 550 / 450 nm STEEL_K = {"r": 3.554, "g": 3.224, "b": 2.894} # Roughness as the standard deviation of surface slope, in radians. # # The groove floor is the important one and it is strongly ANISOTROPIC: the tip # smears metal along its own path, so the floor is nearly smooth along the # groove and rough across it (estimated slope 0.65 deg along, 16 deg across — # a ratio of about 25:1). # # That anisotropy is aligned with the one direction that matters. A normal # perturbed WITHIN the cross-section still satisfies n.t = 0, so the reflected # ray stays on the same cone: across-groove roughness moves light ALONG the arc # it already draws. Only the small along-groove component tilts n out of that # plane and takes light OFF the cone, blurring the arc across its width. So the # arcs survive roughness far better than an isotropic lobe would suggest. SLOPE_MIRROR = math.radians(0.30) # No.8; see the spec-band caveat below SLOPE_GROOVE_ACROSS = math.radians(16.0) # smears along the arc SLOPE_GROOVE_ALONG = math.radians(0.65) # blurs the arc across its width SLOPE_RIDGE = math.radians(18.0) # torn metal, near-diffuse # "No. 8 mirror finish" is a spec band, not one optical state: across its own # Ra 0.01-0.05 um range the true specular spike falls from about 90% to about # 10% of the reflected energy. MIRROR_SPIKE says how much of the polished land's # reflection stays in a sharp lamp image; the rest becomes a wide haze. MIRROR_SPIKE = 0.55 # Hemispherical reflectance relative to the polished land. Roughening a metal # redistributes light rather than absorbing it; the loss is multiple bounces # inside the microtexture. The SPECULAR spike of both the groove floor and the # ridge is effectively zero (their roughness is many wavelengths), which is why # these are hemispherical factors and the lobe widths above carry the shape. REFL_SCALE_GROOVE = 0.90 REFL_SCALE_RIDGE = 0.75 def fresnel_metal(cos_i, n, k): """Unpolarised reflectance of a metal at incidence angle acos(cos_i).""" c = np.clip(np.abs(cos_i), 1e-6, 1.0) s2 = 1.0 - c * c n2k2 = n * n - k * k - s2 a2b2 = np.sqrt(np.maximum(n2k2 * n2k2 + 4.0 * n * n * k * k, 0.0)) a2 = 0.5 * (a2b2 + n2k2) a = np.sqrt(np.maximum(a2, 0.0)) Rs = (a2b2 - 2.0 * a * c + c * c) / (a2b2 + 2.0 * a * c + c * c) t = s2 / c Rp = Rs * (a2b2 - 2.0 * a * t + t * t) / (a2b2 + 2.0 * a * t + t * t) return np.clip(0.5 * (Rs + Rp), 0.0, 1.0) def steel_rgb(cos_i): return np.stack([fresnel_metal(cos_i, STEEL_N[c], STEEL_K[c]) for c in ("r", "g", "b")], axis=-1) # --------------------------------------------------------------------------- # the groove cross-section # --------------------------------------------------------------------------- def cross_section(n=241, load_gf=st.LOAD_GF, with_ridge=True): """The profile a single pass leaves, sampled along its arc. Returns (u, z, nu, nz, ds, kind) with u across the groove and z up, in micrometres. kind: 0 mirror (flat, outside everything), 1 groove floor, 2 ridge. Normals point out of the metal. """ d = st.depth_for_load_um(st.HV_ANNEALED, load_gf) a_c = st.contact_radius_um(d) A0, h0, R = st.A0_UM, st.H0_UM, st.TIP_R_UM ridge_w = st.PILEUP_WIDTH * 2.0 * a_c peak = st.pileup_height_um(d) span = a_c + (ridge_w if with_ridge else 0.0) + 6.0 u = np.linspace(-span, span, n) au = np.abs(u) z = np.zeros_like(u) kind = np.zeros(u.shape, dtype=np.int8) # inside the groove: the swept tip. Cylinder of radius R below the # crossover, straight wall above it. ins = au <= a_c zc = -d + R # cylinder axis height cyl = ins & (au <= A0) z[cyl] = zc - np.sqrt(np.maximum(R * R - u[cyl] ** 2, 0.0)) wall = ins & (au > A0) z[wall] = (-d + h0) + (au[wall] - A0) / st.TAN_A kind[ins] = 1 if with_ridge: # continuous: 0 at the rim, crest a little outside it, 0 again at the far # edge. A step at the rim would put a near-vertical wall 7.9 um high # right next to the groove and block almost everything leaving it. rg = (au > a_c) & (au <= a_c + ridge_w) x = (au[rg] - a_c) / ridge_w c = st.PILEUP_CREST z[rg] = np.where(x < c, peak * x / c, peak * (1.0 - x) / (1.0 - c)) kind[rg] = 2 # outward normal from the profile slope dz = np.gradient(z, u) nrm = np.sqrt(1.0 + dz * dz) nu, nz = -dz / nrm, 1.0 / nrm ds = np.gradient(u) * nrm # arc length per sample return u, z, nu, nz, ds, kind, d, a_c class Profile: """The cross-section, plus exact 2D visibility inside it.""" def __init__(self, n=241, load_gf=st.LOAD_GF, with_ridge=True): (self.u, self.z, self.nu, self.nz, self.ds, self.kind, self.depth, self.half_width) = cross_section(n, load_gf, with_ridge) self.period = 2.0 * (self.u[-1]) # one groove, isolated def horizon(self, eps=1e-9): """Per sample, the angular window through which the sky is visible. The profile is a single-valued valley, so from any point on it the unobstructed directions form ONE angular interval, bounded by the highest thing to the left and the highest thing to the right. Writing an upward direction as beta = atan2(dz, du) in (0, pi), a ray escapes exactly when beta_right < beta < beta_left Precomputing those two bounds turns the visibility test into a scalar comparison, which is what makes the whole-disc integration tractable — and it is exact, not an approximation, because the groove does not vary along its own axis. """ U, Z = self.u, self.z n = U.size dU = U[None, :] - U[:, None] dZ = Z[None, :] - Z[:, None] ang = np.arctan2(dZ, dU) # (-pi, pi] right = dU > eps left = dU < -eps # to the right: the largest elevation blocks everything below it aR = np.where(right, ang, -np.inf).max(axis=1) aR = np.maximum(aR, 0.0) # never below the horizontal # to the left: angles are in (pi/2, pi); the smallest one blocks above it angL = np.where(left, np.mod(ang, 2 * np.pi), np.inf) aL = angL.min(axis=1) aL = np.minimum(np.where(np.isfinite(aL), aL, np.pi), np.pi) return aR, aL def escapes(self, du, dz): """Vectorised escape test. du, dz broadcast against the sample axis.""" if not hasattr(self, "_aR"): self._aR, self._aL = self.horizon() beta = np.arctan2(dz, du) return (dz > 0) & (beta > self._aR) & (beta < self._aL) def escapes_at(self, j, du, dz): """Escape test for one sample index j, vectorised over many directions.""" if not hasattr(self, "_aR"): self._aR, self._aL = self.horizon() beta = np.arctan2(dz, du) return (dz > 0) & (beta > self._aR[j]) & (beta < self._aL[j]) def visible(self, dir_u, dir_z, eps=1e-4): """Can a ray of 2D direction (dir_u, dir_z) leave each sample point? Because the groove is invariant along its own axis, this test in the cross-section is exact for any 3D ray: only the across and vertical components decide whether the profile is in the way. Returns a boolean per sample. Vectorised over samples for one direction. """ du, dz = np.asarray(dir_u, dtype=float), np.asarray(dir_z, dtype=float) m = np.hypot(du, dz) if m == 0: return np.zeros(self.u.shape, dtype=bool) du, dz = du / m, dz / m if dz <= 0: # heading into the metal return np.zeros(self.u.shape, dtype=bool) # march in u, comparing ray height against the profile U, Z = self.u, self.z ok = np.ones(U.shape, dtype=bool) if abs(du) < 1e-9: return ok # straight up: always clear slope = dz / du # ray from (U[j], Z[j]) has height Z[j] + slope*(U-U[j]) at abscissa U ray = Z[:, None] + slope * (U[None, :] - U[:, None]) ahead = ((U[None, :] - U[:, None]) * du) > 0.0 blocked = ahead & (ray < (Z[None, :] - eps)) return ~blocked.any(axis=1) # --------------------------------------------------------------------------- # reflection off a groove of a given direction # --------------------------------------------------------------------------- def groove_reflect(prof, i_vec, t_vec, up=None, cos_gate=1e-4, bounces=2): """Scatter one incident direction off one groove direction. i_vec: unit vector the light travels along, pointing at the surface. t_vec: unit vector along the groove, lying in the surface. up: the SURFACE normal. Defaults to +z, which is only right for a horizontal surface — pass the actual normal for a tilted disc, or the groove gets built in the wrong plane. (The cone property survives either way, since the normals stay perpendicular to t, which is why this was easy to miss.) Returns (dirs, weights, kinds): outgoing unit vectors, the flux each carries per unit groove length (in micrometres of projected width), and which kind of surface produced it. Shadowing and a second bounce are both resolved in the cross-section, exactly. """ i_vec = np.asarray(i_vec, float) i_vec = i_vec / np.linalg.norm(i_vec) t = np.asarray(t_vec, float) t = t / np.linalg.norm(t) up = np.array([0.0, 0.0, 1.0]) if up is None else ( np.asarray(up, float) / np.linalg.norm(up)) t = t - (t @ up) * up # keep t in the surface t = t / np.linalg.norm(t) u_hat = np.cross(t, up) u_hat /= np.linalg.norm(u_hat) # across the groove, in-plane # incident direction resolved onto (u, t, z) iu, it, iz = i_vec @ u_hat, i_vec @ t, i_vec @ up # is each sample lit? the incoming ray arrives along -i, so it can reach a # sample if a ray leaving along -i escapes. lit = prof.visible(-iu, -iz) # 3D normals from the 2D profile N = prof.nu[:, None] * u_hat[None, :] + prof.nz[:, None] * up[None, :] cos_in = -(N @ i_vec) # >0 where facing the light face = cos_in > cos_gate live = lit & face dirs = np.zeros((prof.u.size, 3)) w = np.zeros(prof.u.size) kinds = prof.kind.copy() if not live.any(): return dirs[:0], w[:0], kinds[:0] idx = np.nonzero(live)[0] n_live = N[idx] r = i_vec[None, :] - 2.0 * (n_live @ i_vec)[:, None] * n_live # projected width of each facet as seen by the light, per unit groove length weight = cos_in[idx] * prof.ds[idx] out_dirs = r.copy() out_w = weight.copy() out_k = prof.kind[idx].copy() if bounces > 1: # does the reflected ray escape? test in the cross-section ru = out_dirs @ u_hat rz = out_dirs @ up esc = np.ones(len(idx), dtype=bool) for j in range(len(idx)): v = prof.visible(ru[j], rz[j]) esc[j] = bool(v[idx[j]]) # rays that do not escape hit the opposite facet: reflect once more off # the mirror image of the profile. For a symmetric V the second facet's # normal is the first one's with nu negated. second = ~esc if second.any(): n2 = n_live[second].copy() n2 = n2 - 2.0 * (n2 @ u_hat)[:, None] * u_hat[None, :] d2 = out_dirs[second] d2 = d2 - 2.0 * (d2 @ n2.T).diagonal()[:, None] * n2 out_dirs[second] = d2 # a second bounce costs another reflectance out_w[second] *= 0.55 # keep only what is finally heading away from the surface good = (out_dirs @ up) > 1e-6 return out_dirs[good], out_w[good], out_k[good] # --------------------------------------------------------------------------- # self-test # --------------------------------------------------------------------------- def self_test(): ok = True def check(name, got, want, tol): nonlocal ok good = abs(got - want) <= tol ok &= good print(" %-52s %10.4f want %8.4f %s" % (name, got, want, "ok" if good else "FAIL")) print("groove geometry") prof = Profile() d, a_c = prof.depth, prof.half_width check("depth, um", d, 26.934, 0.01) check("half width, um", a_c, 50.0, 0.05) # the blunting sphere must be TANGENT to the cone, not merely meeting it: # getting this wrong put a 30 degree kink here and the depth 1.93 um too deep Rt = st.TIP_R_UM arc_slope = st.A0_UM / math.sqrt(max(Rt * Rt - st.A0_UM ** 2, 1e-12)) check("kink at the sphere/cone junction, deg", abs(math.degrees(math.atan(arc_slope)) - math.degrees(math.atan(1.0 / st.TAN_A))), 0.0, 1e-6) check("width / depth (blunted 120 deg cone)", a_c * 2.0 / d, 3.71, 0.01) # the cut model and the trajectories must describe the same tool import scribe as _sc check("stylus width vs scribe.GROOVE, mm", st.GROOVE_MM - _sc.GROOVE, 0.0, 1e-12) # the straight wall: 30 deg to the horizontal, normal 30 deg off vertical wall = (np.abs(prof.u) > st.A0_UM + 1.0) & (np.abs(prof.u) < a_c - 1.0) ang = np.degrees(np.arccos(prof.nz[wall])) check("wall normal tilt from vertical, deg", float(np.median(ang)), 30.0, 0.6) print("the cone property: r.t must equal i.t for every facet") worst = 0.0 rng = np.random.default_rng(3) for _ in range(40): v = rng.normal(size=3) v[2] = -abs(v[2]) - 0.15 i_vec = v / np.linalg.norm(v) a = rng.uniform(0, np.pi) t = np.array([math.cos(a), math.sin(a), 0.0]) dirs, w, k = groove_reflect(prof, i_vec, t, bounces=1) if len(dirs) == 0: continue worst = max(worst, float(np.abs(dirs @ t - i_vec @ t).max())) check("max |r.t - i.t| over 40 random cases", worst, 0.0, 1e-9) print("normal incidence off the straight wall") i_vec = np.array([0.0, 0.0, -1.0]) t = np.array([0.0, 1.0, 0.0]) # groove along y up = np.array([0.0, 0.0, 1.0]) u_hat = np.cross(t, up); u_hat /= np.linalg.norm(u_hat) n2 = np.array([-0.5, 0.0, math.sqrt(3) / 2]) # the +u wall's normal r = i_vec - 2.0 * (n2 @ i_vec) * n2 check("reflected angle from vertical, deg", float(np.degrees(np.arccos(r @ up))), 60.0, 1e-6) # and it should run parallel to the opposite wall opp = np.array([-math.cos(math.radians(30)), 0.0, math.sin(math.radians(30))]) check("|sin| between it and the opposite wall", float(abs(np.cross(r, opp)[1])), 0.0, 1e-6) print("a flat mirror must reflect specularly and lose nothing") flat = Profile(n=41, with_ridge=False) flat.z[:] = 0.0 flat.nu[:] = 0.0 flat.nz[:] = 1.0 for ang_deg in (0.0, 30.0, 60.0): a = math.radians(ang_deg) i_vec = np.array([math.sin(a), 0.0, -math.cos(a)]) dirs, w, k = groove_reflect(flat, i_vec, np.array([0.0, 1.0, 0.0]), bounces=1) want = np.array([math.sin(a), 0.0, math.cos(a)]) err = float(np.abs(dirs - want[None, :]).max()) if len(dirs) else 9.9 check("specular error at %2.0f deg incidence" % ang_deg, err, 0.0, 1e-12) print("Fresnel reflectance of polished 304") # A metal does NOT approach unity by 85 degrees: p-polarised reflectance is # still near its pseudo-Brewster minimum there, so the unpolarised mean dips # to about 0.70 and only climbs to 1 in the last degree. Checked by hand # against the a,b parametrisation before trusting it. for ang_deg, lo, hi in ((0.0, 0.50, 0.72), (45.0, 0.50, 0.78), (85.0, 0.60, 0.80), (89.8, 0.95, 1.0)): R = steel_rgb(np.array([math.cos(math.radians(ang_deg))]))[0] m = float(R.mean()) good = lo <= m <= hi ok &= good print(" %-52s %10.4f want %.2f-%.2f %s" % ("mean reflectance at %4.1f deg" % ang_deg, m, lo, hi, "ok" if good else "FAIL")) # and it must be monotonic once past the p-minimum angs = np.radians(np.array([86.0, 87.0, 88.0, 89.0, 89.5, 89.9])) Rs = np.array([steel_rgb(np.array([math.cos(a)]))[0].mean() for a in angs]) mono = bool(np.all(np.diff(Rs) > 0)) ok &= mono print(" %-52s %10s want %8s %s" % ("rises monotonically from 86 to 89.9 deg", "-", "yes", "ok" if mono else "FAIL")) print() print("SELF-TEST %s" % ("PASSED" if ok else "FAILED")) return ok if __name__ == "__main__": ap = argparse.ArgumentParser() ap.add_argument("--self-test", action="store_true") args = ap.parse_args() if args.self_test: raise SystemExit(0 if self_test() else 1) print(__doc__)