Quad curvature correction in bending magnets

[lfarv]

This adds another term to the bending magnets to try and cancel a difference of tracking of combined bending-focusing magnets between AT and Xsuite.

The so-called "K1h" term of the Hamiltonian is included in:

• BndMPoleSymplectic4Pass, BndMPoleSymplectic4RadPass,
• ExactSectorBendPass, ExactSectorBendRadPass.

A successful comparison is necessary before accepting this term. @swhite2401 and @kandre2, can you look at that ?

[lfarv]

Some tests are failing because of small deviations introduced by the new tracking. Please ignore that for the moment, I will adapt the tests if the modification proves successful.

[lfarv]

A few other points appear when looking into that:

• the tests should be done with a small radius of curvature (1 or 2m) to enhance the effect of the K1h term. And preferably a significant bending angle,
• there is another term K0h called "weak focusing", for which I am not sure whether it's included in the curvature of the reference trajectory (which implies that K0 must be removed from the field expansion), or if it should be added. Depending on the test results, we'll see.
• The Xsuite treatment of the exact bend (in the "bend-kick-bend" model) is completely different from what is done in AT, based on E. Forest's formula.
• In the "Xsuite physics manual", there is still another term K2h appearing in equation (1.74) which is not applied neither in Xsuite nor in AT. Should we look at it?
• there is an error in the documentation for the model translation between AT and Xsuite: the Xsuite model similar to the default BndMPoleSymplectic4Pass is "drift-kick-drift-expanded" and not "rot-kick-rot". The code is correct, only the documentation is wrong.

[swhite2401]

I ran my scripts again, and I see several issues:
1- even in the absence of K1 the SectorBend is now different from Xsuite while it was identical before
2-In the presence of K1I still do not get agreement, I have look at this in more details, I suspect an issue with QuadFringeFields but only the Exact integrator has the same model as XSuite and it seems to have a problem so it is difficult to check

I will take a step back and verify that in AT the new integrator give identical results as the old one in the absence of K1.

[swhite2401]

I confirm a problem with the new sector bend, with the following code all particles are lost which is nit the case in the previous implementation

import at
import numpy as np
import pickle

nsteps = 10000
energy = 10e9

sbend_exp = at.Lattice([at.Dipole('D', 1.0, 1.0, k=0.0, PassMethod='BndMPoleSymplectic4Pass',
                                  NumIntSteps=nsteps)], energy = energy)
sbend = at.Lattice([at.Dipole('D', 1.0, 1.0, k=0.0, PassMethod='ExactSectorBendPass', 
                              NumIntSteps=nsteps)], energy = energy)
rbend = at.Lattice([at.Dipole('D', 1.0, 1.0, k=0.0, PassMethod='ExactRectangularBendPass', 
                              NumIntSteps=nsteps)], energy = energy)


n = 11
dp = 0.0
na = np.linspace(-0.1, 0.1, n)

parts = np.zeros((6, n*n))
parts[4, :] += dp

for i in range(n):
    for j in range(n):
        parts[0, i*n+j] = na[i]
        parts[1, i*n+j] = na[j]
        
pout_sbend_exp, *_ = sbend_exp.track(parts.copy())
pout_sbend_exp = np.squeeze(pout_sbend_exp)

pout_sbend, *_ = sbend.track(parts.copy())
pout_sbend = np.squeeze(pout_sbend)

pout_rbend, *_ = rbend.track(parts.copy())
pout_rbend = np.squeeze(pout_rbend)

with open('bend.pkl', 'wb') as fin:
    pickle.dump({'na': na,
                 'sbend_exp': pout_sbend_exp,
                 'sbend': pout_sbend,
                 'rbend': pout_rbend,
                 },
                fin)

[lfarv]

Very strange: the added term has B[1] in factor. Is the problem specific to ExactSectorBendPass or is BndSymplectic4Pass also affected ?

Note that ExactRectangularBendPass is not modified: it's a straight magnet.

I will use you script to further check.

[swhite2401]

Very strange: the added term has B[1] in factor. Is the problem specific to ExactSectorBendPass or is BndSymplectic4Pass also affected ?

Note that ExactRectangularBendPass is not modified: it's a straight magnet.

I will use you script to further check.

