% Author: D. Quartullo
% The code below has been readapted from the Fortran code developed by M. Migliorati in 1993 for his PhD studies on the DAFNE accelerator.
% Some references:
% Migliorati M. "Simulazioni numeriche della dinamica multipacchetto", PhD thesis, University of Rome 'La Sapienza', Rome, 1996.
% Bassetti M. et al. "A time domain simulation code of the longitudinal multibunch instabilities", DAFNE technical note G-19, INFN-LNF Accelerator Division, Frascati, 1993.
% Gallo A. et al. "Simulations of the bunch-by-bunch feedback operation with a broadband RF cavity as longitudinal kicker", DAFNE technical note G-31, INFN-LNF Accelerator Division, Frascati, 1995.
% There are four scripts in total, named Parts 1, 2, 3 and 4. Each Part is indipendent of the others, and no external libraries are needed to run it. Nevertheless, it is advisable to study the four Parts in order, since each builds on top of the previous ones.
% The characters ??? indicate code portions that the reader has to correctly fill in to be able to procede with the script. It is advisible to run the script section after section.
%% Part 1: basic synchrotron motion in DAFNE
%
%
%% DAFNE main parameters
%
% Parameter
% Value
% h
% 120
% C_r
% 97.6 m
% E_0
% 510 MeV
% V_{\textrm{rf}}
% 120 kV
% \alpha_0
% 0.018
% q
% e

% Clear variables and figures
clc
clf
clear all
close all force

% Insert the parameters of the table
h = ???;
Cr = ???; % [m]
E0 = ???; % [eV]
Vrf = ???; % [V]
alpha0 = ???;

%
% Supposing \beta =1, compute the revolution period and the RF period.

c = physconst('LightSpeed'); % [m/s]
T0 = ???; % [s]
Trf = ???; % [s]
f0 = ???; % [Hz]
frf = ???; % [Hz]
fprintf('T0 [ns] = %f', T0*1e9);
fprintf('Trf [ns] = %f', Trf*1e9);
fprintf('f0 [MHz] = %f', f0*1e-6);
fprintf('frf [MHz] = %f', frf*1e-6);

%
% In the code, each bunch is represented by just one macroparticle. The longitudinal equations of motion for each bunch are
% \begin{array}{l}
% {\Delta \varphi }^{\left(n+1\right)} ={\Delta \varphi }^{\left(n\right)} +2\pi h\alpha_0 \delta^{\left(n\right)} \\
% \delta^{\left(n+1\right)} =\delta^{\left(n\right)} +\frac{eV_{\textrm{rf}} \;\cos {\Delta \varphi }^{\left(n+1\right)} }{E_0 }
% \end{array}
% the synchronous phases being in \frac{\pi }{2}+2k\pi \;\left(k=1,2,\ldotp \ldotp \ldotp \right). We can assume that the synchronous phase is \frac{\pi }{2} for each bunch, since the RF voltage is 2\pi-perdiodic.
% The functions to update the phase and energy coordinates from turn n to turn n+1 are

phase_update = @(phase_array,delta_array) phase_array + 2*pi*h*alpha0*delta_array;
delta_update = @(phase_array,delta_array) delta_array + Vrf*cos(phase_array)/E0;

%
% Let's suppose that we have four equispaced bunches in the ring, so the filled buckets are 1, 31, 61, 91.
% As for the initial conditions, the phase of each bunch is the synchronous phase plus a deviation of 1 deg. The energy deviations are zero.

n_bunches = 4;
phi_s = pi/2; % [rad]
init_phase = phi_s + ???;
phase_array = init_phase*ones(1,n_bunches); 
delta_array = zeros(1,n_bunches);

%
% Let's track the four bunches and save the coordinates in a matrix at each turn.

n_turns = 50000; 
phase_matrix = zeros(n_turns, n_bunches);
delta_matrix = zeros(n_turns, n_bunches);
for i = 1:n_turns
    phase_array = phase_update(phase_array,delta_array);
    delta_array = delta_update(phase_array,delta_array);
    phase_matrix(i,:) = phase_array;
    delta_matrix(i,:) = delta_array;
end

%
% Let's plot the coordinates of bunch 1 in phase space.

fig1 = figure();
plot(phi_s*180/pi,0,'.', color = 'r', MarkerSize = 30);
hold on
PLOT = plot(NaN,NaN);
DATASET1 = phase_matrix(:,1)*180/pi;
DATASET2 = delta_matrix(:,1);
axis([min(DATASET1) max(DATASET1) min(DATASET2) max(DATASET2)]);
xlabel('Phase [deg]');
ylabel('delta [1]');
grid
plot_every = 10;
for ii = 1:plot_every:1000
    pause(0.1)
    set(PLOT, 'XData', DATASET1(ii), 'YData', DATASET2(ii), ...
        marker ='.', color = 'b', MarkerSize = 30);
    title(strcat('Turn =  ', num2str(ii)))
end
% Create the same type of phase-space animation for bunch 2 and compare the results.
???

%
% We now plot the phase and energy oscillations as a function of turn.

fig2 = figure();
yyaxis left
plot(phase_matrix(:,1)*180/pi, color = 'b');
xlabel('Turn');
ylabel('Phase [deg]');
set(gca,'YColor','b');
yyaxis right
plot(delta_matrix(:,1), color = 'r');
ylabel('delta [1]');
set(gca,'YColor','r');
grid on
xlim([0 500]);
hold on

% Using findpeaks and plotting a vertical line, check that the phase and energy curves are precisely misaligned by 90°.

[pks,locs] = findpeaks(phase_matrix(:,1));
xline(locs(1), color = 'k');

%
% Estimate the synchrotron frequency by using the locations of two consecutive maxima of the phase curve.

Ts_TD = (locs(2) - ???)*T0;
fs_TD = ???;
fprintf('fs_TD [kHz] = %f', fs_TD*1e-3);

%
% Make an FFT of the phase curve and estimate the synchrotron frequency using the frequency domain.

fig3 = figure();
freq_array = f0/n_turns*(0:n_turns-1);
spectrum = fft(phase_matrix(:,1), n_turns);
plot(freq_array*1e-3, abs(spectrum), color = 'r');
xlabel('Frequency [kHz]');
ylabel('Abs(spectrum) [rad]');
xlim([0 f0/2*1e-3]);
set(gca, 'YScale', 'log')
[pks,locs] = findpeaks(abs(spectrum));
fs_FD = freq_array(locs(1));
xline(fs_FD*1e-3, color = 'k');
fprintf('fs_FD [kHz] = %f', fs_FD*1e-3);
grid

% Check that these synchrotron frequency values are close to the synchrotron frequency obtained with the analytical formula
% f_s =\sqrt{\frac{{\textrm{eV}}_{\textrm{rf}} \;\alpha_0 }{2\pi E_0 T_0 T_{\textrm{rf}} }}

fs_AN = sqrt(???*Vrf*alpha0/(2*pi*E0*T0*Trf));
fprintf('fs_AN [kHz] = %f', fs_AN*1e-3);