def interpolate_linear(self, xi, yi, zdata):
        """
        Piecewise linear interpolation/extrapolation to arbitrary point(s).
        The method is fast, but has only C^0 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

        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)
        """

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

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

        size = xi.size

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

        zdata = self._shuffle_field(zdata)

        # iterate
        for i in range(0, size):
            ist = np.abs(self._x - xi[i]).argmin() + 1
            zi[i], zierr[i] = _srfpack.intrc0(xi[i], yi[i], self._x, self._y, zdata,\
                                       self.lst, self.lptr, self.lend, ist)

        return zi, zierr