BndSymplectic4Pass seems to be fine at first sight

[lfarv]

I see problems with ExactSectorBendPass even on the master branch… I have to look more closely.

[swhite2401]

I see problems with ExactSectorBendPass even on the master branch… I have to look more closely.

Do you? It looked fine for me providing the FringeFields are correctly set (full, or dipole-only in XSuite).

[lfarv]

A bug in ExactSectorBendPass is fixed, and I now get exactly the same result as in the master branch for K == 0.

Note: if K == 0, ExactSectorBendPass results as independent of NumIntSteps. You can even use NumIntSteps = 0, special case which makes no kick at all: a single "bend". You can check with your script that the results are exactly identical, with a huge speed increase.

Now for K != 0, it's up to you !

[swhite2401]

Now for K != 0, it's up to you !

Great thanks! May have to wait until next week (after the restart) before getting back it it though

[swhite2401]

I confirm that now in the absence of K1, the master and this banch give the same results.
However, I still see some substantial differences between AT and XSuite for Dipole-Quadrupole configuration, and this for both SBEND and RBEND.

Anything I can try?

[swhite2401]

Without K1, there is some "pattern" in the RBEND but agreement is within 1.0-13

[swhite2401]

With K1=1, still some discrepancies....

[swhite2401]

This is what I had on the master.... identical in the H plane and better in the V plane...

[swhite2401]

Other problem, in ExactRectangular bend the case where FringeBend*=1 and QuandBend*=0 now translates to 'linear' instead of 'dipole-only'.

[swhite2401]

So it is strange now... some things that used to work don't anymore.... are you sure X0ref is correctly handle? Don't you need to change rbendtune() as well?

[lfarv]

You introduced a bug in xsuite.py, the keyword for fringe in xsuite is "suppressed" and not "suppress"

corrected

in ExactRectangular bend the case where FringeBend*=1 and QuandBend*=0 now translates to 'linear' instead of 'dipole-only'.

Yes, here is what is happening at the moment: according to the doc:

    /*     method 0 no fringe field
     *     method 1 legacy version Brown First Order
     *     method 2 SOLEIL close to second order of Brown
     *     method 3 THOMX
     */

So 1 should indeed translate to "linear". But at the moment, the AT exact methods ignore the value of FringeBend* and always use another method (Forest), even when specifying 0. In the future, this will be method 4, and 1 will be indeed the "linear" method (this is how I tested the agreement with Xsuite). For the time being, don't specify FringeBend* on exact methods: by default it is translated to "dipole-only".

I'm presently working to implement all 4 methods on all the dipole passmethods.

[lfarv]

are you sure X0ref is correctly handle? Don't you need to change rbendtune() as well?

rbendtune just iterates over X0ref until the reference particle exits at 0 angle, whatever x0ref does. So no need to change it. The resulting value may be different from what it was before, but it's correct.

[swhite2401]

are you sure X0ref is correctly handle? Don't you need to change rbendtune() as well?

rbendtune just iterates over X0ref until the reference particle exits at 0 angle, whatever x0ref does. So no need to change it. The resulting value may be different from what it was before, but it's correct.

Still I still get strange results:

nsteps = 10000
energy = 10e9

at_sbend = at.Lattice([at.Dipole('D', 1.0, 1.0, k=0.0, PassMethod='ExactRectangularBendPass', 
                              NumIntSteps=nsteps, EntranceAngle=0.5, ExitAngle=0.5,
                              PolynomB=[0.0, 0.0], PolynomA=[0.0, 0.0], dx = 0.0)],
                              energy = energy)

fbend = 1
fquad = 1
at_sbend[0].FringeBendEntrance = fbend
at_sbend[0].FringeBendExit = fbend
at_sbend[0].FringeQuadEntrance = fquad
at_sbend[0].FringeQuadExit = fquad

p0, *_ = at_sbend.track(np.zeros(6))
print(np.squeeze(p0))
at_sbend[0].rbendtune()
p0, *_ = at_sbend.track(np.zeros(6))
print(np.squeeze(p0))

For me the reference particle is np.zeros(6) before rbendtune() it already exists at zero, so why is rbendtune() moving it to some other coordinates? ... maybe I miss someting?

In any case not running rbendtune in this case I get perfect agreement using:

