Simple Ring modification to solve energy loss computation bug

[lcarver]

This is proposed as a solution to #694 .

Old simple ring contained: ring = [RF Cavity, M66 no rad, Nonlinear element, SimpleQuantDiff] where SimpleQuantDiffhandled all radiation effects (damping, diffusion and U0) This setup was problematic due to the fact that get_energy_loss could not compute the U0 without the damping and diffusion, which had some other knock-on problems.

Now simple ring contains: ring = [RF Cavity, M66 no rad, SimpleRadiation, Nonlinear element, SimpleQuantDiff] where SimpleRadiation is a new PassMethod that only performs U0 and radiation damping and SimpleQuantDiff has been modified to only act as the diffusion source for emittance.

One possible problem that I want to draw to your attention, is that now we have 2 diffusion elements (QuantumDiffusion and SimpleQuantDiff) and I wonder if it is now more convenient to create a Diffusion class/group (not sure of the correct noun) and then each element has the property is_diffusion to enable cleaner code (for example you can see in get_energy_loss I have to mention both elements specifically.

If this is merged, then the structure can become the basis for the fast_ring reorganisation, which will be done in a separate PR.

[lfarv]

Hello @lcarver and @swhite2401

The idea is nice but I think that the implementation is wrong. More exactly, SimpleRadiation is wrong, its more complicated that that.

In the general 6D case; the 6x6 one-turn transfer matrix M_rad can be expressed using the A matrix as:

M_rad = A . R . Λ . A^-1

with R block diagonal rotation matrix (phase advance) and Λ diagonal damping matrix.

so we have:

M_rad = (A . R . A^-1) . (A . Λ . A^-1)

M_rad = M_norad . DAMP

The damping matrix DAMP is A . Λ . A^-1 !

Λ is similar to what you wrote in SimpleRadiation, but:

1. The horizontal damping applies to both x and x' (and same for the other planes),
2. It has to be applied on the "normalised" coordinates. That's what the A matrice does. In the uncoupled case, A is the classical [ [1/sqrt(β) 0][ α/sqrt(β) sqrt(β)]] matrix, but in the end, the SimpleRadiation is a 6x6 block diagonal matrix.

In conclusion, this option (once corrected) costs an additional matrix multiplication.

[swhite2401]

@lfarv, this algorithm is implemented in many well know codes, is well benchmarked and gives the correct result in term of damping time and equilibrium emittance. This is exactly what we are looking in multi-particle simulation where the full ring is represented by a matrix.

We could go for your solution but then the algorithm becomes more complicated since you have to derive the 6x6 damping matrix, providing the huge approximation we are already making in modeling the ring I do not see much interest in this.

[lfarv]

@swhite2401 : I have no doubt that you get the correct damping times and emittances, but for sure the transverse motion is wrong. The damping is applied as if the betas were 1.0. You maybe don't care about transverse optics: you can run multi particle tracking were the turn is defined by simply a phase advance, no need for beta functions, the one turn map is a simple rotation, and then your formula is correct. You just define tunes, chromaticities, and that's it. But if for some reason you want to introduce betas and alphas, it's different.

I'm 100% sure the formulas above are correct, so anything else is wrong. But I understand the need for speed, so I had a look numerically: in my notations,

Λ = diag([exp(-1/τ_x), exp(-1/τ_x), exp(-1/τ_y), ...]) with τ in turns. Λ

For a sample ring (hmba.mat), the damping matrix DAMP= A . Λ . A^-1 is still close to the diagonal matrix Λ. The main terms are dispersive ones, their order of magnitude is about 10^-8. So that's OK for the approximation. Ignoring them,

DAMP ≈ diag([exp(-1/τ_x), exp(-1/τ_x), exp(-1/τ_y), ...])

DAMP ≈ diag([(1-1/τ_x), (1-1/τ_x), (1-1/τ_y), (1-1/τ_y), (1-1/τ_z), (1-1/τ_z)]) for τ >> 1

Compared to what is in SimpleRadiationPass,

1. multiplying by the exact exponential (computed once and stored) would be faster (r6[1] *= xdamp compared to r6[1] -= r6[1]/taux), with xdamp = exp(-1.0/taux)
2. the damping should be applied to both x and p_x instead of doubling the effect on p_x. Why is it done like that?

[lfarv]

Coming back to more practical problems: I agree that a Diffusion class would help. But at the same time, the SimpleRadiation element should also be handled automatically in enable_6d, without explicitly specifying it. So it should be declared as

class SimpleRadiation(_DictLongtMotion, Radiative, Element):

in that order, so that _DictLongtMotion provides the right methods, and Radiative is just used for grouping. That's not very elegant, and there could be some reorganisation of these mixin classes, to separate functionality and grouping. But that's another topic…

[lcarver]

I modified SimpleRadiation in the following way:

Added Radiative to element definition. Passed beta, alpha, dispand disp' to the SimpleRadiationPass.c In the python, use taux,y,z to construct 6d damping matrix.

Assuming the following formula to move from real to normalised coordinates:

X = x/sqrt(beta) - Dx*dpp
XP = xp*sqrt(beta) + x * alpha / sqrt(beta) - DPx*dpp

The particles are converted, damping applied, then reverted. U0 is applied at the end after all changes.

[lcarver]

Below is a test file for your convenience:

import at
from at.physics.fastring import simple_ring
import numpy as np
import matplotlib.pyplot as plt
from scipy.signal import hilbert
from at.constants import clight

dic = {
       'alphax': 1e-3,
       'alphay': 1e-3,
       'betax': 20,
       'betay': 10,
       'dispx': 0.1,
       'dispxp': 1e-3,#0.0,
       'dispy': 0.1,
       'dispyp': 1e-3,#,
       'Qx': 74.1235,
       'Qy': 28.3,
       'Qpx': 0,
       'Qpy': 0,
       'A1': 0.,
       'A2': 0.,
       'A3': 0.,
       'circumference': 843.977,
       'alpha': 8.5e-5,
       'U0': 2.5e6,
       'Vrf': 6e6,
       'energy': 6e9,
       'harmonic_number': 992,
       'emitx': 100e-12,
       'emity': 10e-12,
       'espread': 9.5e-4,
       'taux': 3000,
       'tauy': 5000,
       'tauz': 1000
}

ring = simple_ring(dic['energy'], dic['circumference'], dic['harmonic_number'], dic['Qx'], dic['Qy'], dic['Vrf'], dic['alpha'], taux=dic['taux'], emitx=dic['emitx'], tauy=dic['tauy'], emity=dic['emity'], tauz=dic['tauz'], espread=dic['espread'], U0=dic['U0'], A1=dic['A1'], A2=dic['A2'], A3=dic['A3'], Qpx=dic['Qpx'], Qpy=dic['Qpy'], alphax=dic['alphax'], alphay=dic['alphay'], betax=dic['betax'], betay=dic['betay'], dispx=dic['dispx'], dispy=dic['dispy'], dispxp=dic['dispxp'], dispyp=dic['dispyp'])

r_in = np.zeros((6,1))
r_in[0] = 1e-1
r_in[2] = 1e-1
r_in[4] = 1e-3

# Check radiation damping
nt = 10000
r_out = at.lattice_pass(ring, r_in, nturns=nt)

xdat = r_out[0,0,0,:]
fitx = np.polyfit(np.arange(2500), np.log(np.abs(hilbert(xdat)))[500:3000], 1)
print('taux', -1/fitx[0], dic['taux'])

ydat = r_out[2,0,0,:]
fity = np.polyfit(np.arange(2500), np.log(np.abs(hilbert(ydat)))[500:3000], 1)
print('tauy', -1/fity[0], dic['tauy'])

#spikes at start and finish are artefacts of hilbert transform
plt.suptitle('envelope from hilbert transform')
plt.semilogy(np.abs(hilbert(xdat)),'r', label='x')
plt.semilogy(np.abs(hilbert(ydat)),'b', label='y')
plt.xlabel('turn')
plt.ylabel('env [m]')
plt.legend()
plt.show()

r_in = np.zeros((6,1))
r_in[0] = 0.
r_in[2] = 0.
r_in[4] = 1e-3
r_out = at.lattice_pass(ring, r_in, nturns=nt)
r_in[4] += 5e-2
r_out = at.lattice_pass(ring, r_in, nturns=nt)
dpdat = r_out[4,0,0,:]
zdat = r_out[5,0,0,:]

fitz = np.polyfit(np.arange(500), np.log(np.abs(hilbert(dpdat)))[250:750], 1)
print('tauz', -1/fitz[0], dic['tauz'])

plt.suptitle('abs(hilbert(dp))')
plt.semilogy(np.abs(hilbert(dpdat)))
plt.show()

[lcarver]

Ready for re-review.

[lcarver]

On my side, for betax and betay provided, alphax=alphay=0 and dispersion 0. Damping times are perfect. For an enormous alphax=alphay=1, I see perfect damping in V and 2997 turns (compared to 3000) for H.

For alphax=alphay=0 and Dx=Dy=0.1, I see a deviation in H. Hdamping = 3215 turns (instead of 3000) Vdamping = 5031 turns (instead of 5000).

So I think there is a mistake in how I implement the dispersion.

[lcarver]

Fixed dispersion problem. Should be ok now. Ready for review.

[lcarver]

all changes implemented, but one outstanding question to be resolved (I made a comment in the code)

[lfarv]

@lcarver : I just found out that the linear matrix is missing the dispersion. One must add:

    M04 = (1.0 - M00) * dispx - M01 * dispxp
    M14 = -M10 * dispx + (1.0 - M11) * dispxp

    M24 = (1.0 - M22) * dispy - M23 * dispyp
    M34 = -M32 * dispy + (1.0 - M33) * dispyp

[lcarver]

@lcarver : I just found out that the linear matrix is missing the dispersion. One must add:

    M04 = (1.0 - M00) * dispx - M01 * dispxp
    M14 = -M10 * dispx + (1.0 - M11) * dispxp

    M24 = (1.0 - M22) * dispy - M23 * dispyp
    M34 = -M32 * dispy + (1.0 - M33) * dispyp

Well spotted. I added your modifications to the M66.

[lfarv]

...but one outstanding question to be resolved (I made a comment in the code)

If this was a test, I admit that I failed… What is this question ?