Doubtfull implementation of solenoid

[kyrsjo]

The variable ejf0v is used, but not updated, in the solenoid kick blocks kickvso1 and kickvso2

[rdemaria]

It indeed an error in the formula, ejf0v has been mistaken for P0c since Pc is ejfv. One should replace the ratio in the formula using oidpsv

[vikasnt]

After examining the code in kickvso1 and comparing with thin lens map:

1. ejf0v(j) should be replaced by e0f(j) as commented already in code.
2. crkve=yv(1,j)-(((xv(1,j)*strackx(i))*strackz(i))*ejf0v(j))/ejfv(j) !hr02 cikve=yv(2,j)-(((xv(2,j)*strackx(i))*strackz(i))*ejf0v(j))/ejfv(j) !hr02 and later yv(1,j)=crkve*cos_rn((strackz(i)*ejf0v(j))/ejfv(j))+ &!hr02 &cikve*sin_rn((strackz(i)*ejf0v(j))/ejfv(j)) !hr02 yv(2,j)=cikve*cos_rn((strackz(i)*ejf0v(j))/ejfv(j))- &!hr02 &crkve*sin_rn((strackz(i)*ejf0v(j))/ejfv(j)) Here it looks like yv(1,j) correspond to P_x when compared to map but Twiki writes yv(1,j) equal to P_x/P.
3. yv(1,j)=yv(1,j)-xv(2,j)*strackx(i) yv(2,j)=yv(2,j)+xv(1,j)*strackx(i) this code at start of block doesn't seem to exist in map which at the end of block is negated using following code yv(1,j)=yv(1,j)+xv(2,j)*strackx(i) yv(2,j)=yv(2,j)-xv(1,j)*strackx(i) I am unable to understand the reason of doing this.
4. cos_rn((strackz(i)ejf0v(j))/ejfv(j) . So, theta =(strackz(i)*ejf0v(j))/ejfv(j)) from code, but in map it is R/(1+delta) and 1/(1+delta)= e0f(j)/efjv(j) and assuming point no1. is correct, which makes R=strackz(i). But, from crkve=yv(1,j)-(((xv(1,j)*strackx(i))*strackz(i))*ejf0v(j))/ejfv(j) !hr02 cikve=yv(2,j)-(((xv(2,j)*strackx(i))*strackz(i))*ejf0v(j))/ejfv(j) !hr02, and when compared to map and assuming problem in no 2 is in Twiki, Px- theta * R * x= crkve so R=strackx(i) does that mean R can have different values? for theta calculation strackz(j), otherwise strackx(j).

EDIT by Kyrre: Formatting, adding code markers.

[kyrsjo]

Hi, thanks for posting this here.

First, regarding the first comment: e0f is a scalar, not a particle array - see Riccardo's comment above yours.

I think it is also important to be clear where the input variables strackx and strackz come from -- I find this block of code in trauthin and trauthck:

+cd solenoid
        strack(i)=zero
        strackx(i)=ed(IX)
        strackz(i)=ek(IX)

Here ed and ek are column 3 and 4 in the SINGLE ELEMENTS part of DATEN. These would be the inputs...

I am not sure what the cikve temp varialbe is supposed to mean, but I agree that it is very strange that it is defined first from yv(1,j), and then overwritten from yv(2,j)...

The yv arrays in SixTrack is Px/P and Py/P, which is equivalent to the angle the particle is traveling on relative to the axis of the element, i.e. dx/ds and dy/ds, in milli-radians.

Riccardo, you're very welcome to join the discussion!

[vikasnt]

About the e0f, by mistake i wrote array index while writing comment but i agree it should be scalar. As for yv(1,j)=Px/P, then the equations weren't matching, but if you are sure it should be the ratio, then i'll check again. And i don't think cikve or crkve are being overwritten without being used first, and I think their behavior is ok ( we can discuss the code if you are still unsure about it). Also, ed(IX).. still doesnt explain why i get different values for R ( I would assume it to be constant?)

[kyrsjo]

Right, cikveand crkve are obviously different variables. I just missed the change of i to r.

[vikasnt]

Update(On my first comment): I checked point no. 2 again regarding yv(1,j)=Px/P. From Physics manual, and while from Twiki yv(1,j)=px/p. and in Sixtrack, crkve=yv(1,j)-(((xv(1,j)*strackx(i))*strackz(i))*ejf0v(j))/ejfv(j) cikve=yv(2,j)-(((xv(2,j)*strackx(i))*strackz(i))*ejf0v(j))/ejfv(j) yv(1,j)=crkve*cos_rn((strackz(i)*ejf0v(j))/ejfv(j))+ cikve*sin_rn((strackz(i)*ejf0v(j))/ejfv(j)) Assuming ejf0v(j) is replaced by e0f, 1.Twiki should use yv(1,j)=Px/P inplace of yv(1,j)=px/p to avoid confusion since in manual px=Px/P0.

2. Although the above Sixtrack code is correct, on both side a factor of P0/P is missing as Manual uses px=Px/P0 while code uses yv(1,j)=Px/P(which anyway would get canceled as a constant). Last 2 issues ( in my first comment) are still unresolved.

[kyrsjo]

Working through your comments and the math now - one quick comment is that the comment saying " TODO: Check if ejf0v should be e0f?? or oidpsv=ejf0v(j)/ejfv(j)=1/(1+delta)" is not entirely correct - the last equation should be "oidpsv(j)=e0f/ejfv(j)=1/(1+delta)".

[kyrsjo]

Regarding your comment no. 3, note that the modified 'yv' variables are used to calculated 'crkve' and 'cikve', which is then used to update the yv variables, before "undoing" the transformation of the first two lines (but now using updated x and y).

Regarding problem no. 4, I suspect there may be a problem in the physics manual. The length of the solenoid (for the integrated strength) should be in the equations somewhere, but it is not. Maybe if the "theta * R" in the physics manual should actually be theta * L?, with L=strackx? Alternatively, if theta=L_R/(1+\delta), with L_R = kz, and R=strackx, we would be pretty similar to the solenoid code in MadX, with R = K_s/(2*c).

The relevant kick code in MadX, from file trrun.f90:

subroutine trsol(track,ktrack)
  use math_constfi, only : one, two, half
  implicit none
  !----------------------------------------------------------------------*
  ! Purpose:                                                             *
  !   Track a set of particles through a thin solenoid.                  *
  ! Input/output:                                                        *
  !   TRACK(6,*)(double)    Track coordinates: (X, PX, Y, PY, T, PT).    *
  !   KTRACK    (integer)   number of surviving tracks.                  *
  ! Output:                                                              *
  !   EL        (double)    Length of solenoid.                          *
  !----------------------------------------------------------------------*
  double precision :: track(6,*)
  integer :: ktrack

  integer :: i
  double precision :: beta
  double precision :: sk, skl, sks, sksl, cosTh, sinTh, Q, R, Z
  double precision :: xf, yf, pxf, pyf, sigf, psigf, bvk
  double precision :: onedp, fpsig, fppsig

  double precision :: get_value, node_value

  !---- Initialize.
  ! dtbyds  = get_value('probe ','dtbyds ')
  ! gamma   = get_value('probe ','gamma ')
  beta = get_value('probe ','beta ')
  ! deltap  = get_value('probe ','deltap ')
  !
  !---- Get solenoid parameters
  ! elrad   = node_value('lrad ')
  sksl = node_value('ksi ')
  sks  = node_value('ks ')

  !---- BV flag
  bvk = node_value('other_bv ')
  sks  = sks  * bvk
  sksl = sksl * bvk

  !---- Set up strengths
  ! sk    = sks / two / (one + deltap)
  sk    = sks / two
  skl   = sksl / two

  !---- Loop over particles
  do  i = 1, ktrack
     !     Ripken formulae p.28 (3.35 and 3.36)
     xf    = track(1,i)
     yf    = track(3,i)
     psigf = track(6,i) / beta

     !     We do not use a constant deltap!!!!! WE use full 6D formulae!
     onedp   = sqrt( one + 2*psigf + (beta**2)*(psigf**2) )
     fpsig   = onedp - one
     fppsig  = ( one + (beta**2)*psigf ) / onedp

     !     Set up C,S, Q,R,Z
     cosTh = cos(skl/onedp)
     sinTh = sin(skl/onedp)
     Q = -skl * sk / onedp
     R = fppsig / (onedp**2) * skl * sk
     Z = fppsig / (onedp**2) * skl

     pxf  = track(2,i) + xf*Q
     pyf  = track(4,i) + yf*Q
     sigf = track(5,i)*beta - half*(xf**2 + yf**2)*R

     !       Ripken formulae p.29 (3.37)
     track(1,i) =  xf  * cosTh  +  yf  * sinTh
     track(2,i) =  pxf * cosTh  +  pyf * sinTh
     track(3,i) = -xf  * sinTh  +  yf  * cosTh
     track(4,i) = -pxf * sinTh  +  pyf * cosTh
     track(5,i) =  (sigf + (xf*pyf - yf*pxf)*Z) / beta
     ! track(6,i) =  psigf*beta

  enddo