xs_sbend = at.line_from_lattice(at_sbend, match_model=True)
xs_sbend[0].rbend_shift += xs_sbend[0]._x0_in  - getattr(at_sbend[0], 'X0ref', 0.0)

I would have think that running rbendtune and shifting the xsuite element by the same amount... so something is strange....and clearly the problem lies there

[lfarv]

@swhite: you magnet has K1 == 0. So x0Ref has no effect: you can shift a uniform magnetic field by any amount, the result is the same. if K1 == 0, the result should be the same with or without running rbendtune.

Similarly, if K1 == 0, the Xsuite result should be independent of rbend_shift. If not, there is a problem there !

[lfarv]

@swhite2401: I checked that before iterating, rbendtune shifts the magnet by half the sagitta, so that the trajectory spans equally on each side of the magnetic axis. This gives a good approximation. If K1 == 0, it stops there.

[swhite2401]

@swhite: you magnet has K1 == 0. So x0Ref has no effect: you can shift a uniform magnetic field by any amount, the result is the same. if K1 == 0, the result should be the same with or without running rbendtune.

Similarly, if K1 == 0, the Xsuite result should be independent of rbend_shift. If not, there is a problem there !

I agree, I did not expect this! This is why I am surprised, but maybe I am doing something wrong.... let me check and get back to you...

[lfarv]

This new release makes 3 changes:

1. The point where the x0ref is applied is moved: after the axis rotation and before the fringe field. With this, the agreement with Xsuite for a rectangular magnet with K1 != 0 is much improved,
2. I suppressed the quadrupole wedge when FringeQuad* is 0: this allows a direct comparison with Xsuite for the "dipole-only" model. We can restore it later when everything is understood,
3. I rewrote the "exact bend fringe", bringing a small improvement in the comparison

[lfarv]

A comment on dx processing: a dipole with parallel faces (face_angle == 0.5 * bending_angle) has no effect in the horizontal plane: it's a drift. So its horizontal displacement has no effect. This is true in Xsuite, but not in AT, whatever the passmethod (BndMPoleSymplectic4Pass, ExactSectorBendPass, ExactRectangularBendPass).

This looks like a problem in transformations.

dx works perfectly on straight magnets !

[swhite2401]

@lfarv I am having a look back and here is a summary of the present situation:

Multipole: all looks good

SBEND is ok for any value of FringeQuad*, face angles or K1, dy as long as dx, B0 and A0 are equal to zero
-with dx!=0 I get 1.0e-7 disagreement (instead of 1.,0e-14) that can degrade to 1.0e-6 depending on the value of K1 or the face angles
-with A0!=0 and FringeQuad*!=1 I get 1.0e-7 disagreement that can go up to 1.0e-5 depending of face angles and K1. For FringeQuad*=0 I get good agreement in all cases
-with B0!=0 strong disagreement in all cases

RBEND
-with K1=0 everything looks fine
-with K1!=0 and FringeQuad*=0 good agreement is obtained if rbendtune() is not called
-with K1!=0 and FringeQuad*=1 strong disagreement except for the case where face angles=0.5*bendingangle
-same observation as SBEND for dx, dy, A0, B0

Few questions / remarks come:
-Did you look into the possibility of integrating B0 into the fringe field calculation as it is done in XSuite?
-There is a problem with dx as you pointed out
-X0ref is still problematic for face angle != 0.5*bendingangle

[lfarv]

@swhite2401

RBEND
-with K1=0 everything looks fine

I don't see that: even in that case, I see a strong disagreement as soon as the entrance/exit angles are not half the bending angle

[lfarv]

@swhite2401 : after another review, my conclusions on RBEND are different from yours:

With faceangles == 0.5*bendingangle

Everything works, including

• K1 != 0,
• rbendtune() and the corresponding rbend_shift,
• dipole fringe field with non-zero fringe field integrals,
• multipole fringe field.

Here is my test dipole:

elem = at.Dipole('D', 2.0, 1.0, k=1.0, PassMethod='ExactRectangularBendPass',
                  EntranceAngle=0.5, ExitAngle=0.5,
                  FullGap=0.1, FringeInt1=0.65, FringeInt2=0.65,
                  FringeQuadEntrance=1, FringeQuadExit=1,
                  NumIntSteps=1000)

