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')
+
+
+
+
+

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
+

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])))