#!/usr/bin/env catnip
# Réflexions, miroirs et rotors (algèbre géométrique)
#
# En algèbre géométrique 2D (Clifford Cl(2,0)), réfléchir un vecteur v par le
# miroir de normale unitaire n s'écrit avec le produit sandwich v' = -n v n.
# Le produit géométrique n v se décompose en n·v (scalaire) + n∧v (bivecteur, le
# coefficient de e12 = e1e2, aire orientée engendrée par n et v). Ce sandwich est
# identique à la formule vectorielle classique v' = v - 2(v·n) n.
#
# Composer DEUX réflexions donne une rotation : n2 n1 v n1 n2 = R v R̃, où le
# rotor R = n2 n1 est le produit géométrique des deux normales. Si n1 et n2 sont
# séparées d'un angle α, la rotation vaut 2α — le résultat « produit de deux
# réflexions = rotation » de l'algèbre géométrique.
#
# La sortie est une image PNG (un rayon rebondissant dans un heptagone de miroirs)
# et deux erreurs numériques max : sandwich vs formule vectorielle, et double
# réflexion vs rotation directe.
#
# DEPS: numpy, pillow
numpy = import('numpy')
math = import('math')
Image = import('PIL.Image')
ImageDraw = import('PIL.ImageDraw')
import('pathlib', 'Path')
script_dir = Path(META.file).parent
output_dir = script_dir / 'output'
output_dir.mkdir(exist_ok=True)
WIDTH = 720
HEIGHT = 720
CENTER_X = 360.0
CENTER_Y = 360.0
struct Vec2 {
x: float; y: float;
dot(self, o: Vec2): float => { self.x * o.x + self.y * o.y }
norm(self): float => { math.sqrt(self.dot(self)) }
sub(self, o: Vec2): Vec2 => { Vec2(self.x - o.x, self.y - o.y) }
normalized(self): Vec2 => { s = 1.0 / self.norm(); Vec2(self.x * s, self.y * s) }
}
# Réflexion par le produit sandwich v' = -n v n (n unitaire). On calcule le
# produit géométrique n v (partie scalaire n·v + partie bivecteur n∧v), on le
# multiplie par n à droite, puis on nie : le résultat redevient un vecteur.
reflect_ga = (v: Vec2, n: Vec2): Vec2 => {
scalar = n.x * v.x + n.y * v.y # n · v
bivector = n.x * v.y - n.y * v.x # n ∧ v (composante e12)
Vec2(
-(scalar * n.x + bivector * n.y),
-(scalar * n.y - bivector * n.x),
)
}
# Réflexion spéculaire vectorielle classique : on retranche deux fois la
# projection sur la normale. Sert de chemin de contrôle indépendant du sandwich.
reflect_vec = (v: Vec2, n: Vec2): Vec2 => {
proj = 2.0 * (v.x * n.x + v.y * n.y)
Vec2(v.x - proj * n.x, v.y - proj * n.y)
}
# Rotation directe d'angle a, référence pour la composition de deux réflexions.
rotate = (v: Vec2, a: float): Vec2 => {
Vec2(math.cos(a) * v.x - math.sin(a) * v.y, math.sin(a) * v.x + math.cos(a) * v.y)
}
# --- Contrôle 1 : sandwich GA vs formule vectorielle ---
# Un échantillon dense de normales (36 angles) croisé avec quelques vecteurs.
# Le broadcast sur une liste de structs passe chaque cas entier, sans risque de
# descente scalaire (contrairement à un ndarray).
struct Sample { v: Vec2; n: Vec2 }
test_vectors = list(Vec2(1.0, 0.0), Vec2(0.7, -1.3), Vec2(-2.1, 0.4), Vec2(0.0, 1.0))
samples = list()
a = 0
while a < 36 {
theta = a * math.pi / 18.0
n = Vec2(math.cos(theta), math.sin(theta))
for tv in test_vectors {
samples.append(Sample(tv, n))
}
a = a + 1
}
reflection_errors = samples.[(s) => { reflect_ga(s.v, s.n).sub(reflect_vec(s.v, s.n)).norm() }]
max_reflection_error = float(numpy.max(numpy.array(reflection_errors)))
print(f"⇒ Contrôle réflexion : sandwich -n v n vs v - 2(v·n)n")
print(f" {len(samples)} cas, erreur max = {max_reflection_error}")
# --- Contrôle 2 : deux réflexions = rotation d'angle 2α ---
# Normales n1 = e1 et n2 à l'angle α ; la double réflexion doit égaler rotate(2α).
alphas = list()
b = 1
while b < 24 {
alphas.append(b * math.pi / 24.0)
b = b + 1
}
probe = Vec2(1.3, -0.7)
rotation_errors = alphas.[(alpha) => {
n1 = Vec2(1.0, 0.0)
n2 = Vec2(math.cos(alpha), math.sin(alpha))
doubled = reflect_ga(reflect_ga(probe, n1), n2) # v'' = R v R̃, R = n2 n1
doubled.sub(rotate(probe, 2.0 * alpha)).norm()
}]
max_rotation_error = float(numpy.max(numpy.array(rotation_errors)))
print()
print(f"⇒ Contrôle rotor : deux réflexions (α) vs rotation directe (2α)")
print(f" {len(alphas)} angles, erreur max = {max_rotation_error}")
# --- Trajectoire : un rayon rebondit dans un heptagone de miroirs ---
# Chaque arête du polygone est un miroir ; la vitesse du rayon est réfléchie par
# la normale de l'arête au moyen du sandwich GA reflect_ga (mêmes maths que ci-
# dessus, appliquées à la direction du rayon).
struct Ray { pos: Vec2; dir: Vec2 }
struct Hit { point: Vec2; normal: Vec2 }
SIDES = 7
RADIUS = 320.0
mirror_poly = list()
k = 0
while k < SIDES {
phi = 2.0 * math.pi * k / SIDES - math.pi / 2.0
mirror_poly.append(Vec2(CENTER_X + RADIUS * math.cos(phi), CENTER_Y + RADIUS * math.sin(phi)))
k = k + 1
}
# Première arête rencontrée par le rayon (intersection segment/rayon la plus
# proche en avant). Le polygone est convexe et le rayon part de l'intérieur, donc
# il existe toujours exactement une arête sortante.
next_hit = (ray: Ray, poly): Hit => {
best_t = 1.0e18
best = Hit(ray.pos, Vec2(1.0, 0.0))
nverts = len(poly)
i = 0
while i < nverts {
va = poly[i]
vb = poly[(i + 1) % nverts]
ex = vb.x - va.x
ey = vb.y - va.y
denom = ray.dir.x * ey - ray.dir.y * ex # cross(dir, arête)
if denom != 0.0 {
apx = va.x - ray.pos.x
apy = va.y - ray.pos.y
t = (apx * ey - apy * ex) / denom # distance le long du rayon
u = (apx * ray.dir.y - apy * ray.dir.x) / denom # abscisse sur l'arête
if t > 1.0e-6 and u >= -1.0e-9 and u <= 1.0 + 1.0e-9 and t < best_t {
best_t = t
hit_pt = Vec2(ray.pos.x + t * ray.dir.x, ray.pos.y + t * ray.dir.y)
edge_normal = Vec2(-ey, ex).normalized() # perpendiculaire à l'arête
best = Hit(hit_pt, edge_normal)
}
}
i = i + 1
}
best
}
BOUNCES = 64
trace = (ray: Ray, poly, bounces: int) => {
points = list(ray.pos)
r = ray
i = 0
while i < bounces {
h = next_hit(r, poly)
points.append(h.point)
r = Ray(h.point, reflect_ga(r.dir, h.normal)) # réflexion GA de la direction
i = i + 1
}
points
}
start_angle = 0.371
start_ray = Ray(
Vec2(CENTER_X + 46.0, CENTER_Y + 8.0),
Vec2(math.cos(start_angle), math.sin(start_angle)),
)
path = trace(start_ray, mirror_poly, BOUNCES)
# Longueur totale parcourue, contrôle simple que le rayon rebondit bien.
total_length = 0.0
seg = 0
while seg < len(path) - 1 {
total_length = total_length + path[seg + 1].sub(path[seg]).norm()
seg = seg + 1
}
print()
print(f"⇒ Trajectoire : {BOUNCES} rebonds dans un heptagone de miroirs")
print(f" longueur totale = {round(total_length, 1)} px")
# --- Rendu PNG ---
image = Image.new('RGB', tuple(WIDTH, HEIGHT), tuple(16, 16, 22))
draw = ImageDraw.Draw(image)
# Miroirs : les arêtes du polygone.
m = 0
while m < SIDES {
va = mirror_poly[m]
vb = mirror_poly[(m + 1) % SIDES]
draw.line(list(tuple(va.x, va.y), tuple(vb.x, vb.y)), fill=tuple(90, 100, 130), width=3)
m = m + 1
}
# Trajectoire : un dégradé cyan → magenta suit l'ordre des rebonds.
segs = len(path) - 1
s = 0
while s < segs {
t01 = s / segs
red = int(60 + 195 * t01)
green = int(220 - 150 * t01)
blue = int(230 - 40 * t01)
p0 = path[s]
p1 = path[s + 1]
draw.line(list(tuple(p0.x, p0.y), tuple(p1.x, p1.y)), fill=tuple(red, green, blue), width=2)
s = s + 1
}
# Marqueur du point de départ.
sp = path[0]
draw.ellipse(list(tuple(sp.x - 5.0, sp.y - 5.0), tuple(sp.x + 5.0, sp.y + 5.0)), fill=tuple(250, 250, 120))
output_path = output_dir / 'ga_reflections.png'
image.save(str(output_path))
print()
print(f"⇒ Image → {output_path}")
# Aperçu navigateur : sert l'image une fois puis rend la main.
# --no-browser garde un chemin headless (le PNG reste écrit ci-dessus).
if '--no-browser' not in import('sys').argv {
http = import('http')
b64 = import('base64').b64encode(output_path.read_bytes()).decode('ascii')
http.serve(f'<!doctype html><meta charset="utf-8"><title>Geometric algebra reflections</title><body style="margin:0;background:#0d1117"><img style="max-width:100%;display:block;margin:0 auto" src="data:image/png;base64,{b64}">',
0, 'text/html; charset=utf-8', True)
}