Here is the translated Xsuite element:

View of RBend(length_straight=1.92, angle=1, length=2, k0='from_h', k1=1, h=0.5, k0_from_h=True, model='drift-kick-drift-exact', knl=[0., 0.], ksl=[0., 0.], knl_rel=[0.], ksl_rel=[0.], edge_entry_active=np.int64(1), edge_exit_active=np.int64(1), edge_entry_model='full', edge_exit_model='full', edge_entry_angle=0, edge_exit_angle=0, edge_entry_angle_fdown=0, edge_exit_angle_fdown=0, edge_entry_fint=0.65, edge_exit_fint=0.65, edge_entry_hgap=0.05, edge_exit_hgap=0.05, shift_x=0, shift_y=0, rot_s_rad=0)

And the result:

With faceangles != 0.5*bendingangle

nothing works

[swhite2401]

Strange let me go case by case in details then.

First here is the script I use:

nsteps = 10000
energy = 10e9

at_sbend = at.Lattice([at.Dipole('D', 1.0, 1.0, k=0.0, PassMethod='ExactRectangularBendPass', 
                              NumIntSteps=nsteps, EntranceAngle=0.5, ExitAngle=0.5,
                              PolynomB=[0.0, 0.0], PolynomA=[0.0, 0.0], dx = 0.0)],
                              energy = energy)

#fquad = 1
at_sbend[0].FringeQuadEntrance = fquad
at_sbend[0].FringeQuadExit = fquad
at_sbend[0].EntranceAngle = 0.0
at_sbend[0].ExitAngle = 0.0
at_sbend[0].K = 0.0
at_sbend[0].dx = 0.0
at_sbend[0].dy = 0.0
at_sbend[0].PolynomB[0] = 0.0
at_sbend[0].PolynomA[0] = 0.0

#at_sbend[0].rbendtune()

xs_sbend = at.line_from_lattice(at_sbend, match_model=True)
xs_sbend[0].rbend_shift += xs_sbend[0]._x0_in - getattr(at_sbend[0], 'X0ref', 0.0)

nax, _, _, diffx = at_xs_track(at_sbend, xs_sbend, axes= [0, 1])
nay, _, _, diffy = at_xs_track(at_sbend, xs_sbend, axes= [2, 3])

plt.subplot(211)
plt.contourf(nax, nax, diffx)
plt.xlabel('x [m]')
plt.ylabel('px')
plt.colorbar()

plt.subplot(212)
plt.contourf(nay, nay, diffy)
plt.xlabel('y [m]')
plt.ylabel('py')
plt.colorbar()

plt.show()

[swhite2401]

With K1=0 and rbendtune() not called: work for all face angle and model

With K1=0 and rbendtune() called:

• FringeQuad = 0, face angles = 0.0: not ok
• FringeQuad = 0, face angles = 0.5: ok
• FringeQuad = 1, face angles = 0.0: not ok
• FringeQuad = 1, face angles = 0.5: ok

With K1=1 and rbendtune() not called :

• FringeQuad = 0, face angles = 0.0: ok
• FringeQuad = 0, face angles = 0.5: ok
• FringeQuad = 1, face angles = 0.0: not ok
• FringeQuad = 1, face angles = 0.5: ok

With K1=1 and rbendtune() called :

• FringeQuad = 0, face angles = 0.0: not ok
• FringeQuad = 0, face angles = 0.5: ok
• FringeQuad = 1, face angles = 0.0: not ok
• FringeQuad = 1, face angles = 0.5: ok

For me there is still a problem with X0ref the agreement should not depend on whether we call rbendtune or not as long as rbend_shift is set correctly (am I doing something wrong?)

The the special case with K1!=0, FringeQuad=1 and face angle !=0.5 never work, so there is a problem there too unrelated to X0ref

[lfarv]

@swhite2401:
For both our tests, there is at least a common conclusion: everything is ok when the face angles are 0.5.

If not, my tests always fail, while you get agreement for your 1st case: K1 == 0 and no rbend_tune(). The only difference is that even in that case, you tune rbend_shift while I let it to the default value. Apparently this makes the difference.

