New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Fix multipole fields #322

Merged
merged 4 commits into from Mar 7, 2018

Conversation

Projects
None yet
2 participants
@danieljeans
Contributor

danieljeans commented Feb 28, 2018

I believe there are bugs in the definition of the multipole fields (both in the code and the documentation).
I've compared to the definitions given in
https://dd4hep.web.cern.ch/dd4hep/reference/classdd4hep_1_1MultipoleField.html

BEGINRELEASENOTES

  • correct coded expressions for quadrupole, sextupole and octopole fields.
  • correct expression for octopole field in documentation
    ENDRELEASENOTES

@petricm petricm added the bug fix label Feb 28, 2018

bx += coefficents[3] * (3.0*x2*y - y2*y) + skews[3]*(x2*x - 3.0*x*y2);
[[fallthrough]];
case 3: // Sextupole momentum:
by += -coefficents[2] * (x2 - y2) + skews[2] * 2.0 * xy;
by += coefficents[2] * (x2 - y2) - skews[2] * 2.0 * xy;
bx += coefficents[2] * 2.0 * xy + skews[2] * (x2 - y2);

This comment has been minimized.

@petricm

petricm Feb 28, 2018

Member

Shouldn't this be

by +=  0.5 * coefficents[2] * (x2 - y2) - skews[2] * xy;
bx +=  coefficents[2] * xy + 0.5 * skews[2] * (x2 - y2);

to have the coefficient[1] vs. coefficient[2] properly normalised?

@@ -99,16 +99,16 @@ void MultipoleField::fieldComponents(const double* pos, double* field) {
double y2 = y*y;
switch(coefficents.size()) {
case 4: // Ocupole momentum
by += coefficents[3] * (x2*x - 3.0*x*y2) + skews[3]*(y2*y - 3.0*x*y2);

This comment has been minimized.

@petricm

petricm Feb 28, 2018

Member

I think this fix is correct, but the same comment as below, shouldn't there be a 1/6 in front for normalization?

@danieljeans

This comment has been minimized.

Contributor

danieljeans commented Mar 1, 2018

I believe that what I implemented in the code is consistent with the definition in the documentation in the header file, which is in turn consistent with the references to which it links: for the nth multipole,
By + i Bx = ( b + i a ) ( x + i y ) ^ (n-1). [b is the "coefficient", a is the "skew".]
(eg see p10 of http://cas.web.cern.ch/sites/cas.web.cern.ch/files/lectures/bruges-2009/wolski-1.pdf)
I'm sure that definitions with different normalisations are possible, but this one seems standard at least in the accelerator community.

@petricm

This comment has been minimized.

Member

petricm commented Mar 1, 2018

I cannot agree with your statement

this one seems standard at least in the accelerator community.

If you look at this lecture from another accelerator school you find an explicit definition for the terms (page 23) which give you different normalization (the form of the expression is however the same).

@danieljeans

This comment has been minimized.

Contributor

danieljeans commented Mar 1, 2018

@petricm thanks for pointing this out to me.
I have no preference as to which convention to use, as long as it's consistent and well documented.

@petricm

This comment has been minimized.

Member

petricm commented Mar 1, 2018

I'll ask someone in the magnets department for clarification and let's see what they say...

@petricm

This comment has been minimized.

Member

petricm commented Mar 1, 2018

I still did not receive a reply but for book reference you can look at this book with the definition the same as the USPAS slides.

@petricm

This comment has been minimized.

Member

petricm commented Mar 5, 2018

I have received a reply that we should consider Field Computation for Accelerator Magnets by Stephan Russenschuck as the reference for all matters accelerator magnets. If you wish you can download the book inside the CERN network here.

The conclusion is that a sextupole is parameterised as

Bx(x,y)= 1/r_0^2 ( A3(r_0) (x^2-y^2) + 2B3(r_0)xy)

where r_0 is an arbitrary reference radius, so the boundary conditions defines you normalization.

If you however do a multipole expansion as our code does, so add a dipole field to a quadropole, then there is a 1/(n-1)! term. Chapter 6 discussion from (6.80) to (6.83) and the definition of normalized gradient.

So I would like to kindly ask you if you can fix the code to include this.

@danieljeans

This comment has been minimized.

Contributor

danieljeans commented Mar 6, 2018

thanks for the info. It's easy to include the (n-1)! scaling, I'll do it right away.

if you want the ability to explicitly include the r_0 normalisation (ie to be able to specify the coefficients and skews at any particular reference radius), we would need an extra parameter.
In my opinion, we can explicitly say that r_0 = 1 and avoid this complication: do you agree?

@danieljeans danieljeans force-pushed the danieljeans:fixMultipoleFields branch from 2f66894 to e0b9630 Mar 6, 2018

@danieljeans

This comment has been minimized.

Contributor

danieljeans commented Mar 6, 2018

by the way, I could not get my hands on the book you referenced: the KEK copy is out on loan, and I cannot get into the cern network: so please check that what I implemented looks as expected.

@petricm petricm merged commit 91a695f into AIDASoft:master Mar 7, 2018

1 check passed

continuous-integration/travis-ci/pr The Travis CI build passed
Details
@petricm

This comment has been minimized.

Member

petricm commented Mar 7, 2018

I went through the equations again and they are correct.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment