diff --git a/.github/workflows/pythonpublish.yml b/.github/workflows/pythonpublish.yml
new file mode 100644
index 0000000..21f2f01
--- /dev/null
+++ b/.github/workflows/pythonpublish.yml
@@ -0,0 +1,26 @@
+name: Upload Python Package
+
+on:
+  release:
+    types: [created]
+
+jobs:
+  deploy:
+    runs-on: ubuntu-latest
+    steps:
+    - uses: actions/checkout@v1
+    - name: Set up Python
+      uses: actions/setup-python@v1
+      with:
+        python-version: '3.x'
+    - name: Install dependencies
+      run: |
+        python -m pip install --upgrade pip
+        pip install setuptools wheel twine
+    - name: Build and publish
+      env:
+        TWINE_USERNAME: ${{ secrets.PYPI_USERNAME }}
+        TWINE_PASSWORD: ${{ secrets.PYPI_PASSWORD }}
+      run: |
+        python setup.py sdist bdist_wheel
+        twine upload dist/*
diff --git a/examples/petra4/track_p1.py b/examples/petra4/track_p1.py
index 4de02ce..e287453 100644
--- a/examples/petra4/track_p1.py
+++ b/examples/petra4/track_p1.py
@@ -34,7 +34,7 @@ def __repr__(self):
 print(mad.table.summ.q1, mad.table.summ.q2)
 
 
-twissout = pysixtrack.Particles.from_twiss(mad.twiss(betx=1, bety=1, x=0.001))
+twissout = pysixtrack.Particles.from_mad_twiss(mad.twiss(betx=1, bety=1, x=0.001))
 
 line = pysixtrack.Line.from_madx_sequence(mad.sequence.ring)
 part = pysixtrack.Particles()
diff --git a/examples/rf_multipole/rf_multipole.madx b/examples/rf_multipole/rf_multipole.madx
index abb342d..bb1172b 100644
--- a/examples/rf_multipole/rf_multipole.madx
+++ b/examples/rf_multipole/rf_multipole.madx
@@ -25,8 +25,6 @@ endsequence;
 use, sequence = sequ_rfmultipole;
 
 
-
-
 TRACK, deltap= 0.0, onepass, dump;
 START, x=0.1, px=0.00, y=0.05, py=0.0, t=0.5, pt=0;
 RUN, track_harmon=1e-16;
diff --git a/examples/rf_multipole/rf_multipole.py b/examples/rf_multipole/rf_multipole.py
index 43be24d..4f75dc4 100644
--- a/examples/rf_multipole/rf_multipole.py
+++ b/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=?)
-""
diff --git a/pysixtrack/elements.py b/pysixtrack/elements.py
index e920768..c6db55d 100644
--- a/pysixtrack/elements.py
+++ b/pysixtrack/elements.py
@@ -66,7 +66,6 @@ def track(self, p):
 
 
 def _arrayofsize(ar, size):
-    print(ar, size)
     ar = np.array(ar)
     if len(ar) == 0:
         return np.zeros(size, dtype=ar.dtype)
@@ -112,7 +111,6 @@ def order(self):
 
     def track(self, p):
         order = self.order
-        print(self.knl)
         length = self.length
         knl = _arrayofsize(self.knl, order + 1)
         ksl = _arrayofsize(self.ksl, order + 1)
@@ -151,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),
@@ -180,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
@@ -211,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):
diff --git a/pysixtrack/loader_mad.py b/pysixtrack/loader_mad.py
index 0fd9eb5..0639be5 100644
--- a/pysixtrack/loader_mad.py
+++ b/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
diff --git a/pysixtrack/particles.py b/pysixtrack/particles.py
index 33810e3..0c3cb1d 100644
--- a/pysixtrack/particles.py
+++ b/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,15 +477,15 @@ 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
 
     @classmethod
-    def from_twiss(cls, twiss):
+    def from_mad_twiss(cls, twiss):
         out = cls(
             p0c=twiss.summary.pc * 1e6,
             mass0=twiss.summary.mass * 1e6,
@@ -497,6 +500,24 @@ def from_twiss(cls, twiss):
         )
         return out
 
+    @classmethod
+    def from_mad_track(cls, mad):
+        tracksumm=mad.table.tracksumm
+        mad_beam=mad.sequence().beam
+        out = cls(
+            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
+
     @classmethod
     def from_list(cls, lst):
         ll = len(lst)
diff --git a/requirements.txt b/requirements.txt
index 24ce15a..a8e8446 100644
--- a/requirements.txt
+++ b/requirements.txt
@@ -1 +1,3 @@
 numpy
+scipy
+dataclasses
diff --git a/setup.py b/setup.py
index 1acbdb5..430c0c2 100644
--- a/setup.py
+++ b/setup.py
@@ -1,13 +1,23 @@
 import setuptools
+import sys
+
+requirements={
+    "install": ['numpy']
+}
+ 
+if sys.version_info < (3,7):
+    requirements['install'].append('dataclasses')
+
+version = open("pysixtrack/__init__.py").readline().split('"')[1]
 
 setuptools.setup(
     name="pysixtrack",
-    version="0.0.4",
+    version=version,
     description="6D Tracking Code",
     author="Riccardo De Maria",
     author_email="riccardo.de.maria@cern.ch",
     url="https://github.com/rdemaria/pysixtrack",
     packages=["pysixtrack", "pysixtrack.be_beamfields"],
     package_dir={"pysixtrack": "pysixtrack"},
-    install_requires=["numpy"],
+    install_requires=requirements['install'],
 )