Tesseroid has split method

cookbook/mesher_tesseroid.py

@@ -0,0 +1,14 @@
+"""
+Meshing: Make and plot a tesseroid with the Earth
+"""
+from fatiando import mesher
+from fatiando.vis import myv
+
+model = mesher.Tesseroid(-10, 10, 50, 60, 500000, 0)
+myv.figure(zdown=False)
+myv.tesseroids([model])
+myv.earth()
+myv.continents()
+myv.meridians(range(-180, 180, 10))
+myv.parallels(range(-90, 90, 10))
+myv.show()

fatiando/gravmag/tesseroid.py

@@ -3,10 +3,8 @@
 """
 import numpy
 
-from fatiando.mesher import Tesseroid
 from fatiando.constants import SI2MGAL, SI2EOTVOS, MEAN_EARTH_RADIUS, G
 
-
 try:
     from fatiando.gravmag import _ctesseroid as _kernels
 except ImportError:
@@ -114,7 +112,7 @@ def _optimal_discretize(tesseroids, lons, lats, heights, kernel, ratio, dens):
     radii = MEAN_EARTH_RADIUS + heights
     # Start the computations
     result = numpy.zeros(ndata, numpy.float)
-    maxsize = 10000
+    #maxsize = 10000
     for tesseroid in tesseroids:
         if (tesseroid is None or
             ('density' not in tesseroid.props and dens is None)):
@@ -138,26 +136,14 @@ def _optimal_discretize(tesseroids, lons, lats, heights, kernel, ratio, dens):
                 #    log.warning("Maximum LIFO size reached")
                 #    dont_divide.extend(need_divide)
                 #else:
-                #    lifo.extend([need_divide, t] for t in _split(tess))
-                lifo.extend([need_divide, t] for t in _split(tess))
+                #    lifo.extend([need_divide, t] for t in tess.split(2, 2, 2)
+                lifo.extend([need_divide, t] for t in tess.split(2, 2, 2))
             if len(dont_divide):
                 result[dont_divide] += G*density*kernel(
                     tess, rlons[dont_divide], rlats[dont_divide],
                     radii[dont_divide], _glq_nodes, _glq_weights)
     return result
 
-def _split(tesseroid):
-    dlon = 0.5*(tesseroid.e - tesseroid.w)
-    dlat = 0.5*(tesseroid.n - tesseroid.s)
-    dh = 0.5*(tesseroid.top - tesseroid.bottom)
-    wests = [tesseroid.w, tesseroid.w + dlon]
-    souths = [tesseroid.s, tesseroid.s + dlat]
-    bottoms = [tesseroid.bottom, tesseroid.bottom + dh]
-    split = [
-        Tesseroid(i, i + dlon, j, j + dlat, k + dh, k, props=tesseroid.props)
-        for i in wests for j in souths for k in bottoms]
-    return split
-
 def _distance(tesseroid, lon, lat, radius, points):
     lons = lon[points]
     lats = lat[points]

fatiando/mesher.py

@@ -486,6 +486,58 @@ def get_bounds(self):
         """
         return [self.w, self.e, self.s, self.n, self.top, self.bottom]
 
+    def split(self, nlon, nlat, nh):
+        """
+        Split the tesseroid into smaller ones.
+
+        The smaller tesseroids will share the large one's props.
+
+        Parameters:
+
+        * nlon, nlat, nh : int
+            The number of sections to split in the longitudinal, latitudinal,
+            and vertical dimensions
+
+        Returns:
+
+        * tesseroids : list
+            A list of nlon*nlat*nh tesseroids that make up the larger one.
+
+
+        Examples::
+
+        >>> tess = Tesseroid(-10, 10, -20, 20, 0, -40, {'density':2})
+        >>> split = tess.split(1, 2, 2)
+        >>> print len(split)
+        4
+        >>> for t in split:
+        ...     print t
+        w:-10 | e:10 | s:-20 | n:0 | top:-20 | bottom:-40 | density:2
+        w:-10 | e:10 | s:-20 | n:0 | top:0 | bottom:-20 | density:2
+        w:-10 | e:10 | s:0 | n:20 | top:-20 | bottom:-40 | density:2
+        w:-10 | e:10 | s:0 | n:20 | top:0 | bottom:-20 | density:2
+        >>> tess = Tesseroid(-15, 15, -20, 20, 0, -40)
+        >>> split = tess.split(3, 1, 1)
+        >>> print len(split)
+        3
+        >>> for t in split:
+        ...     print t
+        w:-15 | e:-5 | s:-20 | n:20 | top:0 | bottom:-40
+        w:-5 | e:5 | s:-20 | n:20 | top:0 | bottom:-40
+        w:5 | e:15 | s:-20 | n:20 | top:0 | bottom:-40
+
+        """
+        wests = numpy.linspace(self.w, self.e, nlon + 1)
+        souths = numpy.linspace(self.s, self.n, nlat + 1)
+        bottoms = numpy.linspace(self.bottom, self.top, nh + 1)
+        dlon = wests[1] - wests[0]
+        dlat = souths[1] - souths[0]
+        dh = bottoms[1] - bottoms[0]
+        tesseroids = [
+            Tesseroid(i, i + dlon, j, j + dlat, k + dh, k, props=self.props)
+            for i in wests[:-1] for j in souths[:-1] for k in bottoms[:-1]]
+        return tesseroids
+
 class Sphere(GeometricElement):
     """
     Create a sphere.