So I wonder if the problem is not just a question of transverse positioning : one should always check that the reference particle comes out at 0.0 angle. If not, x0Ref or rbend_shift is wrong. Could it be that when face_angles == 0, your tuning function rbend_shift += xs_sbend[0]._x0_in - getattr(at_sbend[0], 'X0ref', 0.0) would not work ?

if k1 == 0, normally there is no reason to call rbend_tune(). I told you that however, translating a uniform magnetic field does not change anything: that's true only if the faces are parallel, otherwise the magnet length changes. So I will change rbend_tune() to set X0Ref to 0 if K1 == 0. With this your test case 2 will be identical to test case 1.

Remain the problems for K1 != 0.

[lfarv]

Now rbend_tune will always keep X0Ref at 0 if K1 ==0. This ensures a correct reference trajectory whatever the face angles are.

However I found a problem: whatever K1 and face angles, after calling rbendtune(), the output orbit in AT is in the 1.e-15 range. Applying the conversion formula for rbend_shift, I find that:

• if face angles are 0.5 x bending angle, the output orbit in Xsuite is in the 1.e-12 range: significantly larger, but good enough to get a reasonably good agreement on the plots.
• if the face angles are 0, the output orbit in Xsuite is in the 0.02 range, far from zero: the conversion formula does not work. So disagreement between AT and Xsuite comes from the tuning of rbend_shift. Apart from building a loop to cancel the output orbit, as it is done in AT, I do not know what to do…

[swhite2401]

• if the face angles are 0, the output orbit in Xsuite is in the 0.02 range, far from zero: the conversion formula does not work. So disagreement between AT and Xsuite comes from the tuning of rbend_shift. Apart from building a loop to cancel the output orbit, as it is done in AT, I do not know what to do…

We are progressing!

I did write a loop to cancel the output orbit in XSuite but could not reach agreement either.... however, you have made a lot of changes since then. Let me try to run these tests again.

[swhite2401]

@lfarv, I confirm that even with a numerical loop on the xsuite bend agreement is still not good for K1!=0 and face angles = 0.0

[swhite2401]

@lfarv, I may have some news concerning the issue with dx . Looking at the reference:

https://www.sciencedirect.com/science/article/pii/S0168900222007793?ref=cra_js_challenge&fr=RR-1

we are implementing the linearized map and what we observe is rather good agreement up to a certain amplitude which could indicate that the exact map would be needed for large amplitude.
The paper mentions issue with the exact map relating to sqrt of negative numbers in case of "extreme" misalignment, not sure what extreme means in this case.

I think it would be worth implementing the exact map, this could be done in separate PR, @SebastienJoly when you implemented this transformations, did you by any chance take a look at this exact map?

[swhite2401]

For reference this is the kind of error we get for dx=1mm

[lfarv]

@swhite2401 : I came to the same conclusion. Implementing the exact map for dx, dy, dz is trivial: it is what is done to apply X0Ref ! It does not use the R1/2 and T1/2 matrices. But I do not know what to do with rotations.

However from your plots, it looks like for position within 5 cm (which is huge !), the error is within 1.e-7, so it's not so critical.

[swhite2401]

@lfarv, some more observation with B0 this time. I have tried to simplify the problem as much as possible.
For this I use exact sector bend with face angles=0.0, K1=0 and Fringe*=0.0 so basically only the body is considered, in XSuite this results in edge mode=suppressed.

In this case I still don't get agreement between the 2 codes, which seems very strange to me because I thought the implementation of the body was solid and well benchmarked....to be followed, I will do the same test with the small angle approximation passmethod just to be sure.

[SebastienJoly]

@swhite2401 I did not try implementing the exact map. I chose to implement the linearized map to be able to compare the results with SC/pySC, which also rely on it.
Up to +- 5 cm, what is the order of magnitude of the error compared to Xsuite?

[swhite2401]

@swhite2401 I did not try implementing the exact map. I chose to implement the linearized map to be able to compare the results with SC/pySC, which also rely on it. Up to +- 5 cm, what is the order of magnitude of the error compared to Xsuite?

It is very small, we are talking 1.0e-7 difference on the norm of the phase space coordinates, but since we are going for this one to one agreement, we are going through every possible differences. As @lfarv said this is not a problem for now, it is already quite good. But when on of us as time it should be rather simple to implement the exact map.