Feature: shifted quadrupole support

[i-a-morozov]

Added shifted quadrupole support (shifted flag).
Shifted quadrupole is expected to be modeled as an sbend.
Quad strength or shift update recompute sbend length and angle.

Configuration example:

  quad_bend:
    regex: ^QAB
    shifted: true 
    components:
      - B2L: magnet_calibration
    dx: magnet_dx
    dy: magnet_dy
    roll: magnet_roll

[kparasch]

This is quite a violent change!

It is very important that BendingAngle is never changed because it defines the reference orbit, which in turn defines the positioning of all elements.
This is also the reason for the if-statement in L46-L63 in pySC/configuration/magnets_conf.py .

I think the proper way to model this would be to create a link between the k1 (quadrupole) component to the k0 (dipole) component of the quadrupole magnet.

The shifted quadrupole should be modeled as a dipole with at least BendingAngle and k1 != 0.
Then we can assume that:
k_0 = dx * k_1 => dx = k_0 / k_1
and calculate the "design shift" dx. We should still not actually shift the quadrupole because it is reflected in the bending angle.
There are two case:

1. for AT k_0 = element.BendingAngle / element.Length + PolynomB[0]
2. for XSuite k_0 = element.k0

or in code:

assert SC.lattice.is_dipole(index)
if type(SC.lattice) is ATLattice:
    k_0 = SC.lattice.get_bending_angle(index) / magnet_length + SC.lattice.get_magnet_component(index, "B", 1)
else:
    k_0 = SC.lattice.get_magnet_component(index, "B", 1)
k_1 = SC.lattice.get_magnet_component(index, "B", 2)

dx = k_0/k_1

next we need to create a link between the true value of k_1 and k_0 and attach it to the control of k_1:

control_name = f"{magnet.name}/{component}" # component should be B2 or B2L, depending what was previously configured 
target_name = f"{magnet.name}/B1"
quad_shift_link = ControlMagnetLink(
                link_name=f"{control_name}->{target_name}",
                magnet_name=magnet.name,
                control_name=control_name,
                component="B",
                order=1,
                is_integrated=False
            )

if to_design:
    quad_shift_link.error.factor = dx
    quad_shift_link.error.offset = 0
else: 
    # Here we need to reconstruct k1 from the setpoint according to the error.
    # The error can be found in the link called f"{control_name}->{control_name}"
    error_link = magnet_settings.links[f"{control_name}->{control_name}"]
    error_factor = error_link.error.factor
    error_offset = error_link.error.offset
    # Relation between setpoint and true_k1 is:   true_k1 = error_offset + setpoint * error_factor
    # We need factor and offset such that true_k0 = offset + setpoint * factor
    # By replacing in  true_k0 = dx  * true_k1 eq we get:
    quad_shift_link.error.factor = dx * error_factor
    quad_shift_link.error.offset = dx * error_offset

Finally we need to register the link:

    magnet_settings.add_link(quad_shift_link)

All of this can be added in pySC/configuration/magnets_conf.py, at the end of the function generate_default_magnet_control.

P.S. I will think better if k_0 = dx * k_1 is the best assumption we can make or we should do better.

[i-a-morozov]

Thanks for the feedback, I'll think about improving on this implementation, might be keeping sbend unchanged and adding a corrector at it's center (then link the corrector kick with quad str).

If I've not messed up, this implementation changes both angle and length of the sbend, i.e. links quad str to length and angle.

It is very important that BendingAngle is never changed because it defines the reference orbit, which in turn defines the positioning of all elements.

Right, but this should be done for the design lattice, after init changing should be allowed.

[kparasch]

It is very important that BendingAngle is never changed because it defines the reference orbit, which in turn defines the positioning of all elements.

Right, but this should be done for the design lattice, after init changing should be allowed.

Even after init, it should not be changed. I have made an at example to illustrate the behavior of modifying the BendingAngle:

import at
import numpy as np

r_zero = np.zeros(6).T
dip = at.Dipole('dip', length=1, bending_angle=1e-3)

print(dip.track(r_zero))
# output:
# [0. 0. 0. 0. 0. 0.]

dip.BendingAngle += 1e-6
print(dip.track(r_zero))
# output:
# [0. 0. 0. 0. 0. 0.]

dip.BendingAngle -= 1e-6
dip.PolynomB[0] += 1e-6
print(dip.track(r_zero))
# output:
# [-4.99999958e-07 -9.99999833e-07  0.00000000e+00  0.00000000e+00
#  0.00000000e+00 -1.66499992e-10]

When the BendingAngle is changed, you actually don't see a change in the particle coordinates at the exit. That's because those coordinates are defined with respect to the reference orbit which BendingAngle changes. So if you imagine that there is a BPM right after the "Dipole", then by changing the BendingAngle you also move the BPM to the center of the trajectory.
It can be slightly counter-intuitive when comparing different codes, but only PolynomB[0] should be changed so as to always keep the same reference frame.

Thanks for the feedback, I'll think about improving on this implementation, might be keeping sbend unchanged and adding a corrector at it's center (then link the corrector kick with quad str).

If I've not messed up, this implementation changes both angle and length of the sbend, i.e. links quad str to length and angle.

Adding a corrector is not necessary. Changing PolynomB[0] of the same element does the same thing without needing extra elements.

[i-a-morozov]

Adding a corrector is not necessary. Changing PolynomB[0] of the same element does the same thing without needing extra elements.

I see where the problem is now with the angle, and this solution works, and there is no need to change the arc length.

[kparasch]

We are almost there!
I understand it can be a bit frustrating to develop, thank you for the patience and for the work! I will try to write something that can serve as a developer's guide in the future.

Unfortunately, I am leaving for holiday and will be back next week on Thursday again. I may be able to review something meanwhile but I cannot promise.

[kparasch]

I think we do not need interaction with the support system any more.
Assuming ofcourse that the "design" values of the misalignments are 0, but this I think is a topic for another PR.

[kparasch]

magnet.B and magnet.A are values that are recalculated from all the "control setpoints" and the "links". So if you change it directly it will become invalid next time the setpoint is changing. Instead you need to create a link between the setpoint and magnet.B[0].

[i-a-morozov]

Same for me, I'll be away and unlikely to finalize this one until I return somewhere near the next week end, happy holidays!