Commit e31a2da7 authored by Yoann Sallaz Damaz's avatar Yoann Sallaz Damaz
Browse files

Upload raytracing beamline simu

parent cfd17bb3
Pipeline #41733 canceled with stages
from shapely.geometry import Point, LineString, MultiPoint
from shapely.affinity import rotate
from shapely.ops import nearest_points, substring
from math import cos, sin, atan, atan2, degrees, radians, pi
import matplotlib.pyplot as plt
import numpy as np
#################################################################
# DRAW TO DRAW THE BEAM
DRAW = True
# SCAN TO SCAN AND PLOT DIODE COUNT DURING SCAN
SCAN = True
# NB beam in simu (2 to see beam size, 200 is good for scan)
NB_BEAM = 2
# APPERTURE OF SL0
APERTURE_SL0 = 5 #mm
#################################################################
#################################################################
# ERROR ALIGNEMENT
ERROR_ALT_SL0 = 0 #mm
ERROR_ALT_M1 = 0 #mm
ERROR_ALT_MOVEH = 0 #mm
ERROR_ALT_SOURCE = 0 #mm
ERROR_ANG_SOURCE = 0 #rad
#################################################################
NOT_INTER_DISTANCE = 123456789
def extendline_from_direction(orig, angle, dist):
Xend = orig.x+dist*cos(angle)
Yend = orig.y+dist*sin(angle)
return Point(Xend,Yend)
def angle_reflechi(Psource,Pcontact,Psurf1,Psurf2):
angleSource = atan2(Pcontact.y-Psource.y,Pcontact.x-Psource.x)
angleSurface = atan2(Psurf2.y-Psurf1.y,Psurf2.x-Psurf1.x)
#print("angleSource=",degrees(angleSource))
#print("angleSurface=",degrees(angleSurface))
return 2*angleSurface-angleSource
class objet_optique():
def __init__(self, name, liste_coord, centre_rot, angle, shift=0, reflect=False):
self.name = name
self.geo = LineString(liste_coord)
self.centre = centre_rot
self.geo=rotate(self.geo, angle, origin=centre_rot, use_radians=False)
self.reflect = reflect
def intercept(self, begin, angle):
end = extendline_from_direction(begin, angle, 100000)
line = LineString([begin, end])
intersection = line.intersection(self.geo)
if intersection.is_empty:
return None, NOT_INTER_DISTANCE
if not isinstance(intersection, Point):
intersection = list(intersection)[0]
return intersection, LineString([begin, intersection]).length
def get_surface_facette(self, intersection):
ratio_intersection = (intersection.x-self.geo.bounds[0])/(self.geo.bounds[-2]-self.geo.bounds[0])
nearest_geoms = list(substring(self.geo, ratio_intersection-0.002, ratio_intersection+0.002, normalized=True).coords)
return Point(nearest_geoms[0]), Point(nearest_geoms[-1])
# SOURCE
P_Source = Point(0, ERROR_ALT_SOURCE)
Angle_Source = ERROR_ANG_SOURCE
positionList = []
diodeList = []
if DRAW:
scanlist = [0]
SCAN=False
else:
scanlist = list(range(-40,42,2))
for shift in scanlist:
shift/=10.
AllObject = []
#############################
# M1
#############################
courbeM1 = 14.31E6 #mm
lengthM1 = 1280 #mm
heightM1 = 90 #mm
P_M1 = Point(29601, shift+ERROR_ALT_M1) #mm
#M1 partie reflechissante
demiAngleArc = atan(0.5*lengthM1/courbeM1)
listePoint=[]
altituderef = cos(demiAngleArc)*courbeM1
for angle in np.arange(-demiAngleArc, demiAngleArc+demiAngleArc*2/100, demiAngleArc*2/100):
listePoint.append(Point(P_M1.x+sin(angle)*courbeM1, P_M1.y+altituderef-cos(angle)*courbeM1))
AllObject.append(objet_optique("M1_refl", listePoint, P_M1, 0.237,reflect=True))
#M1 corp non reflechissant
first=listePoint[0]
end = listePoint[-1]
listePoint=[first, Point(first.x, first.y-heightM1),
Point(end.x,end.y-heightM1), end]
AllObject.append(objet_optique("M1_corps", listePoint, P_M1, 0.237,reflect=False))
#############################
# MASQUE M1
#############################
P_MasqueM1 = Point(28882,shift+ERROR_ALT_M1)
appertureMasqueM1 = 6.8 #mm
listePoint = [Point(P_MasqueM1.x, 0.5*appertureMasqueM1), Point(P_MasqueM1.x, 0.5*appertureMasqueM1+100)]
AllObject.append(objet_optique("Masque_M1_haut", listePoint, P_MasqueM1, 0, reflect=False))
listePoint = [Point(P_MasqueM1.x, -0.5*appertureMasqueM1), Point(P_MasqueM1.x, -0.5*appertureMasqueM1-100)]
AllObject.append(objet_optique("Masque_M1_bas", listePoint, P_MasqueM1, 0, reflect=False))
#############################
# C1
#############################
P_C1 = Point(32336, 22.62+shift+ERROR_ALT_MOVEH) #mm
lengthC1 = 50 #mm
heightC1 = 50 #mm
#C1 partie reflechissante
listePoint = [Point(P_C1.x-0.5*lengthC1, P_C1.y), Point(P_C1.x+0.5*lengthC1, P_C1.y)]
AllObject.append(objet_optique("C1_refl", listePoint, P_C1, -8.514, reflect=True))
#C1 corp non reflechissant
listePoint = [Point(P_C1.x-0.5*lengthC1, P_C1.y), Point(P_C1.x-0.5*lengthC1, P_C1.y+heightC1),
Point(P_C1.x+0.5*lengthC1, P_C1.y+heightC1), Point(P_C1.x+0.5*lengthC1, P_C1.y)]
AllObject.append(objet_optique("C1_corps", listePoint, P_C1, -8.514, reflect=False))
#############################
# C2
#############################
P_C2 = Point(P_C1.x+208.70, P_C1.y-33.01) #mm
lengthC2 = 50 #mm
heightC2 = 50 #mm
#C2 partie reflechissante
listePoint = [Point(P_C2.x-0.5*lengthC2, P_C2.y), Point(P_C2.x+0.5*lengthC2, P_C2.y)]
AllObject.append(objet_optique("C2_refl", listePoint, P_C1, -8.514, reflect=True))
#C2 corp non reflechissant
listePoint = [Point(P_C2.x-0.5*lengthC2, P_C2.y), Point(P_C2.x-0.5*lengthC2, P_C2.y-heightC2),
Point(P_C2.x+0.5*lengthC2, P_C2.y-heightC2), Point(P_C2.x+0.5*lengthC2, P_C2.y)]
AllObject.append(objet_optique("C2_corps", listePoint, P_C1, -8.514, reflect=False))
#############################
# DIODE
#############################
P_Diode = Point(33550,0)
listePoint = [Point(P_Diode.x, P_Diode.y-50), Point(P_Diode.x,P_Diode.y+0)]
AllObject.append(objet_optique("Diode", listePoint, P_Diode, 0, reflect=False))
#############################
# MASQUE FE
#############################
P_MasqueFE = Point(23437,0)
appertureMasqueFE = 8 #mm
listePoint = [Point(P_MasqueFE.x, P_MasqueFE.y+0.5*appertureMasqueFE),
Point(P_MasqueFE.x, P_MasqueFE.y+0.5*appertureMasqueFE+100)]
AllObject.append(objet_optique("Masque_FE_haut", listePoint, P_MasqueFE, 0, reflect=False))
listePoint = [Point(P_MasqueFE.x, P_MasqueFE.y-0.5*appertureMasqueFE),
Point(P_MasqueFE.x, P_MasqueFE.y-0.5*appertureMasqueFE-100)]
AllObject.append(objet_optique("Masque_FE_bas", listePoint, P_MasqueFE, 0, reflect=False))
#############################
# MASQUE BE
#############################
P_MasqueBE = Point(30577,0)
appertureMasqueBE = 20 #mm A VERIF
listePoint = [Point(P_MasqueBE.x, P_MasqueBE.y+appertureMasqueBE),
Point(P_MasqueBE.x, P_MasqueBE.y+appertureMasqueBE+100)]
AllObject.append(objet_optique("Masque_BE_haut", listePoint, P_MasqueBE, 0, reflect=False))
listePoint = [Point(P_MasqueBE.x, P_MasqueBE.y),
Point(P_MasqueBE.x, P_MasqueBE.y-100)]
AllObject.append(objet_optique("Masque_BE_bas", listePoint, P_MasqueBE, 0, reflect=False))
#############################
# SL0
#############################
P_SL0 = Point(23437+1000,shift+ERROR_ALT_SL0) #AVERIF
appertureSL0 = APERTURE_SL0 #mm
listePoint = [Point(P_SL0.x, P_SL0.y+0.5*appertureSL0), Point(P_SL0.x, P_SL0.y+0.5*appertureSL0+100)]
AllObject.append(objet_optique("SL0_haut", listePoint, P_SL0, 0, reflect=False))
listePoint = [Point(P_SL0.x, P_SL0.y-0.5*appertureSL0), Point(P_SL0.x, P_SL0.y-0.5*appertureSL0-100)]
AllObject.append(objet_optique("SL0_bas", listePoint, P_SL0, 0, reflect=False))
#draw All object
if DRAW :
for obj in AllObject:
x,y = obj.geo.xy
if obj.reflect:
col = "b"
else:
col = "k"
plt.plot(x,y,col)
vert_div_source = 0.174E-3 # radian
step = vert_div_source/(NB_BEAM-1)
nbHit = 0
for angle in np.arange(-0.5*vert_div_source+Angle_Source, 0.5*vert_div_source+Angle_Source+step, step):
AllObjectCopy = AllObject.copy()
rayon = [P_Source]
dead = False
while(not dead):
liste_intersection=[]
liste_distance=[]
for obj in AllObjectCopy:
intersection, distance = obj.intercept(rayon[-1], angle)
liste_intersection.append(intersection)
liste_distance.append(distance)
distance_min = min(liste_distance)
if distance_min == NOT_INTER_DISTANCE:
dead=True
rayon.append(extendline_from_direction(rayon[-1], angle, 1000))
else:
idx = liste_distance.index(distance_min)
obj = AllObjectCopy[idx]
intersection = liste_intersection[idx]
#print("\trayon hit",obj.name,"at",intersection)
if obj.reflect:
objSurf_P1, objSurf_P2 = obj.get_surface_facette(intersection)
new_angle = angle_reflechi(rayon[-1],intersection,objSurf_P1,objSurf_P2)
angle=new_angle
else:
dead=True
if obj.name=="Diode":
nbHit+=1
rayon.append(intersection)
del(AllObjectCopy[idx])
if DRAW :
full_ray = LineString(rayon)
x,y = full_ray.xy
plt.plot(x,y,"g")
print("SHIFT, NB HIT IN DIODE = ", shift, nbHit)
positionList.append(shift)
diodeList.append(nbHit)
if DRAW :
plt.axis('equal')
plt.show()
if SCAN:
plt.plot(positionList,diodeList, marker = 'o')
plt.xlabel("position du scan (mm)")
plt.ylabel("Nb hit diode")
plt.show()
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment