def getBezierPaths(self,origin=None):
        """
        This function returns array that can be used as a Cubic Bezier
        Path in matplotlib. 
        The function returns two arrays, the first one contains
        the verticies for each particles and has the shape
        (Nvert, Nparticles, 2) where Nvert is the number of verticies.
        The second array returned describes the type of verticies to be
        used with matplotlib's Patch class.

        Arguments
        ---------
        origin : multiple, optional
                 If `origin` is None (default), then none of the 
                 coordinates are shifted. If `origin` is an integer
                 then the particle with that index is used as the
                 origin. if `origin` is equal to `com`, then the 
                 centre of mass is used as the origin. 


        Examples
        --------
        The following example reads in a SimulationArchive and plots
        the trajectories as Cubic Bezier Curves. It also plots the 
        actual datapoints stored in the SimulationArchive. 
        Note that the SimulationArchive needs to have enough
        datapoints to allow for smooth and reasonable orbits.

        >>> from matplotlib.path import Path
        >>> import matplotlib.patches as patches
        >>> sa = rebound.SimulationArchive("test.bin")
        >>> verts, codes = sa.getBezierPaths(origin=0)
        >>> fig, ax = plt.subplots()
        >>> for j in range(sa[0].N):
        >>>     path = Path(verts[:,j,:], codes)
        >>>     patch = patches.PathPatch(path, facecolor='none')
        >>>     ax.add_patch(patch)
        >>>     ax.scatter(verts[::3,j,0],verts[::3,j,1])
        >>> ax.set_aspect('equal')
        >>> ax.autoscale_view()
        
        """
        import numpy as np
        Npoints = len(self)*3-2
        if len(self)<=1:
            raise Runtim
        Nparticles = self[0].N
        verts = np.zeros((Npoints,Nparticles,2))
        xy = np.zeros((len(self),Nparticles,2))

        if origin=="com":
            origin = -2
        elif origin is not None:
            try:
                origin = int(origin)
            except:
                raise AttributeError("Cannot parse origin")
            if origin<0 or origin>=Nparticles:
                raise AttributeError("Origin index out of range")


        for i, sim in enumerate(self):
            if origin is None:
                shift = (0,0,0,0)
            elif origin == -2:
                sp = sim.calculate_com()
                shift = (sp.x,sp.y,sp.vx,sp.vy)
            else:
                sp = sim.particles[origin]
                shift = (sp.x,sp.y,sp.vx,sp.vy)
            for j in range(sim.N):
                p = sim.particles[j]
                if i==0:
                    verts[0,j] = p.x-shift[0],p.y-shift[1]
                    verts[1,j] = p.vx-shift[2], p.vy-shift[3]
                else:
                    dt = sim.t-tlast # time since last snapshot
                    verts[-2+i*3,j] = verts[-2+i*3,j]*dt/3.+verts[-3+i*3,j]

                    verts[ 0+i*3,j] = p.x-shift[0],p.y-shift[1]

                    verts[-1+i*3,j] = -p.vx+shift[2], -p.vy+shift[3]
                    verts[-1+i*3,j] = verts[-1+i*3+0,j]*dt/3.+verts[ 0+i*3,j]

                    if i!=len(self)-1:
                        verts[+1+i*3,j] = p.vx-shift[2], p.vy-shift[3]


                xy[i,j] = p.x,p.y
            tlast = sim.t
        codes = np.full(Npoints,4,dtype=np.uint8) # Hardcoded 4 = matplotlib.path.Path.CURVE4
        codes[0] = 1 # Hardcoded 1 = matplotlib.path.Path.MOVETO
        return verts, codes