def interpolate_cubic(self, xi, yi, zdata, gradz=None, derivatives=False):
        """
        Cubic spline interpolation/extrapolation to arbirary point(s).
        This method has C^1 continuity.

        Parameters
        ----------
         xi : float / array of floats, shape (l,)
            x coordinates on the Cartesian plane
         yi : float / array of floats, shape (l,)
            y coordinates on the Cartesian plane
         zdata : array of floats, shape (n,)
            value at each point in the triangulation
            must be the same size of the mesh
         gradz (optional) : array of floats, shape (2,n)
            derivative at each point in the triangulation in the
            x-direction (first row), y-direction (second row)
            if not supplied it is evaluated using self.gradient
         derivatives (optional) : bool (default: False)
            optionally returns the first derivatives at point(s) (xi,yi)

        Returns
        -------
         zi : float / array of floats, shape (l,)
            interpolated value(s) of (xi,yi)
         err : int / array of ints, shape (l,)
            whether interpolation (0), extrapolation (1) or error (other)
         dzx, dzy (optional) : float, array of floats, shape(l,)
            first partial derivatives in x and y direction at (xi,yi)
        """

        if zdata.size != self.npoints:
            raise ValueError('zdata should be same size as mesh')

        if type(gradz) == type(None):
            gradX, gradY = self.gradient(zdata)
            gradX, gradY = self._shuffle_field(gradX, gradY)
        elif np.array(gradz).shape == (2,self.npoints):
            gradX, gradY = self._shuffle_field(gradz[0], gradz[1])
        else:
            raise ValueError("gradz must be of shape {}".format((2,self.npoints)))

        iflgs = 0
        dflag = 1
        sigma = 0.0


        xi = np.array(xi)
        yi = np.array(yi)

        size = xi.size

        zi = np.empty(size)
        dzx = np.empty(size)
        dzy = np.empty(size)
        zierr = np.empty(size, dtype=np.int)

        gradZ = np.vstack([gradX, gradY])
        zdata = self._shuffle_field(zdata)

        for i in range(0, size):
            ist = np.abs(self._x - xi[i]).argmin() + 1
            zi[i], dzx[i], dzy[i], zierr[i] = _srfpack.intrc1(xi[i], yi[i], self._x, self._y, zdata,\
                               self.lst, self.lptr, self.lend, iflgs, sigma, gradZ, dflag, ist)

        if derivatives:
            return zi, zierr, (dzx, dzy)
        else:
            return zi, zierr