save for merge

Plane-class, change rays-class, add Scinitllator-class

BGO.py

@@ -0,0 +1,114 @@
+import json
+from Struktur import Plane, Ray, Rays
+import numpy as np
+
+#load the data from the json files
+with open('BGO_planes.json') as json_file:
+    data1 = json.load(json_file)
+planes_json=data1['Planes']
+with open('Detector_planes.json') as json_file:
+    data6 = json.load(json_file)
+planes_det_json=data6['Planes Det']
+
+
+def extract_planes(planes_json):
+    plane_list=[]
+    for key in planes_json:
+        plane_list.append(Plane(planes_json[key][0],planes_json[key][1],planes_json[key][2]))
+    return plane_list
+
+#extracts the example_points from the json file 
+#and returns them
+def extract_points(points_json):
+    points=[]
+    for key in points_json:
+        points.append(points_json[key])
+    return points
+
+
+def ray_way(ray,scinti_planes,scinti_shapes,point_list,ind):
+    hit_inc=False
+    num=0
+    new_ray=0
+    for i in range(len(scinti_shapes)):
+        cross_point=Plane.cross(scinti_planes[i],ray)
+        if cross_point!= None and i!=ind:
+            if Shape.hit(scinti_shapes[i],point_list[i],cross_point,len(scinti_shapes[i].g_list))==True:
+                p=np.dot(scinti_planes[i].norm,ray.dir)
+                if p<0:
+                    hit_inc=True
+                    num=i
+                    hit_point=cross_point
+    if hit_inc!=True:
+        ray.dir=-ray.dir
+        for i in range(len(scinti_shapes)):
+            cross_point=Plane.cross(scinti_planes[i],ray)
+            if cross_point!= None and i!=ind:
+                if Shape.hit(scinti_shapes[i],point_list[i],cross_point,len(scinti_shapes[i].g_list))==True:
+                    p=np.dot(scinti_planes[i].norm,ray.dir)
+                    if p<0:
+                        hit_inc=True
+                        num=i
+                        hit_point=cross_point
+    angle=Plane.cross_angle(scinti_planes[num],ray)
+    if angle<27.7177475:
+        new_dir=Plane.scatter(scinti_planes[num],ray)
+        new_ray=Ray(hit_point,hit_point+new_dir)
+        print('scatter')
+    else:
+        new_dir=Plane.reflect(scinti_planes[num],ray)
+        new_ray=Ray(hit_point,hit_point+new_dir)
+        print('reflect')
+    return [new_ray,num,hit_point]
+
+def destroy(det_shapes,det_planes,det_points,det_borders,ray):
+    res=False
+    for i in range(len(det_planes)):
+        cross_point=Plane.cross(det_planes[i],ray)
+        if cross_point!=None:
+            if Shape.hit(det_shapes[i],det_points[i],cross_point,len(det_borders[i]))==True:
+                res=i+1
+    if res!=False:
+        return res
+    else:
+        return ray
+    
+
+def ray_journey(ray,scinti_planes,scinti_shapes,scinti_points,det_shapes,det_planes,det_points,det_borders):
+    ray_list=[]
+    hit_ind=9
+    res=False
+    while res==False:
+        stop=destroy(det_shapes,det_planes,det_points,det_borders,ray)
+        if stop==1 or stop==2:
+            res=True
+            print('Hit Detector'+str(stop))
+        else:
+            if len(ray_list)<100:
+                obj=ray_way(ray,scinti_planes,scinti_shapes,scinti_points,hit_ind)
+                hit_ind=obj[1]
+                print(obj[0])
+                print(obj[1]+1)
+                ray_list.append(obj[0])
+                ray=obj[0]
+            else:
+                res=True
+    num=len(ray_list)
+    stretch=0
+    for i in range(num-1):
+        l=ray_list[i+1].origin-ray_list[i].origin
+        stretch=np.linalg.norm(l)*0.001
+    print(stretch)
+    return ray_list
+    
+def create_dat(filename, list):
+    with open(filename, 'w') as file:
+        for i in range(len(list)):
+            if isinstance(list[i],Ray):
+                a,b,c=list[i].origin
+                d,e,f=list[i].dir
+                file.write(f"{a}  {b}  {c}  {d}  {e}  {f}\n")
+    
+
+
+print(extract_planes(planes_json))

BGO_planes.json

@@ -0,0 +1,9 @@
+{"Planes":{ "Plane 1":[[0,0,0],[0,-25.981,45],[20,0,0]],
+            "Plane 2":[[0,0,0],[20,0,0],[0,51.962,0]],
+            "Plane 3":[[0,51.962,0],[20,51.962,0],[0,77.943,45]],
+            "Plane 4":[[0,51.962,90],[0,77.943,45],[20,51.962,90]],
+            "Plane 5":[[0,0,90],[0,51.962,90],[20,0,90]],
+            "Plane 6":[[0,0,90],[20,0,90],[0,-25.981,45]],
+            "Plane 7":[[0,0,0],[0,51.962,0],[0,0,90]],
+            "Plane 8":[[20,0,0],[20,0,90],[20,51.962,0]]}}
+

