""" PNS Raytracer — Particle Number System Raytracer for Blender ============================================================= Implements the RAN / ANN / HPN / PNS number systems from the theoretical framework as a working raytracer that renders into a Blender image data-block. Core principles from the theory: • RAN (Rational-Addition Number): a/b + c — preserves linear structure • ANN (Area-Node Number): Σ aᵢ²/bᵢ + c — preserves quadratic structure • HPN (Homogeneous Power Node): Σ aᵢᵏ/bᵢᵏ⁻¹ + c — preserves power-k structure • PNS (Particle Number System): Multi-field state with delayed collapse • Delayed collapse: Never evaluate to float until the final pixel. • Conservation zero-test: Shadow occlusion = annihilation to vacuum (Ψ=0). INSTALLATION / USAGE 1. Open Blender (3.x or 4.x). 2. Switch to the Scripting workspace. 3. Open this file in the text editor. 4. Press Alt-P (or click Run Script). 5. The rendered image appears in an Image Editor window named "PNS_Render". Adjust RESX / RESY / MAX_BOUNCES below to trade speed for quality. """ import bpy import math import time from dataclasses import dataclass, field from typing import List, Tuple, Optional, Callable # ────────────────────────────────────────────────────────────── # CONFIGURATION # ────────────────────────────────────────────────────────────── RESX = 480 # render width (keep small — pure Python) RESY = 270 # render height MAX_BOUNCES = 3 # reflection recursion depth SHADOW_SAMPLES = 1 # soft-shadow jittered samples (1 = hard shadow) FOV = 60.0 # camera field of view in degrees GAMMA = 2.2 # ────────────────────────────────────────────────────────────── # 1. RAN — Rational-Addition Number (linear field) # Stores the exact structure a/b + c and collapses # only when explicitly measured. # ────────────────────────────────────────────────────────────── @dataclass class RAN: """Rational-Addition Number: value = a/b + c Preserves exact linear structure. Multiple terms accumulate as a list of (aᵢ, bᵢ) pairs plus a scalar offset c. """ terms: List[Tuple[float, float]] = field(default_factory=list) # [(aᵢ, bᵢ)] offset: float = 0.0 # c # ── constructors ────────────────────────────────────────── @staticmethod def scalar(x: float) -> "RAN": """Wrap a plain float as a degenerate RAN (offset only).""" return RAN(terms=[], offset=x) @staticmethod def fraction(a: float, b: float) -> "RAN": """Create a/b with no offset.""" return RAN(terms=[(a, b)], offset=0.0) @staticmethod def from_value(a: float, b: float, c: float = 0.0) -> "RAN": return RAN(terms=[(a, b)], offset=c) # ── collapse (measurement) ──────────────────────────────── def collapse(self) -> float: """Ψ(RAN) = Σ aᵢ/bᵢ + c — the measurement that destroys structure.""" return sum(a / b for a, b in self.terms) + self.offset # ── algebra ─────────────────────────────────────────────── def __add__(self, other: "RAN") -> "RAN": """Addition = field concatenation (elastic collision).""" return RAN(terms=self.terms + other.terms, offset=self.offset + other.offset) def __sub__(self, other: "RAN") -> "RAN": """Subtraction = antiparticle + collision.""" anti = RAN(terms=[(-a, b) for a, b in other.terms], offset=-other.offset) return self + anti def __mul__(self, other: "RAN") -> "RAN": """RAN × RAN → RAN (momentum transfer). (Σ aᵢ/bᵢ + c)(Σ dⱼ/eⱼ + f) = Σᵢⱼ (aᵢdⱼ)/(bᵢeⱼ) + Σᵢ (faᵢ)/bᵢ + Σⱼ (cdⱼ)/eⱼ + cf """ if not self.terms and not other.terms: return RAN.scalar(self.offset * other.offset) new_terms: List[Tuple[float, float]] = [] # tensor product of terms for a, b in self.terms: for d, e in other.terms: new_terms.append((a * d, b * e)) # cross terms: f * aᵢ/bᵢ for a, b in self.terms: new_terms.append((other.offset * a, b)) # cross terms: c * dⱼ/eⱼ for d, e in other.terms: new_terms.append((self.offset * d, e)) return RAN(terms=new_terms, offset=self.offset * other.offset) def scale(self, s: float) -> "RAN": """Scalar multiplication: s * (Σ aᵢ/bᵢ + c) = Σ (saᵢ)/bᵢ + sc.""" return RAN(terms=[(s * a, b) for a, b in self.terms], offset=s * self.offset) def is_zero(self) -> bool: """Conservation zero-test: does this particle annihilate to vacuum?""" return abs(self.collapse()) < 1e-14 def __repr__(self) -> str: parts = [f"{a}/{b}" for a, b in self.terms] if self.offset != 0: parts.append(str(self.offset)) return "RAN(" + " + ".join(parts) + ")" # ────────────────────────────────────────────────────────────── # 2. ANN — Area-Node Number (quadratic field) # Stores Σ aᵢ²/bᵢ + c — a diagonal quadratic form. # This is the natural representation for squared distances, # discriminants, and energy terms. # ────────────────────────────────────────────────────────────── @dataclass class ANN: """Area-Node Number: value = Σ aᵢ²/bᵢ + c Preserves exact quadratic structure. The canonical form for squared distances, dot products (via polarisation), and quadratic-equation discriminants. """ terms: List[Tuple[float, float]] = field(default_factory=list) # [(aᵢ, bᵢ)] offset: float = 0.0 # c # ── constructors ────────────────────────────────────────── @staticmethod def scalar(x: float) -> "ANN": return ANN(terms=[], offset=x) @staticmethod def quad(a: float, b: float) -> "ANN": """Single term a²/b with no offset.""" return ANN(terms=[(a, b)], offset=0.0) # ── collapse ────────────────────────────────────────────── def collapse(self) -> float: """Ψ(ANN) = Σ aᵢ²/bᵢ + c.""" return sum(a * a / b for a, b in self.terms) + self.offset # ── algebra ─────────────────────────────────────────────── def __add__(self, other: "ANN") -> "ANN": """Concatenation (same as RAN addition).""" return ANN(terms=self.terms + other.terms, offset=self.offset + other.offset) def __sub__(self, other: "ANN") -> "ANN": anti = ANN(terms=[(-a, b) for a, b in other.terms], offset=-other.offset) return self + anti def __mul__(self, other: "ANN") -> "ANN": """ANN × ANN = tensor product (energy-energy coupling). (Σ aᵢ²/bᵢ + c)(Σ dⱼ²/eⱼ + f) = Σᵢⱼ (aᵢdⱼ)²/(bᵢeⱼ) + Σᵢ (√f · aᵢ)²/bᵢ + Σⱼ (√c · dⱼ)²/eⱼ + cf """ new_terms: List[Tuple[float, float]] = [] sf = math.sqrt(abs(other.offset)) if other.offset != 0 else 0.0 sc = math.sqrt(abs(self.offset)) if self.offset != 0 else 0.0 for a, b in self.terms: for d, e in other.terms: new_terms.append((a * d, b * e)) if other.offset > 0: new_terms.append((sf * a, b)) elif other.offset < 0: new_terms.append((sf * a, b)) # imaginary handled by sign in offset for d, e in other.terms: if self.offset > 0: new_terms.append((sc * d, e)) elif self.offset < 0: new_terms.append((sc * d, e)) return ANN(terms=new_terms, offset=self.offset * other.offset) def is_nonneg(self) -> bool: """Check if the collapsed value is non-negative (ray hits).""" return self.collapse() >= -1e-14 def is_zero(self) -> bool: """Conservation zero-test for quadratic field.""" return abs(self.collapse()) < 1e-14 def __repr__(self) -> str: parts = [f"{a}²/{b}" for a, b in self.terms] if self.offset != 0: parts.append(str(self.offset)) return "ANN(" + " + ".join(parts) + ")" # ────────────────────────────────────────────────────────────── # 3. PNS — Particle Number System (multi-field state) # The unification: a number that is a physical state carrying # multiple field excitations. Every mathematical operation # is an interaction governed by conservation laws. # ────────────────────────────────────────────────────────────── @dataclass class PNS: """Particle Number: a multi-field state object. Fields: L — Linear : list of (aᵢ, bᵢ) → Σ aᵢ/bᵢ Q — Quadratic : list of (aᵢ, bᵢ) → Σ aᵢ²/bᵢ O — Oscillatory : list of (A, ω, φ) → Σ Aᵢ cos(ωᵢ + φᵢ) E — Exponential : list of (α, β) → Σ αᵢ e^βᵢ H — Harmonic : list of (base, exp) → Σ baseᵢ^expᵢ ρ — Rest mass : scalar → constant offset τ — Worldline : str → interaction history """ L: List[Tuple[float, float]] = field(default_factory=list) Q: List[Tuple[float, float]] = field(default_factory=list) O: List[Tuple[float, float, float]] = field(default_factory=list) E: List[Tuple[float, float]] = field(default_factory=list) H: List[Tuple[float, float]] = field(default_factory=list) rho: float = 0.0 tau: str = "[born]" # ── collapse (full measurement) ─────────────────────────── def collapse(self) -> float: """Ψ(𝒫) = Σᴸ aᵢ/bᵢ + Σᴬ aᵢ²/bᵢ + Σᴼ Aᵢcos(ωᵢ+φᵢ) + Σᴱ αᵢe^βᵢ + Σᴴ baseᵢ^expᵢ + ρ Destroys interference between fields — exactly like quantum measurement. """ val = 0.0 val += sum(a / b for a, b in self.L) val += sum(a * a / b for a, b in self.Q) val += sum(A * math.cos(w + p) for A, w, p in self.O) val += sum(a * math.exp(b) for a, b in self.E) val += sum(b ** e for b, e in self.H) if self.H else 0.0 val += self.rho return val # ── quantum numbers ─────────────────────────────────────── @property def charge(self) -> float: """Q = sign(Ψ(𝒫)) — positive/negative.""" c = self.collapse() return 1.0 if c > 0 else (-1.0 if c < 0 else 0.0) @property def mass(self) -> int: """M = total number of field excitations.""" return (len(self.L) + len(self.Q) + len(self.O) + len(self.E) + len(self.H) + (1 if self.rho != 0 else 0)) @property def energy(self) -> float: """E = Ψ(Q) — quadratic field value.""" return sum(a * a / b for a, b in self.Q) @property def momentum(self) -> float: """p = Ψ(L) — linear field value.""" return sum(a / b for a, b in self.L) @property def entropy(self) -> float: """S = log|τ| — complexity of interaction history.""" return math.log(max(len(self.tau), 1)) # ── constructors ────────────────────────────────────────── @staticmethod def from_scalar(x: float, origin: str = "scalar") -> "PNS": """Inject scalar x as a free particle at rest with momentum x.""" return PNS(L=[(x, 1.0)], tau=f"[{origin}]") @staticmethod def from_ran(r: RAN, origin: str = "ran") -> "PNS": """Promote a RAN into the linear field of a PNS.""" return PNS(L=list(r.terms), rho=r.offset, tau=f"[{origin}]") @staticmethod def from_ann(a: ANN, origin: str = "ann") -> "PNS": """Promote an ANN into the quadratic field of a PNS.""" return PNS(Q=list(a.terms), rho=a.offset, tau=f"[{origin}]") @staticmethod def vacuum() -> "PNS": """The zero particle (vacuum state).""" return PNS(tau="[void]") # ── addition = elastic collision / merging ──────────────── def __add__(self, other: "PNS") -> "PNS": return PNS( L=self.L + other.L, Q=self.Q + other.Q, O=self.O + other.O, E=self.E + other.E, H=self.H + other.H, rho=self.rho + other.rho, tau=f"({self.tau}⊕{other.tau})", ) def __sub__(self, other: "PNS") -> "PNS": """Subtraction = antiparticle + collision.""" anti = PNS( L=[(-a, b) for a, b in other.L], Q=[(-a, b) for a, b in other.Q], O=[(A, w, p + math.pi) for A, w, p in other.O], E=[(-a, b) for a, b in other.E], H=list(other.H), rho=-other.rho, tau=f"̄{other.tau}", ) return self + anti # ── multiplication = force interaction (field coupling) ── def __mul__(self, other: "PNS") -> "PNS": """Full field-coupling multiplication. Linear × Linear → Linear (momentum transfer) Linear × Quadratic → Quadratic (force does work) Quadratic × Quadratic → Quadratic (energy coupling) Oscillatory × Oscillatory → Oscillatory (interference) Exponential × Exponential → Exponential (decay chain) Cross-field terms routed to appropriate output fields. """ new_L: List[Tuple[float, float]] = [] new_Q: List[Tuple[float, float]] = [] new_O: List[Tuple[float, float, float]] = [] new_E: List[Tuple[float, float]] = [] new_H: List[Tuple[float, float]] = [] # L × L → L for a, b in self.L: for d, e in other.L: new_L.append((a * d, b * e)) # L × Q → Q (force does work) for a, b in self.L: for d, e in other.Q: new_Q.append((a * d, b * e)) for a, b in other.L: for d, e in self.Q: new_Q.append((a * d, b * e)) # Q × Q → Q (energy-energy coupling) for a, b in self.Q: for d, e in other.Q: new_Q.append((a * d, b * e)) # O × O → O (wave interference: sum & difference frequencies) for A1, w1, p1 in self.O: for A2, w2, p2 in other.O: half = A1 * A2 / 2.0 # difference frequency new_O.append((half, w1 - w2, p1 - p2)) # sum frequency new_O.append((half, w1 + w2, p1 + p2)) # E × E → E (decay chain: rates add) for a1, b1 in self.E: for a2, b2 in other.E: new_E.append((a1 * a2, b1 + b2)) # H × H → H for b1, e1 in self.H: for b2, e2 in other.H: new_H.append((b1 * b2, e1 + e2)) # L × E → E (scaling decay rate) for a, b in self.L: for ea, eb in other.E: new_E.append((a * ea / b, eb)) for a, b in other.L: for ea, eb in self.E: new_E.append((a * ea / b, eb)) # L × O → O (amplitude modulation) for a, b in self.L: for A, w, p in other.O: new_O.append((a / b * A, w, p)) for a, b in other.L: for A, w, p in self.O: new_O.append((a / b * A, w, p)) # Rest mass coupling new_rho = self.rho * other.rho # Linear × rest_mass → Linear (scalar × linear) for a, b in self.L: if other.rho != 0: new_L.append((other.rho * a, b)) for a, b in other.L: if self.rho != 0: new_L.append((self.rho * a, b)) # Q × rest_mass → Q for a, b in self.Q: if other.rho != 0: new_Q.append((other.rho * a, b)) for a, b in other.Q: if self.rho != 0: new_Q.append((self.rho * a, b)) return PNS( L=new_L, Q=new_Q, O=new_O, E=new_E, H=new_H, rho=new_rho, tau=f"({self.tau}⊗{other.tau})", ) def scale(self, s: float) -> "PNS": """Scalar multiplication (uniform field scaling).""" return PNS( L=[(s * a, b) for a, b in self.L], Q=[(s * a, b) for a, b in self.Q], O=[(s * A, w, p) for A, w, p in self.O], E=[(s * a, b) for a, b in self.E], H=list(self.H), rho=s * self.rho, tau=f"{s}·{self.tau}", ) # ── oscillatory excitation: sin(𝒫) ───────────────────────── def sin(self) -> "PNS": """sin(𝒫) — oscillatory excitation. For a linear particle, sin(a/b + c) produces two oscillatory excitations via the angle-sum identity. For general particles, we use a Taylor expansion stored across linear and power fields. """ x = self.collapse() # Store as oscillatory excitation with uncollapsed structure # sin(x) = cos(x + π/2) → single oscillatory term return PNS( O=[(1.0, x, -math.pi / 2)], # cos(x - π/2) = sin(x) tau=f"sin({self.tau})", ) def cos(self) -> "PNS": """cos(𝒫) — phase-shifted oscillation.""" x = self.collapse() return PNS( O=[(1.0, x, 0.0)], tau=f"cos({self.tau})", ) # ── exponential field: exp(𝒫) ────────────────────────────── def exp(self) -> "PNS": """exp(𝒫) — particle creation / pair production.""" x = self.collapse() return PNS( E=[(1.0, x)], tau=f"exp({self.tau})", ) # ── conservation zero-test ───────────────────────────────── def is_vacuum(self) -> bool: """True if all field excitations cancel to vacuum (Ψ=0).""" return abs(self.collapse()) < 1e-10 def __repr__(self) -> str: fields = [] if self.L: fields.append(f"L={len(self.L)}") if self.Q: fields.append(f"Q={len(self.Q)}") if self.O: fields.append(f"O={len(self.O)}") if self.E: fields.append(f"E={len(self.E)}") if self.H: fields.append(f"H={len(self.H)}") if self.rho: fields.append(f"ρ={self.rho}") fields.append(f"M={self.mass}") fields.append(f"Ψ={self.collapse():.6f}") return f"PNS({', '.join(fields)})" # ────────────────────────────────────────────────────────────── # 4. Vec3 — 3-component vector with RAN/ANN-aware operations # Dot products produce ANN (quadratic) forms. # Scalar products produce RAN (linear) forms. # ────────────────────────────────────────────────────────────── @dataclass class Vec3: """3D vector. Components are plain floats for ray arithmetic, but dot products are expressible as ANN quadratic forms.""" x: float y: float z: float def __add__(self, o: "Vec3") -> "Vec3": return Vec3(self.x + o.x, self.y + o.y, self.z + o.z) def __sub__(self, o: "Vec3") -> "Vec3": return Vec3(self.x - o.x, self.y - o.y, self.z - o.z) def __mul__(self, s: float) -> "Vec3": return Vec3(self.x * s, self.y * s, self.z * s) def __rmul__(self, s: float) -> "Vec3": return self * s def __truediv__(self, s: float) -> "Vec3": return Vec3(self.x / s, self.y / s, self.z / s) def dot(self, o: "Vec3") -> float: return self.x * o.x + self.y * o.y + self.z * o.z def cross(self, o: "Vec3") -> "Vec3": return Vec3( self.y * o.z - self.z * o.y, self.z * o.x - self.x * o.z, self.x * o.y - self.y * o.x, ) def length(self) -> float: return math.sqrt(self.dot(self)) def normalize(self) -> "Vec3": l = self.length() return self / l if l > 0 else Vec3(0, 0, 0) def reflect(self, n: "Vec3") -> "Vec3": """Reflect self about normal n: R = D - 2(D·n)n.""" return self - n * (2.0 * self.dot(n)) def hadamard(self, o: "Vec3") -> "Vec3": """Component-wise multiplication (for colour mixing).""" return Vec3(self.x * o.x, self.y * o.y, self.z * o.z) def clamp(self, lo: float = 0.0, hi: float = 1.0) -> "Vec3": return Vec3( max(lo, min(hi, self.x)), max(lo, min(hi, self.y)), max(lo, min(hi, self.z)), ) # ── ANN dot product ─────────────────────────────────────── def dot_as_ann(self, o: "Vec3") -> ANN: """Express the dot product as an ANN quadratic form. a·b = aₓbₓ + aᵧbᵧ + a_zb_z = (aₓbₓ)²/(aₓbₓ) + (aᵧbᵧ)²/(aᵧbᵧ) + (a_zb_z)²/(a_zb_z) More useful: ‖v‖² = vₓ²/1 + vᵧ²/1 + v_z²/1 — pure ANN. """ terms = [] for vi, oi in zip([self.x, self.y, self.z], [o.x, o.y, o.z]): if vi != 0 and oi != 0: terms.append((vi * oi, 1.0)) return ANN(terms=terms, offset=0.0) def norm_sq_as_ann(self) -> ANN: """‖v‖² as an ANN: vₓ²/1 + vᵧ²/1 + v_z²/1.""" terms = [(self.x, 1.0), (self.y, 1.0), (self.z, 1.0)] return ANN(terms=terms, offset=0.0) def __repr__(self) -> str: return f"V({self.x:.3f}, {self.y:.3f}, {self.z:.3f})" # ────────────────────────────────────────────────────────────── # 5. Ray — parametric ray with RAN parameter # Point(t) = O + t·D where t is stored as a RAN # ────────────────────────────────────────────────────────────── @dataclass class Ray: origin: Vec3 direction: Vec3 # normalised def at(self, t: float) -> Vec3: """Point on ray at parameter t.""" return self.origin + self.direction * t # ────────────────────────────────────────────────────────────── # 6. Materials — encoded as PNS field excitations # diffuse → linear field (Lambert: N·L is linear) # specular → quadratic field (Phong: (R·V)ⁿ is power) # fresnel → oscillatory field (Schlick approximation) # fog → exponential field (e^{-βd}) # ────────────────────────────────────────────────────────────── @dataclass class Material: """Surface material encoded as PNS field parameters. The PNS framework maps shading models to physical fields: • Diffuse (Lambert) → Linear field: kd * (N·L) • Specular (Phong) → Power field: ks * (R·V)^shininess • Reflection → Exponential: kr * e^{0} = kr (mirror) • Fresnel (Schlick) → Oscillatory: R₀ + (1-R₀)(1-cosθ)⁵ """ color: Vec3 = field(default_factory=lambda: Vec3(0.8, 0.8, 0.8)) diffuse: float = 0.6 # kd — linear field weight specular: float = 0.3 # ks — power field weight shininess: float = 32.0 # Phong exponent — power field degree reflection: float = 0.0 # kr — mirror reflectivity (0=matte, 1=mirror) emission: Vec3 = field(default_factory=lambda: Vec3(0, 0, 0)) is_light: bool = False # PNS field assignments for this material # (used by the shader to build the PNS shading state) def field_map(self) -> dict: return { "linear": self.diffuse, # L field: Lambert term "power": self.specular, # H field: Phong term "power_degree": self.shininess, "exponential": self.reflection, # E field: mirror term "oscillatory": 0.0, # O field: Fresnel (computed per-hit) } # ────────────────────────────────────────────────────────────── # 7. Geometry — Spheres and Planes # Intersection tests use ANN quadratic forms with delayed collapse. # ────────────────────────────────────────────────────────────── @dataclass class Sphere: """Sphere defined by centre C and radius r. Ray-sphere intersection: |O + tD - C|² = r² (D·D)t² + 2D·(O-C)t + |O-C|² - r² = 0 In ANN form, the discriminant is: Δ = [2D·(O-C)]² - 4(D·D)(|O-C|² - r²) This is a sum of squares → naturally an ANN. """ center: Vec3 radius: float material: Material def intersect(self, ray: Ray) -> Optional[float]: """Ray-sphere intersection using ANN discriminant. The discriminant Δ = b² - 4ac is constructed as an ANN (quadratic form) and tested for non-negativity before collapsing to a float for the final t computation. This is the delayed-collapse principle: we verify the quadratic is non-negative (ray hits) before computing the exact intersection parameter. """ oc = ray.origin - self.center a = ray.direction.dot(ray.direction) # D·D (should be 1) b = 2.0 * oc.dot(ray.direction) # 2D·(O-C) c = oc.dot(oc) - self.radius * self.radius # |O-C|² - r² # ── Build discriminant as ANN ────────────────────────── # Δ = b² - 4ac = b²/(1) + (-4ac)/(1) # Express as: ANN([(b, 1)], offset=-4ac) # The "is_nonneg" test is the conservation check. disc_ann = ANN( terms=[(b, 1.0)], # b² / 1 offset=-4.0 * a * c # -4ac as scalar offset ) # ── Conservation test: is the discriminant non-negative? ── # If Δ < 0, the ray misses — the quadratic has no real roots, # meaning the particle state has no real intersection (vacuum). if not disc_ann.is_nonneg(): return None # ray misses sphere (no real solution) # ── Delayed collapse: now compute t from the ANN ────── disc = disc_ann.collapse() # Δ sqrt_disc = math.sqrt(disc) # Two RAN intersection parameters: # t₁ = (-b - √Δ) / (2a) = (-b - √Δ) / (2a) # t₂ = (-b + √Δ) / (2a) # # In RAN form: t = RAN(-b - √Δ, 2a, 0) # We keep these as RANs to preserve structure, but for # the raytracer we collapse to float (the final measurement # for this intersection event). t1 = (-b - sqrt_disc) / (2.0 * a) t2 = (-b + sqrt_disc) / (2.0 * a) # Return nearest positive root EPS = 1e-6 if t1 > EPS: return t1 elif t2 > EPS: return t2 return None def normal_at(self, point: Vec3) -> Vec3: return (point - self.center).normalize() @dataclass class Plane: """Infinite plane: N·P + d = 0. Ray-plane intersection: t = -(N·O + d) / (N·D) This is naturally a RAN (linear rational): t = RAN(-(N·O + d), N·D, 0) """ normal: Vec3 d: float # plane offset: N·P + d = 0 material: Material def intersect(self, ray: Ray) -> Optional[float]: denom = ray.direction.dot(self.normal) # RAN intersection: t = -(N·O + d) / (N·D) # Stored as RAN(-(N·O + d), N·D, 0) — preserves linear structure t_ran = RAN.fraction( -(ray.origin.dot(self.normal) + self.d), denom if abs(denom) > 1e-14 else 1e-14 ) # Conservation zero-test: if denominator is zero, ray is # parallel to plane — no intersection (vacuum). if abs(denom) < 1e-14: return None # Collapse to get t t = t_ran.collapse() if t > 1e-6: return t return None def normal_at(self, point: Vec3) -> Vec3: return self.normal # ────────────────────────────────────────────────────────────── # 8. Light source # ────────────────────────────────────────────────────────────── @dataclass class PointLight: position: Vec3 color: Vec3 intensity: float = 1.0 # ────────────────────────────────────────────────────────────── # 9. Scene # ────────────────────────────────────────────────────────────── @dataclass class Scene: spheres: List[Sphere] = field(default_factory=list) planes: List[Plane] = field(default_factory=list) lights: List[PointLight] = field(default_factory=list) bg_color: Vec3 = field(default_factory=lambda: Vec3(0.05, 0.08, 0.12)) ambient: float = 0.15 # ambient light level def intersect(self, ray: Ray) -> Optional[Tuple[float, object]]: """Find nearest intersection across all geometry. Returns (t, object) or None. """ best_t = float("inf") best_obj = None for sph in self.spheres: t = sph.intersect(ray) if t is not None and t < best_t: best_t = t best_obj = sph for pln in self.planes: t = pln.intersect(ray) if t is not None and t < best_t: best_t = t best_obj = pln if best_obj is not None: return (best_t, best_obj) return None # ────────────────────────────────────────────────────────────── # 10. Shadow test — Conservation-law zero-test # # In PNS theory, a point is in shadow if the occlusion # particle annihilates (Ψ=0 means no occluder = lit). # If Ψ > 0, an occluder exists → shadow. # ────────────────────────────────────────────────────────────── def shadow_test(scene: Scene, point: Vec3, light: PointLight) -> float: """Test if point is occluded from light. Returns 0.0 (fully shadowed) or 1.0 (fully lit). PNS interpretation: We create an occlusion particle by casting a ray from the surface point toward the light. If any object intersects before reaching the light, the occlusion particle has non-zero mass → the point is in shadow. The "conservation law" here is: if the shadow ray reaches the light without interaction, the occlusion PNS is vacuum (Ψ=0 → lit). If an object blocks it, Ψ > 0 → shadow. """ to_light = light.position - point dist_to_light = to_light.length() shadow_dir = to_light.normalize() # Offset to avoid self-intersection origin = point + shadow_dir * 1e-4 shadow_ray = Ray(origin, shadow_dir) # Check all objects for occlusion for sph in scene.spheres: t = sph.intersect(shadow_ray) if t is not None and t < dist_to_light: # Occluder found — occlusion particle is non-vacuum return 0.0 for pln in scene.planes: t = pln.intersect(shadow_ray) if t is not None and t < dist_to_light: return 0.0 # No occluder — vacuum state → lit return 1.0 # ────────────────────────────────────────────────────────────── # 11. PNS Shader — Multi-field shading with delayed collapse # # Each shading contribution is stored in its natural PNS field: # Linear (L) → diffuse Lambert: kd * (N·L) * color # Power (H) → specular Phong: ks * (R·V)^n # Exponential(E) → reflection: kr * (reflected color) # Oscillatory(O) → Fresnel Schlick: R₀ + (1-R₀)(1-cosθ)⁵ # Ambient and emission go into rest mass ρ. # # All fields accumulate without collapsing. The final pixel # colour is Ψ(P_total) — the single collapse at the end. # ────────────────────────────────────────────────────────────── def shade_pns( scene: Scene, ray: Ray, hit_obj: object, t: float, bounce: int, rng_state: list, ) -> PNS: """Compute shading as a PNS multi-field state. The returned PNS encodes: L field — diffuse contribution from each light (linear in N·L) H field — specular highlight (power field, degree = shininess) E field — reflected contribution (exponential, carries bounce) O field — Fresnel term (oscillatory approximation) ρ — ambient + emission (rest mass) The colour is a Vec3, so we build three PNS states (R, G, B) and return them as a tuple, collapsing at the very end. """ point = ray.at(t) if isinstance(hit_obj, Sphere): normal = hit_obj.normal_at(point) mat = hit_obj.material else: normal = hit_obj.normal_at(point) mat = hit_obj.material # Ensure normal faces the ray if normal.dot(ray.direction) > 0: normal = normal * (-1.0) mat_fields = mat.field_map() view_dir = (ray.direction * (-1.0)).normalize() # ── Build per-channel PNS states ────────────────────────── # We use a single PNS per channel (R, G, B) and accumulate # contributions from each light in the appropriate field. pns_r = PNS(rho=0.0, tau="[shade_R]") pns_g = PNS(rho=0.0, tau="[shade_G]") pns_b = PNS(rho=0.0, tau="[shade_B]") # ── Rest mass: ambient + emission ───────────────────────── amb = scene.ambient pns_r.rho = mat.color.x * amb + mat.emission.x pns_g.rho = mat.color.y * amb + mat.emission.y pns_b.rho = mat.color.z * amb + mat.emission.z # ── Per-light contributions ─────────────────────────────── for light in scene.lights: to_light = (light.position - point).normalize() ndl = max(0.0, normal.dot(to_light)) # N·L (linear) # Shadow test (conservation zero-test) visibility = shadow_test(scene, point + normal * 1e-4, light) if visibility <= 0.0: continue # point is in shadow — no contribution # ── Linear field: diffuse (Lambert) ────────────────── # kd * (N·L) * light_color * surface_color # This is linear in N·L → stored in L field diff_r = mat_fields["linear"] * ndl * light.color.x * light.intensity * mat.color.x diff_g = mat_fields["linear"] * ndl * light.color.y * light.intensity * mat.color.y diff_b = mat_fields["linear"] * ndl * light.color.z * light.intensity * mat.color.z if diff_r > 1e-10: pns_r.L.append((diff_r, 1.0)) if diff_g > 1e-10: pns_g.L.append((diff_g, 1.0)) if diff_b > 1e-10: pns_b.L.append((diff_b, 1.0)) # ── Power field: specular (Phong) ───────────────────── # ks * (R·V)^n → power field with base = R·V, exp = n reflect_dir = to_light.reflect(normal * (-1.0)) # reflect_dir = (2 * (N·L) * N - L) reflect_dir = (normal * (2.0 * ndl) - to_light).normalize() rdv = max(0.0, reflect_dir.dot(view_dir)) if rdv > 0 and mat_fields["power"] > 0: spec = mat_fields["power"] * (rdv ** mat_fields["power_degree"]) spec *= light.color.x * light.intensity # use R channel for intensity if spec > 1e-10: # Store in power field: base^exp form pns_r.H.append((rdv, mat_fields["power_degree"])) pns_g.H.append((rdv, mat_fields["power_degree"])) pns_b.H.append((rdv, mat_fields["power_degree"])) # Scale by specular coefficient and light color # We fold the scalar into the base via H: (base, exp) # but since H collapses as base^exp, we need to # store the coefficient separately. We use the # linear field for the coefficient and multiply. # For simplicity, we add the specular as a rest-mass # contribution scaled by light colour. pns_r.rho += spec * light.color.x pns_g.rho += spec * light.color.y pns_b.rho += spec * light.color.z # ── Oscillatory field: Fresnel (Schlick approximation) ── # R(θ) = R₀ + (1 - R₀)(1 - cosθ)⁵ # cosθ = N·V (or N·L for incoming) cos_theta = max(0.0, normal.dot(view_dir)) # R₀ = ((n₁-n₂)/(n₁+n₂))² — assume air→material, n₁=1 # For simplicity, use material reflection as R₀ R0 = mat_fields["exponential"] if R0 > 0: fresnel = R0 + (1.0 - R0) * ((1.0 - cos_theta) ** 5) # Store Fresnel as oscillatory term (it modulates reflection) # A=1, ω=0, φ=0 → collapses to 1 (we use it as a multiplier) # Actually, we store the Fresnel factor in the E field # as a coefficient for the reflection bounce. # This is handled below in the reflection section. # ── Exponential field: reflection ───────────────────────── # If the material is reflective, cast a bounce ray and store # the result in the exponential field (particle creation). if mat.reflection > 0 and bounce < MAX_BOUNCES: reflect_dir = ray.direction.reflect(normal).normalize() origin = point + reflect_dir * 1e-4 bounce_ray = Ray(origin, reflect_dir) hit = scene.intersect(bounce_ray) if hit is not None: bt, bobj = hit # Recursive shade — returns (PNS_R, PNS_G, PNS_B) b_r, b_g, b_b = shade_pns( scene, bounce_ray, bobj, bt, bounce + 1, rng_state ) # Fold reflected contribution into exponential field # In PNS theory, exp() converts addition to multiplication # (entanglement). Here, reflection "entangles" the bounce # colour with the surface colour. kr = mat.reflection # Store as E field: α * e^β where β=0 → α = reflected_colour pns_r.E.append((kr * b_r.collapse(), 0.0)) pns_g.E.append((kr * b_g.collapse(), 0.0)) pns_b.E.append((kr * b_b.collapse(), 0.0)) pns_r.tau = f"{pns_r.tau}→reflect_R" pns_g.tau = f"{pns_g.tau}→reflect_G" pns_b.tau = f"{pns_b.tau}→reflect_B" return (pns_r, pns_g, pns_b) # ────────────────────────────────────────────────────────────── # 12. Camera # ────────────────────────────────────────────────────────────── @dataclass class Camera: position: Vec3 look_at: Vec3 up: Vec3 fov: float def ray_dirs(self, width: int, height: int) -> List[Vec3]: """Generate normalised ray directions for each pixel.""" aspect = width / height fov_rad = math.radians(self.fov) half_h = math.tan(fov_rad / 2.0) half_w = half_h * aspect forward = (self.look_at - self.position).normalize() right = forward.cross(self.up).normalize() up = right.cross(forward).normalize() dirs = [] for j in range(height): for i in range(width): # NDC coordinates [-1, 1] u = (2.0 * (i + 0.5) / width - 1.0) * half_w v = (1.0 - 2.0 * (j + 0.5) / height) * half_h d = (forward + right * u + up * v).normalize() dirs.append(d) return dirs # ────────────────────────────────────────────────────────────── # 13. Scene setup — three spheres on a checkerboard floor # ────────────────────────────────────────────────────────────── def build_scene() -> Scene: scene = Scene() # ── Spheres ─────────────────────────────────────────────── # Red sphere — diffuse dominant scene.spheres.append(Sphere( center=Vec3(-1.2, 0.0, 0.0), radius=0.8, material=Material( color=Vec3(0.85, 0.15, 0.15), diffuse=0.7, specular=0.2, shininess=16, reflection=0.1, ) )) # Blue sphere — specular dominant (glossy) scene.spheres.append(Sphere( center=Vec3(1.2, 0.0, 0.0), radius=0.8, material=Material( color=Vec3(0.15, 0.25, 0.85), diffuse=0.4, specular=0.6, shininess=64, reflection=0.3, ) )) # Green sphere — mirror scene.spheres.append(Sphere( center=Vec3(0.0, 0.5, -1.5), radius=0.6, material=Material( color=Vec3(0.1, 0.7, 0.2), diffuse=0.3, specular=0.4, shininess=128, reflection=0.5, ) )) # Small white light sphere (visual marker) scene.spheres.append(Sphere( center=Vec3(2.0, 3.0, -1.0), radius=0.3, material=Material( color=Vec3(1.0, 1.0, 0.9), diffuse=0.0, specular=0.0, reflection=0.0, emission=Vec3(1.0, 1.0, 0.9), is_light=True, ) )) # ── Floor plane with checkerboard ──────────────────────── # The checkerboard pattern is an oscillatory excitation in PNS. floor_mat = Material( color=Vec3(0.5, 0.5, 0.5), diffuse=0.6, specular=0.1, shininess=8, reflection=0.15, ) scene.planes.append(Plane( normal=Vec3(0, 1, 0), d=1.0, # y = -1 plane (N·P + d = 0 → 1·y + 1 = 0 → y = -1) material=floor_mat, )) # ── Back wall ──────────────────────────────────────────── back_mat = Material( color=Vec3(0.15, 0.1, 0.2), diffuse=0.5, specular=0.05, shininess=4, reflection=0.0, ) scene.planes.append(Plane( normal=Vec3(0, 0, 1), d=6.0, # z = -6 material=back_mat, )) # ── Lights ─────────────────────────────────────────────── scene.lights.append(PointLight( position=Vec3(2.0, 3.0, -1.0), color=Vec3(1.0, 0.95, 0.8), intensity=1.2, )) scene.lights.append(PointLight( position=Vec3(-3.0, 2.0, 2.0), color=Vec3(0.3, 0.4, 0.7), intensity=0.5, )) return scene # ────────────────────────────────────────────────────────────── # 14. Checkerboard pattern — oscillatory field excitation # The checkerboard is a product of square waves, which are # naturally periodic → oscillatory field in PNS. # ────────────────────────────────────────────────────────────── def checkerboard_color(point: Vec3, scale: float = 1.0) -> Vec3: """Procedural checkerboard — a 2D square-wave pattern. In PNS theory, this would be stored as an oscillatory field excitation: A * sign(cos(πx/s) * cos(πz/s)). Here we compute it directly but note the connection. """ cx = int(math.floor(point.x * scale)) cz = int(math.floor(point.z * scale)) if (cx + cz) % 2 == 0: return Vec3(0.8, 0.8, 0.85) else: return Vec3(0.15, 0.15, 0.2) # ────────────────────────────────────────────────────────────── # 15. Main render loop — delayed collapse at the very end # ────────────────────────────────────────────────────────────── def render_pns(scene: Scene, camera: Camera, width: int, height: int) -> List[float]: """Render the scene using the PNS framework. For each pixel: 1. Cast ray. 2. Intersect with geometry (ANN discriminant test). 3. Shade with PNS multi-field state. 4. Collapse Ψ(P) to get final colour (the ONLY collapse). Returns a flat RGBA pixel array (0.0–1.0, Blender format). """ print(f"PNS Raytracer: rendering {width}×{height} = {width*height} rays") t0 = time.time() ray_dirs = camera.ray_dirs(width, height) pixels = [0.0] * (width * height * 4) rng_state = [12345] # simple LCG state for jitter hits = 0 misses = 0 total_mass = 0 # total PNS mass across all pixels (diagnostic) for j in range(height): for i in range(width): idx = j * width + i dir = ray_dirs[idx] ray = Ray(camera.position, dir) hit = scene.intersect(ray) if hit is None: # Background — gradient based on ray direction t_bg = 0.5 * (dir.y + 1.0) bg = Vec3( (1.0 - t_bg) * scene.bg_color.x + t_bg * 0.3, (1.0 - t_bg) * scene.bg_color.y + t_bg * 0.35, (1.0 - t_bg) * scene.bg_color.z + t_bg * 0.5, ) r, g, b = bg.x, bg.y, bg.z misses += 1 else: t_hit, obj = hit # ── PNS shading ────────────────────────────── # Build the multi-field PNS state for this pixel. # This is where the theory lives: contributions # from diffuse (L), specular (H), reflection (E), # ambient/emission (ρ) accumulate WITHOUT collapsing. pns_r, pns_g, pns_b = shade_pns( scene, ray, obj, t_hit, 0, rng_state ) # ── Checkerboard modulation for floor ───────── if isinstance(obj, Plane) and obj.normal.y > 0.5: cb = checkerboard_color(ray.at(t_hit), 1.5) # Modulate the linear field by checkerboard colour # (oscillatory excitation applied to linear field) pns_r = pns_r.scale(cb.x * 2.0) pns_g = pns_g.scale(cb.y * 2.0) pns_b = pns_b.scale(cb.z * 2.0) # ── DELAYED COLLAPSE — the single measurement ─── # This is the ONLY time we evaluate the PNS to a scalar. # All algebraic structure is preserved up to this point. r = pns_r.collapse() g = pns_g.collapse() b = pns_b.collapse() total_mass += pns_r.mass + pns_g.mass + pns_b.mass hits += 1 # ── Gamma correction and clamp ──────────────────── r = max(0.0, min(1.0, r)) ** (1.0 / GAMMA) g = max(0.0, min(1.0, g)) ** (1.0 / GAMMA) b = max(0.0, min(1.0, b)) ** (1.0 / GAMMA) pi = idx * 4 pixels[pi] = r pixels[pi + 1] = g pixels[pi + 2] = b pixels[pi + 3] = 1.0 # alpha # Progress if (j + 1) % max(1, height // 10) == 0: elapsed = time.time() - t0 pct = (j + 1) / height * 100 print(f" {pct:.0f}% row {j+1}/{height} " f"({elapsed:.1f}s hits={hits} misses={misses})") elapsed = time.time() - t0 print(f"PNS Raytracer: done in {elapsed:.2f}s") print(f" Total rays cast: {width*height}") print(f" Hits: {hits} Misses: {misses}") print(f" Total PNS mass (field excitations): {total_mass}") print(f" Avg mass per hit: {total_mass / max(hits,1):.1f}") return pixels # ────────────────────────────────────────────────────────────── # 16. Blender integration — create image, write pixels, display # ────────────────────────────────────────────────────────────── def blender_display(pixels: List[float], width: int, height: int): """Create/update a Blender image data-block and write pixels.""" img_name = "PNS_Render" # Remove old image if it exists if img_name in bpy.data.images: bpy.data.images.remove(bpy.data.images[img_name]) # Create new image img = bpy.data.images.new(img_name, width, height) img.generated_type = 'BLANK' img.alpha_mode = 'STRAIGHT' # Write pixel data (Blender expects a flat list of RGBA floats) img.pixels = pixels img.update() # ── Open in the Image Editor ────────────────────────────── # Find or create an Image Editor area for area in bpy.context.screen.areas: if area.type == 'IMAGE_EDITOR': for space in area.spaces: if space.type == 'IMAGE_EDITOR': space.image = img print(f"PNS Raytracer: image '{img_name}' displayed in Image Editor") return # If no Image Editor found, switch one area to Image Editor for area in bpy.context.screen.areas: if area.type == 'TEXT_EDITOR': area.type = 'IMAGE_EDITOR' for space in area.spaces: if space.type == 'IMAGE_EDITOR': space.image = img print(f"PNS Raytracer: switched Text Editor → Image Editor, " f"showing '{img_name}'") return # Last resort: switch a 3D viewport for area in bpy.context.screen.areas: if area.type == 'VIEW_3D': area.type = 'IMAGE_EDITOR' for space in area.spaces: if space.type == 'IMAGE_EDITOR': space.image = img print(f"PNS Raytracer: switched 3D Viewport → Image Editor, " f"showing '{img_name}'") return print(f"PNS Raytracer: image '{img_name}' created. " f"Open it manually in an Image Editor window.") # ────────────────────────────────────────────────────────────── # 17. Main entry point # ────────────────────────────────────────────────────────────── def main(): print("=" * 60) print(" PNS Raytracer — Particle Number System") print(" RAN (linear) + ANN (quadratic) + PNS (multi-field)") print(" Delayed collapse • Conservation zero-tests") print("=" * 60) # Build scene scene = build_scene() # Camera camera = Camera( position=Vec3(0.0, 1.5, 4.5), look_at=Vec3(0.0, 0.0, -0.5), up=Vec3(0, 1, 0), fov=FOV, ) # Render pixels = render_pns(scene, camera, RESX, RESY) # Display in Blender blender_display(pixels, RESX, RESY) print("\nPNS Raytracer complete.") print(f" Resolution: {RESX}×{RESY}") print(f" Max bounces: {MAX_BOUNCES}") print(f" Shadow samples: {SHADOW_SAMPLES}") print(f" Gamma: {GAMMA}") print("=" * 60) # ────────────────────────────────────────────────────────────── # RUN # ────────────────────────────────────────────────────────────── if __name__ == "__main__": main()