ocelot-collab / ocelot

OCELOT is a multiphysics simulation toolkit designed for studying FEL and storage ring-based light sources.
GNU General Public License v3.0
85 stars 58 forks source link

Octupole element for particle tracking #201

Closed vitben closed 8 months ago

vitben commented 11 months ago

Hello,

thanks for developing such a good code! I am performing particle tracking of a 6D distribution. The beamline have a set of sextupole and octupole elements. When performing the tracking, I noticed that there is not difference in the case when I include the octupoles and when I do not. Here below you can find the python code for the tracking with and without octupoles.

`# initialization of tracking method
%autoreload
from my_lattice import *
method = {"global": SecondTM}

# for first order tracking uncomment next line
# method = {"global": TransferMap}

lat_t1 = MagneticLattice(cell1, start=None, stop=None, method=method)

navi = Navigator(lat_t1)
p_array = deepcopy(p_array_init)
start = time.time()
tws_track1, p_array1 = track(lat_t1, p_array, navi, calc_tws=True)
print("\n time exec:", time.time() - start, "sec")

lat_t2 = MagneticLattice(cell2, start=None, stop=None, method=method)

navi = Navigator(lat_t2)
p_array = deepcopy(p_array_init)
start = time.time()
tws_track2, p_array2 = track(lat_t2, p_array, navi, calc_tws=True)
print("\n time exec:", time.time() - start, "sec")

tau1 = np.array([p.tau for p in p_array1])*1000
dp1 = np.array([p.p for p in p_array1])*1000
x1 = np.array([p.x for p in p_array1])*1000
xp1 = np.array([p.px for p in p_array1])*1000
y1 = np.array([p.y for p in p_array1])*1000
yp1 = np.array([p.py for p in p_array1])*1000

tau2 = np.array([p.tau for p in p_array2])*1000
dp2 = np.array([p.p for p in p_array2])*1000
x2 = np.array([p.x for p in p_array2])*1000
xp2 = np.array([p.px for p in p_array2])*1000
y2 = np.array([p.y for p in p_array2])*1000
yp2 = np.array([p.py for p in p_array2])*1000

The octupoles elements I use are defined as follow:

O_1=  Octupole(l = 0.2 , k3 =-328.29496694078165, tilt = 0, eid = 'O_1')
O_2= Octupole(l = 0.2, k3 =-208.57373589036763, tilt = 0, eid = 'O_2')
O_3=   Octupole(l = 0.2, k3 =2844.666880229662, tilt = 0, eid = 'O_3')

Do you have any ideas of what the issue could be?

Thanks a lot for your help!

sergey-tomin commented 11 months ago

hi, the issue is that Octuple is 3d order element and by default we activate first order maps for all elements. To make Octuple work you need to specify for it a different map. method = {"global": SecondTM, Octupole: KickTM}

see more details in tutorial. https://nbviewer.org/github/ocelot-collab/ocelot/blob/dev/demos/ipython_tutorials/small_useful_features.ipynb#TM

Sergey.

vitben commented 11 months ago

Hi Sergey,

Thanks for the very quick answer. Now it works perfectly.

Best,

Vittorio