end subroutine trsol

I'll dig out the madx->sixtrack conversion code in a minute

[kyrsjo]

More MadX code, this time from mad_6track.c:

  else if ((strcmp(p->base_name,"solenoid") == 0))
  {
    c6t_elem = new_c6t_element(11,t_name,p->base_name);
    clean_c6t_element(c6t_elem);
    strcpy(c6t_elem->org_name,t_name);
    c6t_elem->value[0] = el_par_value_recurse("l",p->p_elem);
    c6t_elem->value[2] = el_par_value_recurse("ks",p->p_elem);
    c6t_elem->value[3] = el_par_value_recurse("ksi",p->p_elem);
  }

and

static void
att_solenoid(struct c6t_element* el)
{
  el->out_1 = 25;
  el->out_2 = el->value[2];
  el->out_3 = el->value[3];
  el->out_4 = el->value[0];
}

I'm posting this code here, because the most typical way of making a fort.2 is to take a machine description from MadX (our main beam dynamics program), and call a madx to sixtrack conversion routine from there which converts it and writes a fort.2. This routine is not perfect, but it's often a good idea to get an idea of what the data flow looks like.

This is the code that actually writes the fort.2 (except for rf multipoles):

static void
write_c6t_element(struct c6t_element* el)
{
  if (strcmp(el->name, "CAV") != 0) {
    if (strcmp(el->base_name, "rfmultipole")==0) { 
      write_rfmultipole(el);
    } else {
      fprintf(f2, "%-16s %3d  %16.9e %17.9e  %17.9e  %17.9e  %17.9e  %17.9e\n",
          el->name, el->out_1, el->out_2, el->out_3, el->out_4, el->out_5, el->out_6, el->out_7);
    }
  }
  el->w_flag = 1;
}

So we see that it writes the name, element type (kz), ed, ek, el, bbbx, bbby, bbbz by comparison to the line in sixtrack.s/subroutine daten which actually reads the line:

      read(ch1,*) idat,kz(i),ed(i),ek(i),el(i),bbbx(i),bbby(i),bbbs(i)!read fort.2 (or fort.3), idat -> bez = single element name, kz = type of element, ed,ek,el = strength, random error on strenght,length (can be anything),bbbx,bbby,bbbs = beam-beam, beam-beam parameters will be removed soon.

Thus, strackx(i)=ed(IX)= el->out_2 = el->value[2] = el_par_value_recurse("ks",p->p_elem); i.e. strackx is "ks" in MadX, which in the manual is defined as the bending strength B0/B_rho, where B_rho is the rigidity of the beam B*rho = p_s/q.

Further, strackz(i)=ek(IX)= el->out_3 = el->value[3] -> el_par_value_recurse("ksi",p->p_elem), which is Ks*L in MadX.

MadX user manual: http://cern.ch/madx/releases/last-dev/madxuguide.pdf

[kyrsjo]

My notes when writing out the equations in kickvso1:

The kickvso1 block can be in mathematical notation be written as follows. I have here assumed that ejf0v(j)/ejfv(j) should actually be oidpsv(j) = 1/(1+delta). Further, I have simplified the variable names strackx -> kx, strackz -> kz, crkve -> cr, and cikve ->ci. The first part of the transformation is: Second part: Third part, first within the crlibm/not crlibm block: Forth part, still within crlibm/not crlibm: Then outside of the crlibm/not crlibm: And then, finally:

In kickvso2, this is followed an update of sigmv (but not the other longitudinal variables...).

[kyrsjo]

OK, just to make it clear: The "Ripken formulae" referred to in the MadX code are from this document: http://arxiv.org/pdf/acc-phys/9510005.pdf

Some more info about the MadX implementation: https://ab-dep-abp.web.cern.ch/ab-dep-abp/LCU/LCU_meetings/2005/251005/AK.pdf

If we can verify that the MadX implementation is the same as written in "A symplectic six-dimensional thin-lens formalism for tracking", and maybe also tease out any typos in the physics manual, that would be great. The next step would then be to compare it with the other document http://arxiv.org/pdf/acc-phys/9510005.pdf

If they all agree, and the SixTrack code is different, then we shall contact Frank etc. to discuss what has happened.

[kyrsjo]

More info on the MadX solenoid implementation: http://accelconf.web.cern.ch/AccelConf/e06/PAPERS/WEPCH043.PDF

[rdemaria]

Just to add that in eq. 9

delta=f(p_sigma) f'(p_sigma)=beta_0/beta=rvv in sixtrack

On Thu, Aug 11, 2016 at 5:24 PM, Kyrre Ness Sjøbæk <notifications@github.com

wrote:

More info on the MadX solenoid implementation: http://accelconf.web.cern.ch/AccelConf/e06/PAPERS/WEPCH043.PDF

— You are receiving this because you commented. Reply to this email directly, view it on GitHub https://github.com/SixTrack/SixTrack/issues/4#issuecomment-239195306, or mute the thread https://github.com/notifications/unsubscribe-auth/AADUYdjALwrbMK91WtYTBggz0tnj2PDiks5qez6cgaJpZM4Fyev_ .

[vikasnt]

After comparing equations in http://accelconf.web.cern.ch/AccelConf/e06/PAPERS/WEPCH043.PDF (this one has everything nicely put in a single place) and https://cds.cern.ch/record/281283/files/sl-95-012.pdf and http://arxiv.org/pdf/acc-phys/9510005.pdf with Sixtrack Physics Manual,

1. theta=R*L/(1+delta) while manual currently states theta=R/1+delta
2. equations are not correct if px=Px/Po ( which is written in Sixtrack manual but not in other files) unless other files are also assuming this but haven't written it explicitly.Even if the equations are same, Sixtrack code implementation is still wrong since px(manual) to x'(sixtrack yv(1,j) variable) transformation is not correct.( I earlier wrote 2. Although the above Sixtrack code is correct, on both side a factor of P0/P is missing as Manual uses px=Px/P0 while code uses yv(1,j)=Px/P(which anyway would get canceled as a constant). but this was a mistake, as it won't cancel out)
3. R=ks/2 and R*L=ks*L/2=ks_i/2 from both http://accelconf.web.cern.ch/AccelConf/e06/PAPERS/WEPCH043.PDF and Sixtrack Manual but from @kyrsjo earlier comment, Thus, strackx(i)=ed(IX)= el->out_2 = el->value[2] = el_par_value_recurse("ks",p->p_elem); i.e. strackx is "ks" in MadX, which in the manual is defined as the bending strength B0/Brho, where Brho is the rigidity of the beam B*rho = p_s/q. Further, strackz(i)=ek(IX)= el->out_3 = el->value[3] -> el_par_value_recurse("ksi",p->p_elem), which is Ks*L in MadX. Thus it looks like inputs are off by a factor of 1/2, since we should have strackx(i)=R=ks/2 and strackz(i)=R*L=ks*L/2.
4. yv(1,j)=yv(1,j)-xv(2,j)*strackx(i) and yv(2,j)=yv(2,j)+xv(1,j)*strackx(i) still doesn't make sense and if we simply neglect these equation then sixtrack implementation looks almost correct (ofcourse if we also replace ejf0v-->e0f and point 2 is taken care of).

[kyrsjo]

Hi! Great that you could look more at this.

Regarding 1), I agree. This should be fixed in the manual.

2) Agree. I think the easiest way to fix this would be to transform the x' -> Px, do the calculation, then transform it back again. And then there is the ej0v vs. e0f confusion as you mention in 4). Note that the ratio is already pre-calculated and stored in oidpsv(j)=e0f/ejfv(j)=1/(1+delta), so you could use that.

3) I suspected also that some factor 1/2 was missing somewhere. This should be put in (maybe in +cd solenoid?) and documented in the physics manual. Then how to use the element should be documented in the user manual.

4) (After discussing with Riccardo on my end) These may be due to that the canonical momentum Px includes the vector potential of the field \vec A. So when you enter the solenoid, the field changes, meaning that the canonical momentum also changes. However, as discussed in section 3.5 in http://www.worldscientific.com/doi/pdf/10.1142/9781783262786_0003 [1], this dependence may cancel due to the effect of the fringe field. This should be discussed in the physics manual.

It seems that thing are getting clearer :)

[1] Its a chapter from the book "Beam Dynamics in High Energy Particle Accelerators" by Andrzej Wolski; you should have access to it through CERN. If you have any problems, please contact me on slack. http://www.worldscientific.com/worldscibooks/10.1142/p899#t=toc

[JamesMolson]

@tpersson I assume this was fixed by #543 and can be closed?