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...

[lfarv]

@swhite2401 : I reproduced your test setup and here are my findings, taking into account these 2 problems in Xsuite:

1. The quadrupole wedge is not taken into account in Xsuite when edge_*_model == "dipole_only", while it is in AT (unless we decide to change that)
2. The multipole fringe field is wrong in Xsuite for electrons, while it should be independent of particle.

For sector magnets (ExactSectorBendPass) with K1 == 0.0

• the magnet core agrees within 1.0e-13,
• the face angle agrees within 1.0e-13,
• the dipole fringe ("dipole-only") agrees within 1.0e-13, including with fringe field integrals,
• the multipole rings field agrees within 1.0e-13,
• the dipole fringe field ("linear") also agrees within 1.0e-13. This is not normally available in AT, I made a specific branch to test that.

For sector magnets (ExactSectorBendPass) with K1 != 0.0

• For the magnet core, the agreement goes down to 1.0e-7 instead of 1.0e-13, but with a much better large central area.
• with face angles, the agreement falls to 1.0e-5, with the same square central area

For rectangular magnets (ExactRectangularBendPass)

Agreement only with K1 == 0.0 and face_angle == 0.5*bending_angle. There is clearly a problem there.

For sector magnets (BndMPoleSymplectic4Pass)

No result yet, I have to compare with the Xsuite "drift-kick-drift-expanded" model.

[swhite2401]

@lfarv, thanks for the summary, sorry I had forgotten about this charge issue... I agree with your findings!
To complement, the problem with the rectangular bend faces is most likely from XSuite, to probe this I have compared XSuite sector and rectangular bends with different face angles and only the case with 0.5*angle matches...

[swhite2401]

@lfarv see issue #1076 from @gubaidulinvadim maybe it could be integrated in there?

[swhite2401]

@lfarv we had a discussion with CERN colleagues yesterday and made some prgress see below.

@swhite2401 : I reproduced your test setup and here are my findings, taking into account these 2 problems in Xsuite:

1. The quadrupole wedge is not taken into account in Xsuite when edge_*_model == "dipole_only", while it is in AT (unless we decide to change that)
   The quadrupole wedge effect was recently implemented and this is a way for them to turn it off, they admit this is not very clean. A dedicated switch should be better, shall we implement it as well?

2. The multipole fringe field is wrong in Xsuite for electrons, while it should be independent of particle.

They agree on that and will fix this

For sector magnets (ExactSectorBendPass) with K1 == 0.0

• the magnet core agrees within 1.0e-13,
• the face angle agrees within 1.0e-13,
• the dipole fringe ("dipole-only") agrees within 1.0e-13, including with fringe field integrals,
• the multipole rings field agrees within 1.0e-13,
• the dipole fringe field ("linear") also agrees within 1.0e-13. This is not normally available in AT, I made a specific branch to test that.

For sector magnets (ExactSectorBendPass) with K1 != 0.0

• For the magnet core, the agreement goes down to 1.0e-7 instead of 1.0e-13, but with a much better large central area.
• with face angles, the agreement falls to 1.0e-5, with the same square central area

We could not reproduce this, for us the sector bend agrees well for all face angles

For rectangular magnets (ExactRectangularBendPass)

Agreement only with K1 == 0.0 and face_angle == 0.5*bending_angle. There is clearly a problem there.

Here there is some subtelties related to X0ref, to match XSuite and AT.

Without running rbendtune (X0ref=0): xs_sbend[0].rbend_shift += xs_sbend[0]._x0_in
With this we have agreement for all face angle and K1=0

With running rbendtune (X0ref!=0): xs_sbend[0].rbend_shift += xs_sbend[0]._x0_in - at_sbend[0].X0ref
With this we have agreement for parallel faces and K1!=0

For sector magnets (BndMPoleSymplectic4Pass)

No result yet, I have to compare with the Xsuite "drift-kick-drift-expanded" model.

[swhite2401]

Remaining issue:

B0, A0: They include these components in the multipole fringe field calculation, while we have the skip_b0 option activated. When using dipole fringe only I get agreement with A0 but not B0 (as expected)

dx: this is strange we get agreement only with paralell faces, there is something fishy there either for us or for them...do we neglect the change in dipole length when applying a horizontal offset in a sector bend?

X0ref: we apply the shift only on the body of the magnet while it seem that rbendshift applies an offset on everything including edges (I checked this by using parallel faces and turning off multipole fringe, in this case the 2 codes agree)

I think this is all... there are so many cases I am getting confused myself with all these tests.... but we are slowly converging!

[lfarv]

@swhite2401 : very nice, we are progressing. In the meantime, I made some additional tests:

For sector magnets (BndMPoleSymplectic4Pass)

Perfect agreement with the drift-kick-drift-expanded model, in all conditions of K1, face angles and "linear" fringe field.

For quadrupoles (StrMPoleSymplectic4Pass)

Perfect agreement with the drift-kick-drift-expanded model. A very small deviation in the fringe field model can be easily corrected. dx and dy are also perfectly reproduced

For quadrupoles (ExactMultipolePass)

Perfect agreement with the bend-kick-bend model.

So I found no additional problems compared to your list above.

[swhite2401]

Perfect agreement with the bend-kick-bend model.

Great! Just for my understanding, I would have expected this to match with kick-drift-kick-exact, why do you need to bend in a multipole?

[lfarv]

Great! Just for my understanding, I would have expected this to match with kick-drift-kick-exact, why do you need to bend in a multipole?

Sorry, it's a typo, you are right!

[lfarv]

X0ref: we apply the shift only on the body of the magnet

You are right. it should be applied to the whole magnet. I just pushed a correction for that. Let me know what you find.

[swhite2401]

X0ref: we apply the shift only on the body of the magnet

You are right. it should be applied to the whole magnet. I just pushed a correction for that. Let me know what you find.

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

[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