BGO_point_in_plane.json

@@ -0,0 +1,9 @@
+{"Example_point":{  "Plane 1":[10,-10,17.321],
+                    "Plane 2":[10,20,0],
+                    "Plane 3":[10,61.962,17.321],
+                    "Plane 4":[10,61.962,62.321],
+                    "Plane 5":[10,20,90],
+                    "Plane 6":[10,-10,62.321],
+                    "Plane 7":[0,20,45],
+                    "Plane 8":[20,20,45]}}
+

Detector_planes.json

@@ -0,0 +1,2 @@
+{"Planes Det":{ "Det 1":[[0,0,0],[20,0,0],[0,51.962,0]],
+                "Det 2":[[0,0,90],[20,0,90],[0,51.962,90]]}}

Detector_point_in_plane.json

@@ -0,0 +1,3 @@
+{"Example_point Det":{  "Det 1":[10,20,0],
+                        "Det 2":[10,20,90]}}
+                   

Random_gen.py

@@ -0,0 +1,41 @@
+import random
+from Struktur import Plane,Ray
+
+def test_planes(test_size):
+    rand_planes=[]
+    rand_vec=[]
+    vec=[]
+    for i in range(test_size*9):
+        vec.append(random.randrange(-100,100,1))
+        if len(vec)==3:
+            rand_vec.append(vec)
+            vec=[]
+        if len(rand_vec)==3:
+            rand_planes.append(Plane(rand_vec[0],rand_vec[1],rand_vec[2]))
+            rand_vec=[]
+    return rand_planes
+
+def test_lines(test_size):
+    rand_lines=[]
+    rand_vec=[]
+    vec=[]
+    for i in range(test_size*6):
+        vec.append(random.randrange(-100,100,1))
+        if len(vec)==3:
+            rand_vec.append(vec)
+            vec=[]
+        if len(rand_vec)==2:
+            rand_lines.append(Ray(rand_vec[0],rand_vec[1]))
+            rand_vec=[]
+    return rand_lines
+
+p=test_planes(1)
+g=test_lines(1)
+print(p[0].dir1)
+print(p[0].dir2)
+print(p[0].o)
+print(g[0].origin)
+print(g[0].dir)
+print(Plane.cross(p[0],g[0]))
+print(p[0].a)
+print(p[0].norm)

Struktur.py

@@ -0,0 +1,269 @@
+from sys import argv, stdout, stderr
+import numpy as np
+import random
+from numpy.linalg import solve
+
+#size of the simulation
+size=3
+
+def vtuple(v):
+    return tuple(v[:,0])
+
+def matrix(vec1,vec2,vec3):
+    return np.array([vec1,vec2,vec3])
+
+def norm_vec(vec):
+    return vec/np.linalg.norm(vec)
+
+def vector(x=0., y=0., z=0.):
+    return np.array([[x],[y],[z]])
+
+def length(v):
+    x=0
+    for i in range(len(v)):
+        x=x+v[i]**2
+    res=np.sqrt(x)
+    return res
+
+def scale(x=1., y=1., z=1.):
+    return np.array(((x,0,0),
+                  (0,y,0),
+                  (0,0,z)))
+
+def rotate(axis, φ):
+    m = scale()
+    c = np.cos(φ)
+    s = -np.sin(φ)
+    for i in range(3):
+        if i != axis:
+            m[i,i]=c
+            for j in range(3):
+                if j != i and j != axis:
+                    m[i,j] = s
+                    s = -s
+    return m
+
+#performs a rotation around any axis in the 3D space defined by a vector
+#using the Rodriguez formula
+def rotate_any(axis, vec, angle):
+    K=matrix([0,-axis[2],axis[1]],
+            [axis[2],0,axis[0]],
+            [-axis[1],axis[0],0])
+    A=np.identity(3) + np.sin(angle) * K + (1-np.cos(angle)) * np.dot(K,K)
+    return np.dot(A,vec)
+    
+class Ray:
+    #self.o is the starting point. self.d the direction vector of the line
+    #o and p are two points on the line
+    def __init__(self, o, p):
+        self.origin=o
+        self.dir=np.subtract(o,p) #calculates o-p
+
+    def __str__(self):
+        return str(self.origin)+str(self.dir)
+
+class rays:
+    "Fast reimplementation of a beam of rays"
+    
+    def __init__(self, r, Δr, origin):
+        self.origin=origin
+        nr = int(r/Δr + 1.001)
+        self.n = 0
+        rpsi = 0
+        for ir in range(nr):
+            nnpsi = 2*np.pi*ir + 0.5 + rpsi
+            npsi = int(nnpsi + 0.5)
+            rpsi = nnpsi - npsi
+            self.n += npsi
+        self.xy = np.empty((self.n,3,1))
+        self.r = np.empty((self.n,))
+        self.psi = np.empty((self.n,))
+        i=0
+        rpsi = 0
+        for ir in range(nr):
+            r = ir*Δr
+            nnpsi = 2*np.pi*ir + 0.5 + rpsi
+            npsi = int(nnpsi + 0.5)
+            rpsi = nnpsi - npsi
+            dpsi = 2*np.pi/npsi
+            for ipsi in range(npsi):
+                psi = ipsi * dpsi + dpsi/2
+                self.r[i] = r
+                self.psi[i] = psi
+                self.xy[i,:,0] = (r*np.cos(psi), r*np.sin(psi), 0)
+                i = i+1
+
+    def bend(self, th, phi):
+        self.th = th
+        self.phi = phi
+        self.v = rotate(2,phi) @ rotate(1,-th)
+        self.dir = self.v @ vector(z=1)
+        self.origin = (self.v @ self.xy)[...,0].T
+        #print(self.origin)
+    
+class Plane:
+
+    #o, v, w are points on the plane
+    def __init__(self, o, v, w):
+        self.o=o
+        self.v=v
+        self.w=w
+        self.dir1=o-v
+        self.dir2=v-w
+        self.norm=np.cross(self.dir1,self.dir2,axis=0).T
+        self.a=np.dot(self.norm,o)
+        self.tolerance=1e-9
+        self.A=np.empty((135,3,3))
+        self.A[...,1:2]=o-v
+        self.A[...,2:3]=w-v
+    
+    #calculates the matrix for the rotate function
+    def rotate_matrix(self):
+        return matrix([0,-self.norm[2],self.norm[1]],
+                        [self.norm[2],0,-self.norm[0]],
+                        [-self.norm[1],self.norm[0],0])
+    
+    #claculates the cross point of a plane and a line and returns it
+    #also checks wether the line and plane are parallel
+    def cross(self,g):
+        res=True
+        self.A[...,0:1]=-g.xy
+        print(self.A)
+        try:
+            t = np.linalg.solve(self.A, g.origin-self.v)
+        except np.linalg.LinAlgError:  
+            return False,0
+        return g.origin+t*g.xy
+    
+    #calculates the angle between the normal vector of the plane 
+    #and the direction vector of the line
+    def cross_angle(self,g):
+        n_len=length(self.norm)
+        dir_len=length(g.dir)
+        sc=np.dot(g.dir, self.norm)
+        arg=sc/(n_len*dir_len)
+        res=np.rad2deg(np.arccos(arg))
+        if sc<0:
+            res=360-res-180
+        return res
+    
+    #calculates the cross line of two planes
+    #returns the origin and direction vector of the line
+    def cross_planes(self,pl2):
+        cv=np.cross(self.norm,pl2.norm)
+        if np.allclose(cv,0):
+            return 'The planes are parallel'
+        else:
+            A=np.array([[self.norm[0],self.norm[1]],[pl2.norm[0],pl2.norm[1]]])
+            b=np.array([-self.a, -pl2.a])
+            cross_point=np.linalg.solve(A,b)
+            origin=cross_point[0] * self.norm
+        return origin, cv
+    
+    #calculates the reflection vector of a line g and a plane 
+    def reflect(self,g):
+        rot_vec=-rotate_any(self.norm,g.dir,np.pi)
+        if Plane.cross_angle(self,g)==0:
+            rot_vec=-g.dir
+        if np.dot(rot_vec,self.norm)<0:
+            rot_vec=-rot_vec
+        return rot_vec
+        
+    
+    #calcultes the refracted vector of a line g, a plane and the refraction index 
+    #of the scintillator n1 and the surrounding medium n2
+    def refract(self,g,n1,n2):
+        norm_dir=norm_vec(g.dir)
+        refrac_ind=n2/n1
+        angle=Plane.cross_angle(self,g)
+        l=length(g.dir)
+        vec_parallel=np.dot(norm_vec(self.norm),norm_dir)
+        refrac_angle=np.arcsin(refrac_ind * np.sin(angle))
+        vec_orth=[0,0,0]
+        for i in range(3):
+            vec_orth[i]=np.sin(refrac_angle)*norm_dir[i]
+        new_vec_norm=vec_parallel+vec_orth
+        new_vec=-l*new_vec_norm
+        return new_vec
+        
+    def scatter(self,g):
+        u = np.random.uniform()
+        v = np.random.uniform()
+        phi = 2 * np.pi * u
+        theta = np.arccos(2 * v - 1)
+        new_vec = np.array([(np.sin(theta) * np.cos(phi))*40,(np.sin(theta) * np.sin(phi))*40,(np.cos(theta))*40])
+        if np.dot(new_vec, self.norm) < 0:
+            new_vec = -new_vec
+        return new_vec
+ 
+    def line_in_plane(self,g):
+        res=True
+        prod=np.dot(g.dir,self.norm)
+        if prod!=0 or self.point_in_plane(g.origin)!=True:
+            res=False
+        return res
+            
+        
+    def point_in_plane(self,point):
+        distance=np.dot(self.norm,self.o - point)
+        if abs(distance)<self.tolerance:
+            return True
+        else:
+            return False
+    #Attention: the tolerance might need to be higher for certain examples
+    #source of error!
+
+    def __str__(self):
+        return str(self.norm)+str(self.a)+str(self.o).str(self.dir_vec1)
+    
+class Scintillator(list):
+
+    def real_hit(self,g,i):
+        cp=self[i].cross(g)
+        hit=True
+        j=0
+        while j<len(self) and hit==True:
+            if np.dot(self[j].norm,cp-self[j].o)<0:
+                hit=False
+        return hit
+    
+
+#p=Scintillator([Plane(vector(0,0,0),vector(0,-25.981,45),vector(20,0,0)),
+#               Plane(vector(0,0,0),vector(20,0,0),vector(0,51.962,0)),
+#               Plane(vector(0,51.962,0),vector(20,51.962,0),vector(0,77.943,45)),
+#               Plane(vector(0,51.962,90),vector(0,77.943,45),vector(20,51.962,90)),
+#               Plane(vector(0,0,90),vector(0,51.962,90),vector(20,0,90)),
+#               Plane(vector(0,0,90),vector(20,0,90),vector(0,-25.981,45)),
+#               Plane(vector(0,0,0),vector(0,51.962,0),vector(0,0,90)),
+#               Plane(vector(20,0,0),vector(20,0,90),vector(20,51.962,0))])
+
+pl=Plane(vector(0,0,0),vector(20,0,0),vector(0,51.962,0))
+r=rays(3,0.5,vector(10,10,10))
+print(Plane.cross(pl,r))
+
+
+def create_dat(filename, list):
+    with open(filename, 'w') as file:
+        for i in range(len(list)):
+            d=list[i][0][0]
+            e=list[i][1][0]
+            f=list[i][2][0]
+            a,b,c=0,0,0
+            file.write(f"{a}  {b}  {c} {d} {e} {f}\n")
+    
+
+    
+
+
+
+  
+
+        
+
+
+
+
+
+
+
+ 

Visualize.gp

@@ -0,0 +1,29 @@
+#program to visualize the BGO, using the planes norm equation
+#Plane 1: 900*y+519.62*z=0
+#Plane 2: z=0
+#Plane 3: -900*y+519.62*z=-46765.8
+#Plane 4: -900*y-519.62*z=-93531.6
+#Plane 5: -1039.24*z=-93531.6
+#Plane 6: 900*y-519.62*z=-46765.8
+#Plane 7: x=0
+#Plane 8: -4676.58*x=-93531.6
+#[  0.   900.   519.62]0.0
+#[   0.      0.   1039.24]0.0
+#[   0.   -900.    519.62]-46765.8
+#[   0.   -900.   -519.62]-93531.59999999999
+#[    0.       0.   -1039.24]-93531.6
+#[   0.    900.   -519.62]-46765.8
+#[4676.58    0.      0.  ]0.0
+#[-4676.58     0.       0.  ]-93531.6
+
+
+#splot -(900/519.62)*y title "Plane 1", 0 title "Plane 2", (900/519.62)*y-(46765.8/519.62) title "Plane 3", -(900/519.62)*y+(93531.6/519.62) title "Plane 4", (93531.6/1039.24) title "Plane 5", (900/519.62)*y+(46765.8/519.62) title "Plane 6", '/Users/clarapittschellis/Desktop/Light/Light/journey.dat' using 1:2:3:4:5:6 with vectors lw 3
+splot '/Users/clarapittschellis/Desktop/Light/Light/journey.dat' using 1:2:3:4:5:6 with vectors lw 3
+
+set title "BGO"
+set xlabel "X"
+set ylabel "Y"
+set zlabel "Z"
+set xrange[-10:20]
+set yrange[-40:80]
+set zrange[0:100]

Visualize.py

@@ -0,0 +1,64 @@
+import matplotlib.pyplot as plt
+import numpy as np
+from Struktur import Rays, Plane, Ray
+from mpl_toolkits.mplot3d import Axes3D
+
+#vektor = list(map(float, input("Bitte geben Sie die drei Zahlen für den Vektor ein (getrennt durch Leerzeichen): ").split()))
+
+def vis_vector(vektor):
+
+    #if len(vektor) != 3:
+        #print("Bitte geben Sie genau drei Zahlen ein.")
+    #else:
+        #print("Der erstellte Vektor lautet:", vektor)
+
+    ray=Rays(vektor, 1)
+    vec=ray.get_data()
+    ax = plt.figure().add_subplot(projection='3d')
+    u=[]
+    v=[]
+    w=[]
+    for i in range(len(vec)-1):
+        u.append(vec[i][0])
+        v.append(vec[i][1])
+        w.append(vec[i][2])
+    ax.quiver(vektor[0],vektor[1],vektor[2],u,v,w, length=3, normalize=True)
+    ax.set_xlim([-5, 5])
+    ax.set_ylim([-5, 5])
+    ax.set_zlim([-5, 5])
+ 
+
+
+def vis_plane(plane):
+
+    # Define a plane in the form ax + by + cz + d = 0
+    a=plane.norm[0]
+    b=plane.norm[1]
+    c=plane.norm[2]
+    d=plane.a
+
+    # Generate grid points for x and y
+    x = np.linspace(-10, 10, 100)
+    y = np.linspace(-10, 10, 100)
+    x, y = np.meshgrid(x, y)
+
+    # Calculate corresponding z values for the plane
+    z = (-a * x - b * y - d) / c
+
+    # Create a 3D plot
+    fig = plt.figure()
+    ax = fig.add_subplot(111, projection='3d')
+
+    # Plot the surface
+    ax.plot_surface(x, y, z, alpha=0.5)
+
+    # Set labels and title
+    ax.set_xlabel('X')
+    ax.set_ylabel('Y')
+    ax.set_zlabel('Z')
+    ax.set_title('Plane in 3D')
+
+
+
+
+

