add --equivalent

find equivalent planes in all detektor volumes and replace duplicates

geometryfactor.py

@@ -121,6 +121,7 @@ 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
@@ -143,15 +144,29 @@ 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)))
-            return -1, None
+            self.solution = (-1, None)
-        return c, r.o - c*r.v
+        return self.solution
 
     def plane_intersection(self, other, limit=1e10, delta=1e-5):
         p = [ self.intersection(r)[1] for r in other.rays]
@@ -292,7 +307,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
+        n2 = n-1
     else:
         c = cos(2*π/n * n1)
         s = sin(2*π/n * n1)
@@ -315,6 +330,25 @@ 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)
@@ -370,7 +404,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):
@@ -406,7 +440,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=" ,
+            "radius=", "resolution=", "Δr=", "equivalent",
             "detector", "copy=", "select=", "box=", "prism=",
             "plane=", "ray=",
             "rx=", "ry=", "rz=", "move=", "scale=",
@@ -446,6 +480,7 @@ do_test = False
 do_run = False
 do_corners = False
 do_gnuplot = False
+do_equiv = False
 rot = rotate(0,0)
 ori = vector()
 
@@ -508,6 +543,8 @@ 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":
@@ -548,13 +585,17 @@ for o,v in options:
             if k in vv[1]:
                 raise ValueError("recursive macro %s" % v)
             macros[k]=vv[1]
-            stderr.write("define macro %s=%s\n" % (k,vv[1]))
+            if verbose:
+                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: