% N = number of particles
% Q = tune
% R = transfer matrix
% X = gaussian bunch

% lattice
Q = 0.2516;
phi = 2*pi*Q;
R = [ cos(phi) sin(phi) ; -sin(phi) cos(phi) ];

% bunch with N = 10000 particles
N = 10000;
X = 0.2 * randn(2,N);

% CASE 1: simulate 1000 turns with only linear optics
%for i=1:1000
%  X = R * X;
%end

% CASE 2: simulate 1000 turns with a sexupolar kick each turn
for i=1:1000
  X = R * (X + [ zeros(1,N)
                 X(1,:) .** 2 ] );
end

% plot
f = figure;
set(f, "visible", "off");

plot(X(1,:), X(2,:), '.');
axis([ -1 1 -1 1 ]);
xlabel('x [mm]');
ylabel('x'' [mrad]');

% save the plot as a PDF file
print(f, 'SextupolarKick.pdf', '-dpdf', '-FTimes New Roman:16');
