Add ProportionalIntegral Loop and reorg for cavity feedbacks

[lcarver]

Until now, only a proportional loop was possible for the RF feedbacks. This was fine in the single harmonic case, but for harmonic cavities where the beam is sensitive to the harmonic voltage, it was not possible to avoid driving large voltage oscillations.

This PR performs a direct implementation from the nice work done by Naoto Yamamoto, where he derives a ProportionalIntegral loop, with the basis of generator current. By using generator current, he also considers the cavity response in the feedback. His work was implemented into Mbtrack2 ( https://gitlab.synchrotron-soleil.fr/PA/collective-effects/mbtrack2/-/blob/develop/mbtrack2/tracking/rf.py?ref_type=heads )

I have done a complete integration of his MBTRACK2 implementation into AT. This version is now working with an example script shown at the end of this post.

Some things still to do:

1. Test all previous beam loading functionalities (passive cavity tuning, buffer averaging).
2. Fix the shapes of the long complex arrays when displayed in python.
3. Try to improve the clarity and organisation of the code.
4. Test and benchmarking or more complex cases including with harmonic cavities.
5. Write the PyTests

For now it is a WIP, I am fully open to any suggestions. I worked hard to try to be efficient from a memory point of view. There are many buffers that needed to be initialised in the beginning. If you think some are unneeded, I am happy to hear suggestions.

import numpy as np
from at.constants import clight
import at
import matplotlib.pyplot as plt

ring_dict = {
            'circumference': 843.977345496621,
            'harmonic_number': 992,
            'ac': 8.511728693902269e-05,
            'energy': 6e9,
            'sigma_e': 0.0009356632506279087,
            'U0': 2532505.857673822
            }

Vc = 5.5e6
Q0 = 37500
beta = 2.8
QL = Q0/(1+beta)
Rsl = 145*11*QL
   
phi_s = np.pi - np.arcsin((ring_dict['U0']/Vc))


Nbunches = 992
I0 = 50e-3
Nslice = 1#00

nparts_per_bunch = 1
Nparts = Nbunches

Nturns = 5000
sigma_matrix = at.sigma_matrix(betax=1, betay=1, alphax=0, alphay=0, emitx=100e-12, emity=10e-12, espread=ring_dict['sigma_e'], blength=11e-3)



simple_ring = at.simple_ring(ring_dict['energy'], ring_dict['circumference'], ring_dict['harmonic_number'], 0.1, 0.1, Vc, ring_dict['ac'], U0=ring_dict['U0'], tauz=3229.114299199956, espread=ring_dict['sigma_e'])
simple_ring.set_fillpattern(Nbunches)

simple_ring.set_beam_current(I0)

#Setup main cavity
mc_elem = simple_ring[0]
tl_main = clight*(np.pi - phi_s)/(2*np.pi*mc_elem.Frequency)
mc_elem.TimeLag = tl_main

# add beamloading
#gain = [0.0,0.0]
tunergain = 1e-4
gain = [0.1,1e4]

at.add_beamloading(simple_ring, QL, Rsl, cavpts=[0], Nslice=1, Gain=gain, TunerGain=tunergain, cavitymode=at.CavityMode.ACTIVE, feedback_angle_offset=0.0, fbmode=at.FeedbackMode.PROP_INTEGRAL, sample_num=1, every=992, delay=1)

parts = at.beam(Nparts, sigma_matrix)

all_Vg = np.zeros((4, Nturns))
allI = np.zeros(Nturns,dtype=complex)
for iturn in np.arange(Nturns):

    _ = simple_ring.track(parts, in_place=True, refpts=None, nturns=1, keep_counter=True)
    all_Vg[:,iturn] = simple_ring[0].Vgen
    print('I_record', simple_ring[0]._I_record)
    allI[iturn] = simple_ring[0]._I_record[0] + 1j*simple_ring[0]._I_record[1]    
    
fig, ax = plt.subplots(4,1)
ax[0].plot(all_Vg[0,:])
ax[1].plot(all_Vg[1,:])
#ax[2].plot(all_Vg[2,:])
ax[2].plot(np.real(allI))
ax[3].plot(np.imag(allI))
plt.show()