geometryfactor.py

@@ -121,7 +121,6 @@ class plane:
         self.v = v
         self.w = w
         self.normalize()
-        self.cached_r = None
 
     def normalize(self):
         self.n = cross(self.v, self.w, axis=0).T
@@ -144,29 +143,15 @@ class plane:
     def inside(self, p):
         return (self.n @ (p-self.o))[0] >= -0.00001
 
-    def equivalent(self, other):
-        if length(cross(self.n, other.n)) > 1e-10:
-            return False
-        if (self.n.T @ other.n)[0,0] <= 0:
-            return False
-        if length(self.n @ (other.o - self.o)) > 1e-10:
-            return False
-        return True
-
     def intersection(self, r):
-        if self.cached_r is r and self.cached_v is r.v:
-            return self.solution
-        self.cached_r = r
-        self.cached_v = r.v
         self.M[:,2:3] = r.v
         try:
             c=solve(self.M, r.o-self.o)[...,2,:]
-            self.solution = (c, r.o - c*r.v)
         except Exception as e:
             if verbose:
                 stderr.write("solve: %s: %s\n" % (str(self.M), repr(e)))
-            self.solution = (-1, None)
+            return -1, None
-        return self.solution
+        return c, r.o - c*r.v
 
     def plane_intersection(self, other, limit=1e10, delta=1e-5):
         p = [ self.intersection(r)[1] for r in other.rays]
@@ -307,7 +292,7 @@ def prism(n=6, a=1.0, z=1.0, φ=0.0, n1=0, n2=None):
         plane(vector(z= z), vector(y=1), vector(x=1)),
     ]
     if n2 is None:
-        n2 = n-1
+        n2 = n
     else:
         c = cos(2*π/n * n1)
         s = sin(2*π/n * n1)
@@ -330,25 +315,6 @@ def prism(n=6, a=1.0, z=1.0, φ=0.0, n1=0, n2=None):
 
 µM = box(41., 41., 13.)
 
-def find_equivalent_planes(world):
-    done = []
-    for i,d in enumerate(world):
-        for j,p in enumerate(d):
-            if (i,j) in done:
-                continue
-            for ii,dd in enumerate(world):
-                if ii < i:
-                    continue
-                for jj,pp in enumerate(dd):
-                    if ii == i and jj <= j:
-                        continue
-                    if p.equivalent(pp):
-                        dd[jj] = p
-                        done.append((ii,jj))
-                        if do_size or verbose:
-                            print(f"equiv panes({len(done)}) {i}/{j} → {ii}/{jj}",
-                                  file=stderr)
-
 def run(world, r, Δr, Δθ, θ1=0, θ2=π/2, φ1=0, φ2=2*π):
     fmt = "%g %g %g %g "+len(world)*" %g"+"\n"
     nr = int(r/Δr + 1.001)
@@ -404,7 +370,7 @@ def run_beam(world, r, Δr, Δθ, θ1=0, θ2=π/2, φ1=0, φ2=2*π):
     n = int(N*(cos(θ1) - cos(θ2)) + 0.5)
     t = arange(n) + t1
     for θ in arccos(1 - t/N):
-        φ = (2*π * θ/Δθ) % (2*π)
+        φ = 2*π * θ/Δθ
         beam.bend(θ, φ)
         pl = [d.pathlength_beam(beam) for d in world]
         for i in range(beam.n):
@@ -440,7 +406,7 @@ short_opt = "hSNTREGU:D:C:B:P:µMr:d:a:p:x:"
 long_opt = ["help",
             "φ1=", "φ2=", "θ1=", "θ2=", "φ₁=", "φ₂=", "θ₁=", "θ₂=",
             "Δθ=", "µMustang", "test", "size", "run", "corners", "gnuplot",
-            "radius=", "resolution=", "Δr=", "equivalent",
+            "radius=", "resolution=", "Δr=" ,
             "detector", "copy=", "select=", "box=", "prism=",
             "plane=", "ray=",
             "rx=", "ry=", "rz=", "move=", "scale=",
@@ -480,7 +446,6 @@ do_test = False
 do_run = False
 do_corners = False
 do_gnuplot = False
-do_equiv = False
 rot = rotate(0,0)
 ori = vector()
 
@@ -543,8 +508,6 @@ for o,v in options:
         do_gnuplot=True
     if o in "-E --corners":
         do_corners=True
-    if o in "--equivalent":
-        do_equiv=True
     if o in "--rx":
         world[-1] = world[-1].transform(rotate(0, aval(v)))
     if o in "--ry":
@@ -585,17 +548,13 @@ for o,v in options:
             if k in vv[1]:
                 raise ValueError("recursive macro %s" % v)
             macros[k]=vv[1]
-            if verbose:
+            stderr.write("define macro %s=%s\n" % (k,vv[1]))
-                stderr.write("define macro %s=%s\n" % (k,vv[1]))
 
 if not world[-1]:
     world[-1:]=[]
 if not world:
     world=[µM]
 
-if do_equiv:
-    find_equivalent_planes(world)
-
 if not r or do_corners:
     max_r, c = corners(world)
     if not r: