diff --git a/particle.py b/particle.py
index a7d9484..c3af2a8 100644
--- a/particle.py
+++ b/particle.py
@@ -123,9 +123,12 @@ class Particles:
                 key = ('x','y','z')[axis]
                 self.attr[key] %= self.period[axis]
         return
-    def position_duplicates(self,ipart,padding=0.0):
-        pos  = np.array((px[ipart],py[ipart],pz[ipart]))
-        rp   = pr[ipart]+padding
+    def position_with_duplicates(self,ipart,padding=0.0):
+        pos  = np.array(
+            (self.attr['x'][ipart],
+             self.attr['y'][ipart],
+             self.attr['z'][ipart]))
+        rp = self.attr['r'][ipart]+padding
         posd = [pos.copy()]
         for axis in range(3):
             if self.period[axis] is not None:
@@ -133,41 +136,29 @@ class Particles:
                 if pos[axis]-rp<0.0:
                     for ii in range(nposd):
                         tmp = posd[ii].copy()
-                        tmp[axis] = np.mod(tmp[axis]-rp,self.period[axis])
+                        tmp[axis] = tmp[axis]+self.period[axis]
                         posd.append(tmp)
                 if pos[axis]+rp>self.period[axis]:
                     for ii in range(nposd):
                         tmp = posd[ii].copy()
-                        tmp[axis] = np.mod(tmp[axis]+rp,self.period[axis])
+                        tmp[axis] = tmp[axis]-self.period[axis]
                         posd.append(tmp)
         return posd
     def mask_field(self,fld,cval=np.nan,padding=0.0):
-      '''Fills grid points which lie inside of solid phase with values.'''
-      for ipart in range(0,part.num):
-        # Slice a box from the field around the particle
-        xp,yp,zp,rp = (px[ipart],py[ipart],pz[ipart],pr[ipart])
-        rp += padding
-        idxlo = np.array(fld.nearest_gridpoint(np.array(xp,yp,zp)-rp,lower=True))
-        idxhi = idxlo+2*rp/fld.spacing
-        
-
-        # Get bounding box of particle
-        idx_x = np.nonzero((xg>=xp-rp) & (xg<=xp+rp))[0]
-        idx_y = np.nonzero((yg>=yp-rp) & (yg<=yp+rp))[0]
-        idx_z = np.nonzero((zg>=zp-rp) & (zg<=zp+rp))[0]
-        # Triple for loop
-        for ii in range(idx_x[0],idx_x[-1]+1):
-          Dx = xg[ii]-xp
-          for jj in range(idx_y[0],idx_y[-1]+1):
-            Dy = yg[jj]-yp
-            for kk in range(idx_z[0],idx_z[-1]+1):
-              Dz = zg[kk]-zp
-              isInside = Dx*Dx+Dy*Dy+Dz*Dz <= rp*rp
-              if isInside:
-                if reconstruct:
-                  self.field[key][ii,jj,kk] = coeff_lin + coeff_rotx*Dx + coeff_roty*Dy + coeff_rotz*Dz
-                else:
-                  self.field[key][ii,jj,kk] = cval
+        '''Fills grid points which lie inside of solid phase with values.'''
+        assert self.has_attribute('r'), "Attribute 'r' required."
+        for ipart in range(self.num):
+            # Slice a box from the field around the particle
+            rp = self.attr['r'][ipart]+padding
+            for pos in self.position_with_duplicates(ipart,padding=padding):
+                idxlo = np.array(fld.nearest_gridpoint(pos-rp,lower=True))
+                sfdim = np.ceil(2*rp/fld.spacing+1).astype('int')
+                sf = fld.extract_subfield(idxlo,sfdim,deep=False,strict_bounds=False)
+                xsf = sf.x().reshape((-1,1,1))
+                ysf = sf.y().reshape((1,-1,1))
+                zsf = sf.z().reshape((1,1,-1))
+                dist = (xsf-pos[0])**2 + (ysf-pos[1])**2 + (zsf-pos[2])**2
+                sf.data[dist<=rp*rp] = cval
         return
     def to_vtk(self,deep=False):
         import pyvista as pv