#Calculate the normal tunes for a one-turn map and beam-beam kick

# the beam-beam parameters

xi1 = -0.0107  # beam 1-2
xi2 = -0.0107  # beam 2-1

beta1 = 0.28  # beam 1 beta twiss
beta2 = 0.28  # beam 2 beta twiss
alpha1 = 0.01 # beam1 alpha twiss
alpha2 = -alpha1 # beam 2 alpha twiss

tune1 = 0.428;  # beam 1 base tune
tune2 = 0.428;  # beam 2 base tune

mu1 = pi * tune1;  # beam 1 half  phase advance
mu2 = pi * tune2;  # beam 2 half phase advance

				# thin beambeam kick

k1 = -2 * pi * xi1/beta1;
k2 = -2 * pi * xi2/beta2;

bbkick = \
    [ 1, 0, 0, 0; \
     -k1, 1, k1, 0; \
     0, 0, 1, 0; \
     k2, 0, -k2, 1]

gamma1 = (1 + alpha1**2)/beta1;
gamma2 = (1 + alpha2**2)/beta2;

basemap = \
    [ cos(mu1) + alpha1 * sin(mu1), beta1 * sin(mu1), 0, 0; \
     -gamma1 * sin(mu1), cos(mu1) - alpha1 * sin(mu1), 0, 0; \
     0, 0, cos(mu2) + alpha2 * sin(mu2), beta2 * sin(mu2); \
     0, 0, -gamma2 * sin(mu2), cos(mu2) - alpha2 * sin(mu2)]


oneturnmap = basemap*bbkick * basemap;

				# get eigenvalues

[v,l] = eig(oneturnmap)

# Extract tunes
printf("Tunes: \n");
tunes = zeros(4,1);
for i =[1:4]
  ll = l(i,i);
  tunes(i) = imag(log(ll))/(2*pi);
  printf("    %g \n", tunes(i));
endfor