Geometry
The code presented in this section is a Python3 re-implementation of the methods described in the book of D.Sunday Pratical Geometry Algorithms. It implements class that represent and compute operations on vectors, lines, points and rays.
This algorithm is used in Zgoubidoo where we must “re-align” the tracks on a reference one to compute Twiss parameters. The main function in this module is:
project_on_reference(ref: ReferenceTrajectory, trajectories: list):
This method takes as input a Reference trajectory and tracks to interpolate and returns the new coordinates for the tracks. A Reference trajectory is a class that take as input a numpy.ndarray where columns are x,y,z,t,p and s.
The example below shows an example where, at each point of the reference tracks (in blue, we find the perpendicular points on given tracks (in red and dark orange).
%matplotlib inline
import matplotlib.pyplot as plt
import numpy as _np
import georges_core
import georges_core.geometry as ggeom
x1 = _np.arange(0,11, 0.7)
y1 = _np.sin(x1/2)/2
t1 = _np.cos(x1/2)/4
dist = _np.linalg.norm(_np.array((x1, y1)).T - _np.roll(_np.array((x1, y1)).T, -1, axis=0),axis=1)
s1 = _np.concatenate((_np.array([0]), _np.cumsum(dist[0:-1])))
x2 = _np.arange(0,11, 0.3)
y2 = _np.sin(x2)
t2 = _np.cos(x2)
dist = _np.linalg.norm(_np.array((x2, y2)).T - _np.roll(_np.array((x2, y2)).T, -1, axis=0),axis=1)
s2 = _np.concatenate((_np.array([0]), _np.cumsum(dist[0:-1])))
x3 = _np.arange(0,11, 0.3)
y3 = _np.cos(x3)
t3 = -_np.sin(x3)
dist = _np.linalg.norm(_np.array((x3, y3)).T - _np.roll(_np.array((x3, y3)).T, -1, axis=0),axis=1)
s3 = _np.concatenate((_np.array([0]), _np.cumsum(dist[0:-1])))
ref = _np.array([x1, y1, _np.zeros(len(x1)), t1, _np.zeros(len(x1)), s1]).T
traj2 = _np.zeros([len(x2), 6])
traj2[:,0] = x2
traj2[:,1] = y2
traj2[:,3] = t2
traj2[:,5] = s2
traj3 = _np.zeros([len(x3), 6])
traj3[:,0] = x3
traj3[:,1] = y3
traj3[:,3] = t3
traj3[:,5] = s3
ref_traj = ggeom.ReferenceTrajectory(points=ggeom.Points(ref))
interpolated_points = ggeom.project_on_reference(ref=ref_traj, trajectories = [ggeom.Trajectories(traj2),
ggeom.Trajectories(traj3)])
def get_lines_points(point, normal):
if normal[1] == 0:
return [[point[0], point[0]],[point[1]-1, point[1]+1]]
else:
xa = point[0] - 0.5
xb = point[0] + 0.5
ya = point[1]-((normal[0]/normal[1])*(xa - point[0]))
yb = point[1]-((normal[0]/normal[1])*(xb - point[0]))
return [[xa,xb],[ya,yb]]
fig, axs = plt.subplots(figsize=(10, 10), nrows=2, ncols=1, sharex=False)
axs[0].plot(ref_traj.data[:, 0], ref_traj.data[:, 1], '-b*')
axs[0].plot(x2, y2, ls='-', color='r', marker='x')
axs[0].plot(interpolated_points[0,:, 0], interpolated_points[0,:, 1], 'ok', fillstyle='none')
for k in list(map(ggeom.Plane, ref_traj.frenet_frames().data)):
xl, yl = get_lines_points(k.point, k.normal)
axs[0].plot(xl, yl, ls='--', lw=0.5, color='k')
axs[0].set_xticks(_np.arange(0, 11+1, 1.0))
axs[0].set_yticks(_np.arange(-3, 4, 1.0))
axs[0].axis('equal')
axs[0].set_ylabel('y')
axs[0].set_ylim([-1.5,1.5])
axs[0].grid(True)
axs[1].plot(ref_traj.data[:, 0], ref_traj.data[:, 1], '-b*')
axs[1].plot(x3, y3, ls='-', color='darkorange', marker='x')
axs[1].plot(interpolated_points[1,:, 0], interpolated_points[1,:, 1], 'ok', fillstyle='none')
for k in list(map(ggeom.Plane, ref_traj.frenet_frames().data)):
xl, yl = get_lines_points(k.point, k.normal)
axs[1].plot(xl, yl, ls='--', lw=0.5, color='k')
axs[1].set_xticks(_np.arange(0, 11+1, 1.0))
axs[1].set_yticks(_np.arange(-3, 4, 1.0))
axs[1].axis('equal')
axs[1].set_ylabel('y')
axs[1].set_ylim([-1.5,1.5])
axs[1].grid(True)