% 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 2: adding the energy losses due to synchrotron radiation
%
%
% In DAFNE the average energy lost per turn due to synchrotron radiation is 8.88 keV.
% Parameter
% Value
% h
% 120
% C_r
% 97.6 m
% E_0
% 510 MeV
% V_{\textrm{rf}}
% 120 kV
% \alpha_0
% 0.018
% q
% e
% U_0
% 8.88 keV
% We set the parameters and compute again periods and frequencies as in Part 1.

% 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 = ???;
U0 = ???; % [eV]

% Compute periods and frequencies
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);

%
% Including the synchrotron radiation term, the longitudinal equations of motion become
% \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 }-\frac{U_0 }{E_0 }\left(1+2\delta^{\left(n\right)} \right)
% \end{array}

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 - ???;

%
% We have again four equispaced bunches in the ring, with filled buckets 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 all zero. Note that the synchronous phase is not anymore 90°, but
% \Delta \varphi_s =\arccos \left(\frac{U_0 }{{\textrm{eV}}_{\textrm{rf}} }\right)

n_bunches = 4;
phi_s = acos(???); % [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.

n_turns = 100000; 
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);
bunch_phases = phase_matrix(:,1)*180/pi;
bunch_deltas = delta_matrix(:,1);
axis([min(bunch_phases) max(bunch_phases) min(bunch_deltas) max(bunch_deltas)]);
xlabel('Phase [deg]');
ylabel('delta [1]');
grid
plot_every = 200;
plot_up_to_turn = 50000;
for ii = 1:plot_every:plot_up_to_turn
    pause(0.1)
    set(PLOT, 'XData', bunch_phases(ii), 'YData', bunch_deltas(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(bunch_phases, 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
hold on

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

[pks_phase,locs_phase] = findpeaks(bunch_phases);
xline(locs_phase(1), color = 'k');

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

Ts_TD = (locs_phase(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(bunch_phases, 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_freq,locs_freq] = findpeaks(abs(spectrum));
fs_FD = freq_array(locs_freq(1));
xline(fs_FD*1e-3, color = 'k');
grid
fprintf('fs_FD [kHz] = %f', fs_FD*1e-3);

%
% 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}} }\sin \Delta \varphi_s -\alpha_{r,\textrm{SR}}^2 } where \alpha_{r,\textrm{SR}} =-\frac{U_0 }{E_0 T_0 } is the damping rate of the oscillations.

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

%
% We make an exponential fit of the envelope of the phase oscillations. Check that the obtained damping rate is close to the analytical one.

% isolate the maxima
time_array = T0:T0:n_turns*T0;
time_peaks = time_array(locs_phase);
phase_peaks = bunch_phases(locs_phase);
% compute the initial guesses
c0 = mean(???);
b0 = log((phase_peaks(1)-c0)/(phase_peaks(2)-c0))/(time_peaks(1)-time_peaks(2));
a0 = (phase_peaks(1)-c0)*exp(-b0*time_peaks(1));
% compute the fit
g = fittype('a*exp(b*x)+c','indep','x');
fit_object = fit(time_peaks', phase_peaks, g, 'StartPoint',[a0 b0 c0]);
coefficients = coeffvalues(fit_object);
fprintf('alpha_damp_FIT = %f\n', coefficients(2))
fprintf('alpha_damp_AN = %f\n', alpha_damp_AN)
% plot of the fit
fig4 = figure();
plot(time_array,bunch_phases*180/pi, color = 'b')
hold on
phase_fit = coefficients(1)*exp(coefficients(2)*time_array)+coefficients(3);
plot(time_array,phase_fit*180/pi, 'g', 'LineWidth',5)
xlabel('Time [s]');
ylabel('Phase [deg]');
grid

%