chaos-jr.gf

@@ -1,26 +0,0 @@
-#! /usr/bin/env ./geometryfactor.py
-# CHAOS Jr.  HET without A-detectors
-
-φ1=0
-φ2=π/6
-Δr=0.1
-Δθ=0.5°
-
-#BGO
-define=W=20
-define=D=40
-prism= 6, ${D}/2*sqrt(3/4), ${W}/2
-
-# SSDs
-define=a=24.3
-define=w=0.3
-define=d=36.5
-define=n=36
-
-prism= ${n}, ${d}/2, ${w}/2
-move= 0,0,-${a}/2
-copy=-1
-move= 0,0,${a}
-
-θ₂=atan2(${d},${a}+${w})
-radius=sqrt(${d}**2 + (${a}+${w})**2)/2

chaos.gf

@@ -1,56 +0,0 @@
-#! /usr/bin/env ./geometryfactor.py
-# CHAOS Jr.  HET without A-detectors
-
-φ1=0
-φ2=π/6
-Δr=0.1
-Δθ=0.5°
-
-#BGO
-define=W=20
-define=D=2*51.962
-prism= 6, ${D}/2*sqrt(3/4), ${W}/2
-
-#Aerogel
-#box = 62,62,40
-#move = 0,0,${W}/2+10
-prism= 4,31,20
-move= 0,0,40
-
-# SSDs
-define=w=0.3
-define=d=36.5
-define=n=36
-define=d1=17.4
-
-prism= ${n}, ${d}/2, ${w}/2
-#A
-move= 0,0,${W}/2+55
-#C
-copy=-1
-move= 0,0,-50
-#E
-copy=-1
-move= 0,0,-30
-
-#E123
-copy=-1
-move=-18.5,-10.53,-4
-copy=-1
-move=2*18.5,0,-4
-copy=-1
-move=-18.25,31.6,-4
-
-#A1
-prism= ${n}, ${d1}/2, ${w}/2
-move= 0,0,${W}/2+55
-#C1
-copy=-1
-move= 0,0,-50
-
-#?
-define=a=24.7
-θ₂=atan2(${d},${a}+${w})
-radius=sqrt(${d}**2 + (${a}+${w})**2)/2
-
-

geometryfactor.py

@@ -1,10 +1,9 @@
 #! /usr/bin/python3
 
 from sys import argv, stdout, stderr
-from math import pi as π, floor
-from math import atan, atan2
+from math import pi as π
 import numpy
-from numpy import sqrt, sin, cos, arccos, tan, array, arange, empty, cross, newaxis
+from numpy import sqrt, sin, cos, tan, array, arange, empty, cross, newaxis
 from numpy.linalg import solve
 from getopt import getopt
 
@@ -61,13 +60,7 @@ class ray:
 class rays:
     "Fast reimplementation of a beam of rays"
     
-    def __init__(self, r, Δr, spiral=False):
-        if spiral:
-            self.spiral(r, Δr)
-        else:
-            self.rings(r, Δr)
-
-    def rings(self, r, Δr):
+    def __init__(self, r, Δr):
         nr = int(r/Δr + 1.001)
         self.n = 0
         rψ = 0
@@ -94,16 +87,6 @@ class rays:
                 self.xy[i,:,0] = (r*cos(ψ), r*sin(ψ), 0)
                 i = i+1
 
-    def spiral(self, r, Δr):
-        self.n = int(π * (r/Δr)**2 + 0.5)
-        self.ψ = sqrt(4*π * arange(self.n))
-        self.r = Δr/(2*π) * self.ψ
-        self.ψ %= 2*π
-        self.xy = empty((self.n, 3, 1))
-        self.xy[:,0,0] = self.r * cos(self.ψ)
-        self.xy[:,1,0] = self.r * sin(self.ψ)
-        self.xy[:,2,0] = 0
-
     def bend(self, θ, φ):
         self.θ = θ
         self.φ = φ
@@ -121,7 +104,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 +126,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 self.solution
+            return -1, None
+        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]
@@ -279,7 +247,7 @@ class detector(list):
         if edjes is None:
             edjes = self.edges()
         for e in edjes:
-            stdout.write("%g %g %g\n%g %g %g\n\n\n" % (tuple(e[0][:,0])+tuple(e[1][:,0])))
+            stdout.write("%g %g %g\n%g %g %g\n\n\n" % (tuple(e[0])+tuple(e[1])))
 
     def planes(self):
         # is this useful?
@@ -301,23 +269,12 @@ def box(x=1.0, y=1.0, z=1.0):
         plane(vector(y= y), vector(x=1), vector(z=1)),
     ])
 
-def prism(n=6, a=1.0, z=1.0, φ=0.0, n1=0, n2=None):
+def prism(n=6, a=1.0, z=1.0, φ=0.0):
     p = [
         plane(vector(z=-z), vector(x=1), vector(y=1)),
         plane(vector(z= z), vector(y=1), vector(x=1)),
     ]
-    if n2 is None:
-        n2 = n-1
-    else:
-        c = cos(2*π/n * n1)
-        s = sin(2*π/n * n1)
-        p.append(plane(vector(), vector(c,s), vector(z=-1)))
-        c = cos(2*π/n * n2)
-        s = sin(2*π/n * n2)
-        p.append(plane(vector(), vector(c,s), vector(z= 1)))
- 
-    φ += 2*π/n * arange(int(floor(n1+0.4999)),
-                        int(floor(n2+1.5001)) )
+    φ = 2*π * arange(int(n))/n + φ
     p.extend([
         plane(
             vector(a*c, a*s),
@@ -330,25 +287,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)
@@ -374,7 +312,7 @@ def run(world, r, Δr, Δθ, θ1=0, θ2=π/2, φ1=0, φ2=2*π):
                     pl = [d.pathlength(ry) for d in world]
                     stdout.write(fmt % ((θ, φ, r, ψ)+tuple(pl)))
 
-def run_rings(world, r, Δr, Δθ, θ1=0, θ2=π/2, φ1=0, φ2=2*π):
+def run_beam(world, r, Δr, Δθ, θ1=0, θ2=π/2, φ1=0, φ2=2*π):
     fmt = "%g %g %g %g "+len(world)*" %g"+"\n"
     beam = rays(r, Δr)
     b = None
@@ -394,24 +332,6 @@ def run_rings(world, r, Δr, Δθ, θ1=0, θ2=π/2, φ1=0, φ2=2*π):
                 ll = tuple(bb+[l[i] for l in pl])
                 stdout.write(fmt % ll)
         
-def run_beam(world, r, Δr, Δθ, θ1=0, θ2=π/2, φ1=0, φ2=2*π):
-    if φ2 < 2*π - Δθ or φ1 > Δθ:
-        run_rings(world, r, Δr, Δθ, θ1, θ2, φ1, φ2)
-    fmt = "%g %g %g %g "+len(world)*" %g"+"\n"
-    beam = rays(r, Δr, spiral=True)
-    N = 2*π / Δθ**2
-    t1 = N*(1-cos(θ1))
-    n = int(N*(cos(θ1) - cos(θ2)) + 0.5)
-    t = arange(n) + t1
-    for θ in arccos(1 - t/N):
-        φ = (2*π * θ/Δθ) % (2*π)
-        beam.bend(θ, φ)
-        pl = [d.pathlength_beam(beam) for d in world]
-        for i in range(beam.n):
-            bb=[θ, φ, beam.r[i], beam.ψ[i]]
-            ll = tuple(bb+[l[i] for l in pl])
-            stdout.write(fmt % ll)
-
 test_rays = [
     ("x-axis", ray().ov(vector(),vector(x=1))),
     ("y-axis", ray().ov(vector(),vector(y=1))),
@@ -440,7 +360,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=",
@@ -474,13 +394,12 @@ world = [detector()]
 θ1=0
 θ2=π/2
 φ1=0
-φ2=2*π
+φ2=2*π/8
 do_size = False
 do_test = False
 do_run = False
 do_corners = False
 do_gnuplot = False
-do_equiv = False
 rot = rotate(0,0)
 ori = vector()
 
@@ -543,8 +462,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 +502,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:
@@ -606,10 +519,11 @@ if do_gnuplot:
         d.gnuplot()
 
 if do_size:
-    no = int(π * (r/Δr)**2 + 0.5)
-    nv = int((φ2 - φ1) / Δθ**2 * (cos(θ1) - cos(θ2)) + 0.5)
-    sr = nv * Δθ**2
-    A = π * r**2
+    nr = int(r/Δr+1.001)
+    no = π*nr**2 - (π-0.5)*nr
+    sr = (φ2-φ1)*(cos(θ1)-cos(θ2+Δθ))
+    A = π*((nr-0.5)*Δr)**2
+    nv = sr/Δθ**2
     nn = int(nv)*int(no)
     stderr.write("""Run size estimate:
 detectors: %d  planes: %d

journey.dat

@@ -0,0 +1,65 @@
+0  0  0 -0.0 0.0 0.0
+0  0  0 0.43301270189221935 0.24999999999999997 0.0
+0  0  0 3.061616997868383e-17 0.5 0.0
+0  0  0 -0.43301270189221924 0.25000000000000017 0.0
+0  0  0 -0.4330127018922194 -0.2499999999999999 0.0
+0  0  0 -9.184850993605148e-17 -0.5 0.0
+0  0  0 0.4330127018922192 -0.2500000000000002 0.0
+0  0  0 0.970941817426052 0.23931566428755774 0.0
+0  0  0 0.7485107481711012 0.6631226582407952 0.0
+0  0  0 0.3546048870425358 0.9350162426854147 0.0
+0  0  0 -0.1205366802553229 0.992708874098054 0.0
+0  0  0 -0.5680647467311556 0.8229838658936566 0.0
+0  0  0 -0.8854560256532096 0.46472317204376906 0.0
+0  0  0 -1.0 5.66553889764798e-16 0.0
+0  0  0 -0.8854560256532101 -0.46472317204376806 0.0
+0  0  0 -0.5680647467311559 -0.8229838658936564 0.0
+0  0  0 -0.12053668025532356 -0.992708874098054 0.0
+0  0  0 0.35460488704253473 -0.9350162426854152 0.0
+0  0  0 0.7485107481711009 -0.6631226582407955 0.0
+0  0  0 0.9709418174260518 -0.23931566428755865 0.0
+0  0  0 1.4815325108927067 0.23465169756034632 0.0
+0  0  0 1.336509786282552 0.6809857496093201 0.0
+0  0  0 1.0606601717798214 1.0606601717798212 0.0
+0  0  0 0.6809857496093202 1.3365097862825517 0.0
+0  0  0 0.23465169756034637 1.4815325108927067 0.0
+0  0  0 -0.2346516975603462 1.4815325108927067 0.0
+0  0  0 -0.6809857496093201 1.336509786282552 0.0
+0  0  0 -1.0606601717798212 1.0606601717798214 0.0
+0  0  0 -1.3365097862825517 0.6809857496093203 0.0
+0  0  0 -1.4815325108927064 0.23465169756034648 0.0
+0  0  0 -1.4815325108927067 -0.23465169756034615 0.0
+0  0  0 -1.336509786282552 -0.68098574960932 0.0
+0  0  0 -1.0606601717798214 -1.0606601717798212 0.0
+0  0  0 -0.6809857496093203 -1.3365097862825517 0.0
+0  0  0 -0.2346516975603466 -1.4815325108927064 0.0
+0  0  0 0.234651697560346 -1.4815325108927067 0.0
+0  0  0 0.6809857496093199 -1.336509786282552 0.0
+0  0  0 1.0606601717798212 -1.0606601717798214 0.0
+0  0  0 1.3365097862825517 -0.6809857496093205 0.0
+0  0  0 1.4815325108927064 -0.23465169756034668 0.0
+0  0  0 1.9842294026289558 0.2506664671286085 0.0
+0  0  0 1.8595529717765027 0.736249105369356 0.0
+0  0  0 1.618033988749895 1.1755705045849463 0.0
+0  0  0 1.2748479794973793 1.5410264855515785 0.0
+0  0  0 0.8515585831301453 1.8096541049320392 0.0
+0  0  0 0.3747626291714494 1.9645745014573772 0.0
+0  0  0 -0.1255810390586268 1.9960534568565431 0.0
+0  0  0 -0.6180339887498951 1.902113032590307 0.0
+0  0  0 -1.0716535899579935 1.68865585100403 0.0
+0  0  0 -1.4579372548428233 1.369094211857377 0.0
+0  0  0 -1.7526133600877274 0.9635073482034304 0.0
+0  0  0 -1.9371663222572624 0.49737977432970876 0.0
+0  0  0 -2.0 -6.432490598706546e-16 0.0
+0  0  0 -1.9371663222572622 -0.49737977432971003 0.0
+0  0  0 -1.7526133600877265 -0.9635073482034315 0.0
+0  0  0 -1.4579372548428222 -1.3690942118573781 0.0
+0  0  0 -1.0716535899579926 -1.6886558510040306 0.0
+0  0  0 -0.6180339887498935 -1.9021130325903075 0.0
+0  0  0 -0.1255810390586264 -1.9960534568565431 0.0
+0  0  0 0.37476262917145026 -1.9645745014573772 0.0
+0  0  0 0.8515585831301453 -1.8096541049320392 0.0
+0  0  0 1.27484797949738 -1.541026485551578 0.0
+0  0  0 1.6180339887498958 -1.1755705045849452 0.0
+0  0  0 1.859552971776503 -0.7362491053693556 0.0
+0  0  0 1.9842294026289558 -0.2506664671286076 0.0

norm.dat

@@ -0,0 +1,9 @@
+10 -10 17.321 0 900 519.62
+#10 20 0 0 0 1039.24
+#10 61.962 17.321 0 -900 519.62
+#10 61.962 62.321 0 -900 -519.62
+#10 20 90 0 0 -1039.24
+#10 -10 62.321 0 900 -519.62
+#0 20 45 4676.58 0 0
+#20 20 45 -4676.58 0 0
+

seth.gf

@@ -1,110 +0,0 @@
-#! /usr/bin/env ./geometryfactor.py
-# CHAOS Jr.  HET without A-detectors
-
-Δr=1
-Δθ=5°
-radius= 66
-
-#BGO
-define=W=20
-define=D=2*51.962
-prism= 6, ${D}/2*sqrt(3/4), ${W}/2
-
-move= 0, 0, -15
-copy= -1
-move= 0, 0,  30
-
-box= 10, 5, 0.15
-move= 0, -49.7, -30.25
-copy= -1
-move= 0, 14.2, 0
-copy= -1
-move= 0, 14.2, 0
-copy= -1
-move= 0, 14.2, 0
-copy= -1
-move= 0, 14.2, 0
-copy= -1
-move= 0, 14.2, 0
-copy= -1
-move= 0, 14.2, 0
-copy= -1
-move= 0, 14.2, 0
-
-copy= -7
-move= -27, 0,0
-copy= -7
-move= -27, 0,0
-copy= -7
-move= -27, 0,0
-copy= -7
-move= -27, 0,0
-copy= -7
-move= -27, 0,0
-copy= -7
-move= -27, 0,0
-
-copy= -6
-move= 54
-copy= -6
-move= 54
-copy= -6
-move= 54
-copy= -6
-move= 54
-copy= -6
-move= 54
-copy= -6
-move= 54
-
-copy= -20
-move= 0,0, 60.5
-copy= -20
-move= 0,0, 60.5
-copy= -20
-move= 0,0, 60.5
-copy= -20
-move= 0,0, 60.5
-copy= -20
-move= 0,0, 60.5
-copy= -20
-move= 0,0, 60.5
-copy= -20
-move= 0,0, 60.5
-copy= -20
-move= 0,0, 60.5
-copy= -20
-move= 0,0, 60.5
-copy= -20
-move= 0,0, 60.5
-copy= -20
-move= 0,0, 60.5
-copy= -20
-move= 0,0, 60.5
-copy= -20
-move= 0,0, 60.5
-copy= -20
-move= 0,0, 60.5
-copy= -20
-move= 0,0, 60.5
-copy= -20
-move= 0,0, 60.5
-copy= -20
-move= 0,0, 60.5
-copy= -20
-move= 0,0, 60.5
-copy= -20
-move= 0,0, 60.5
-copy= -20
-move= 0,0, 60.5
-
-# SSDs
-define=w=0.5
-define=d=36.5
-define=n=36
-define=d1=17.4
-
-prism= ${n}, ${d}/2, ${w}/2
-move= 0,0, -${w}/2
-prism= ${n}, ${d1}/2, ${w}/2
-move= 0,0, ${w}/2

test.py

@@ -0,0 +1,34 @@
+import random
+from Random_gen import test_planes,test_lines
+from Struktur import Plane
+
+def test_cross(test_size):
+    rand_planes=test_planes(test_size)
+    rand_lines=test_lines(test_size)
+    cross_point=[]
+    for i in range(len(rand_planes)):
+        cross_point_fun=Plane.cross(rand_planes[i],rand_lines[i])
+        if cross_point!=str:
+            cross_point.append(cross_point_fun)
+        if Plane.point_in_plane(rand_planes[i],cross_point)==False:
+            raise Exception('Catastrophe!')
+    return cross_point
+
+def test_cross_angle(test_size):
+    rand_planes=test_planes(test_size)
+    rand_lines=test_lines(test_size)
+    cross_angle=[]
+    for i in range(len(rand_planes)):
+        angle=Plane.cross_angle(rand_planes[i],rand_lines[i])
+        cross_angle.append(angle)
+    return cross_angle
+
+ 
+print(test_cross_angle(60))
+
+
+
+        
+
+
+

try.py

@@ -0,0 +1,35 @@
+import numpy as np
+
+A=np.array([[[2,2,2],[4,4,4]],[[4,5,6],[7,8,9]]])
+b=np.array([6,6,6])
+c=np.array(6)
+#print(A.shape)
+#print(b.shape)
+#print(c.shape)
+#print(A*b)
+b=b[np.newaxis,np.newaxis,:]
+#print(b.shape)
+#print(b)
+b=np.tile(b,(2,2,1))
+#print(b.shape)
+#print(A*b)
+
+P=np.zeros((3,3,3))
+#print(P)
+n=np.array([[4],[5],[6]])
+n=np.squeeze(n.reshape(1,-1))
+print(n)
+P[:,1]=n.T
+print(P)
+P[...,2:3]=np.array([[7],[8],[9]])
+#print(P[...,0:1])
+n=np.array([[[2],[1],[3]],[[6],[4],[3]],[[8],[4],[0]]])
+#print(n)
+P[...,0:1]=n
+v=np.array([4,9,2])
+v=v[np.newaxis,:]
+v=np.tile(v,(3,1))
+#print(P)
+#print(v.shape)
+#print(P.shape)
+#print(np.linalg.solve(P,v-np.array([2,4,3])))