fix rf multipole [wip]

examples/rf_multipole/rf_multipole.py

@@ -7,33 +7,25 @@
 # run MADX tests
 mad = Madx()
 mad.call("rf_multipole.madx")
-tracksumm_madx = mad.table.tracksumm
 
 # same test in pysixtrack
 mad_sequence = mad.sequence["sequ_rfmultipole"]
 rf_mulitpole_mad = mad_sequence.elements[1]
 freq = rf_mulitpole_mad.freq * 1e6  # MAD units are MHz
 knl = rf_mulitpole_mad.knl
-pnl = np.array(rf_mulitpole_mad.pnl) * 360  # MAD units are 2pi
+pn = np.array(rf_mulitpole_mad.pnl) * 360  # MAD units are 2pi
 lag = rf_mulitpole_mad.lag * 360  # MAD units are 2pi
 
-my_rf_multipole = pysixtrack.elements.RFMultipole(
-    voltage=0, frequency=freq, lag=lag, knl=knl, ksl=[0], pn=pnl, ps=[0]
+rf_multipole = pysixtrack.elements.RFMultipole(
+    voltage=0, frequency=freq, lag=lag, knl=knl, ksl=[0], pn=pn, ps=[0]
 )
 
-p0c = mad_sequence.beam.pc * 1e9
-x = tracksumm_madx.x[0]
-px = tracksumm_madx.px[0]
-y = tracksumm_madx.y[0]
-py = tracksumm_madx.py[0]
-t = tracksumm_madx.t[0]
-pt = tracksumm_madx.pt[0]
+mad_part=Particles.from_mad_track(mad)
+p1=mad_part.copy(0)
+p2=mad_part.copy(1)
+p3=p1.copy()
+rf_multipole.track(p3)
 
-part = Particles(p0c=p0c, x=x, px=px, y=y, py=py, tau=t, pt=pt)
-# print(part)
+p2.compare(p3)
 
-my_rf_multipole.track(part)
-print(part)
 
-# part(p0c=p0c, x=x, px=px, y=y, py=py, zeta=?, delta=?)
-""

pysixtrack/elements.py

@@ -149,6 +149,13 @@ def track(self, p):
 
 
 class RFMultipole(Element):
+    """
+    H= -l sum   Re[ (kn[n](zeta) + i ks[n](zeta) ) (x+iy)**(n+1)/ n ]
+
+    kn[n](z) = k_n cos(2pi w tau + pn/180*pi)
+    ks[n](z) = k_n cos(2pi w tau + pn/180*pi)
+
+    """
     _description = [
         ("voltage", "volt", "Voltage", 0),
         ("frequency", "hertz", "Frequency", 0),
@@ -178,10 +185,6 @@ def track(self, p):
         ps = _arrayofsize(self.ps, order + 1) * deg2rad - ktau
         x = p.x
         y = p.y
-        # fnr = knl[0]
-        # fni = 0
-        # fsr = ksl[0]
-        # fsi = 0
         dpx = 0
         dpy = 0
         dptr = 0
@@ -209,8 +212,8 @@ def track(self, p):
         chi = p.chi
         p.px += -chi * dpx
         p.py += chi * dpy
-        phase = self.lag * deg2rad - ktau
-        p.add_to_energy(chi * (self.voltage * sin(phase) + dptr))
+        dv0 = self.voltage * sin(self.lag * deg2rad - ktau)
+        p.add_to_energy(chi * (dv0 - p.p0c*k*dptr))
 
 
 class Cavity(Element):

pysixtrack/loader_mad.py

@@ -86,6 +86,17 @@ def iter_from_madx_sequence(
                 lag=ee.lag * 360,
             )
 
+        elif mad_etype == "rfmultipole":
+            newele = classes.RFMultipole(
+                voltage=ee.volt * 1e6,
+                frequency=ee.freq * 1e6,
+                lag=ee.lag * 360,
+                knl=ee.knl,
+                ksl=ee.ksl,
+                pn=ee.pnl,
+                ps=ee.psl,
+            )
+
         elif mad_etype == "beambeam":
             if ee.slot_id == 6 or ee.slot_id == 60:
                 # BB interaction is 6D

pysixtrack/particles.py

@@ -86,11 +86,14 @@ def _f4(self, mass0, gamma0):
         beta0 = sqrt(1 - 1 / gamma0 ** 2)
         return self._g2(mass0, beta0, gamma0)
 
-    def copy(self):
+    def copy(self,index=None):
         p = Particles()
         for k, v in list(self.__dict__.items()):
             if type(v) in [np.ndarray, dict]:
-                v = v.copy()
+                if index is None:
+                    v = v.copy()
+                else:
+                    v= v[index]
             p.__dict__[k] = v
         return p
 
@@ -474,10 +477,10 @@ def compare(self, particle, rel_tol=1e-6, abs_tol=1e-15):
             if v1 is not None and v2 is not None:
                 diff = v1 - v2
                 if abs(diff) > abs_tol:
-                    print(kk, v1, v2, diff)
+                    print(f"{kk} {v1} {v2}  diff:{diff}")
                     res = False
                 if abs(v1) > 0 and abs(diff) / v1 > rel_tol:
-                    print(kk, v1, v2, abs(diff) / v1)
+                    print(f"{kk} {v1} {v2} rdiff:{diff/v1}")
                     res = False
         return res
 
@@ -498,18 +501,20 @@ def from_mad_twiss(cls, twiss):
         return out
 
     @classmethod
-    def from_mad_track(cls, tracksumm):
+    def from_mad_track(cls, mad):
+        tracksumm=mad.table.tracksumm
+        mad_beam=mad.sequence().beam
         out = cls(
-            p0c=twiss.summary.pc * 1e6,
-            mass0=twiss.summary.mass * 1e6,
-            q0=twiss.summary.charge,
-            s=twiss.s[:],
-            x=twiss.x[:],
-            px=twiss.px[:],
-            y=twiss.py[:],
-            py=twiss.py[:],
-            tau=twiss.t[:],
-            ptau=twiss.pt[:],
+            p0c=mad_beam.pc * 1e6,
+            mass0=mad_beam.mass * 1e6,
+            q0=mad_beam.charge,
+            s=tracksumm.s[:],
+            x=tracksumm.x[:],
+            px=tracksumm.px[:],
+            y=tracksumm.py[:],
+            py=tracksumm.py[:],
+            tau=tracksumm.t[:],
+            ptau=tracksumm.pt[:],
         )
         return out