docs: def plotMT(T, N, P, size=200, plot_zerotrace=True,
             x0=0, y0=0, xy=(0, 0), width=200):
      """
      Uses a principal axis T, N and P to draw a beach ball plot.
  
      :param ax: axis object of a matplotlib figure
      :param T: :class:`~PrincipalAxis`
      :param N: :class:`~PrincipalAxis`
      :param P: :class:`~PrincipalAxis`
  
      Adapted from ps_tensor / utilmeca.c / `Generic Mapping Tools (GMT)`_.
  
      .. _`Generic Mapping Tools (GMT)`: http://gmt.soest.hawaii.edu
      """
      # check if one or two widths are specified (Circle or Ellipse)
      try:
          assert(len(width) == 2)
      except TypeError:
          width = (width, width)
      collect = []
      colors = []
      res = [value / float(size) for value in width]
      b = 1
      big_iso = 0
      j = 1
      j2 = 0
      j3 = 0
      n = 0
      azi = np.zeros((3, 2))
      x = np.zeros(400)
      y = np.zeros(400)
      x2 = np.zeros(400)
      y2 = np.zeros(400)
      x3 = np.zeros(400)
      y3 = np.zeros(400)
      xp1 = np.zeros(800)
      yp1 = np.zeros(800)
      xp2 = np.zeros(400)
      yp2 = np.zeros(400)
  
      a = np.zeros(3)
      p = np.zeros(3)
      v = np.zeros(3)
      a[0] = T.strike
      a[1] = N.strike
      a[2] = P.strike
      p[0] = T.dip
      p[1] = N.dip
      p[2] = P.dip
      v[0] = T.val
      v[1] = N.val
      v[2] = P.val
  
      vi = (v[0] + v[1] + v[2]) / 3.
      for i in range(0, 3):
          v[i] = v[i] - vi
  
      radius_size = size * 0.5
  
      if np.fabs(v[0] * v[0] + v[1] * v[1] + v[2] * v[2]) < EPSILON:
          # pure implosion-explosion
          if vi > 0.:
              cir = patches.Ellipse(xy, width=width[0], height=width[1])
              collect.append(cir)
              colors.append('b')
          if vi < 0.:
              cir = patches.Ellipse(xy, width=width[0], height=width[1])
              collect.append(cir)
              colors.append('w')
          return colors, collect
  
      if np.fabs(v[0]) >= np.fabs(v[2]):
          d = 0
          m = 2
      else:
          d = 2
          m = 0
  
      if (plot_zerotrace):
          vi = 0.
  
      f = -v[1] / float(v[d])
      iso = vi / float(v[d])
  
      # Cliff Frohlich, Seismological Research letters,
      # Vol 7, Number 1, January-February, 1996
      # Unless the isotropic parameter lies in the range
      # between -1 and 1 - f there will be no nodes whatsoever
  
      if iso < -1:
          cir = patches.Ellipse(xy, width=width[0], height=width[1])
          collect.append(cir)
          colors.append('w')
          return colors, collect
      elif iso > 1 - f:
          cir = patches.Ellipse(xy, width=width[0], height=width[1])
          collect.append(cir)
          colors.append('b')
          return colors, collect
  
      spd = np.sin(p[d] * D2R)
      cpd = np.cos(p[d] * D2R)
      spb = np.sin(p[b] * D2R)
      cpb = np.cos(p[b] * D2R)
      spm = np.sin(p[m] * D2R)
      cpm = np.cos(p[m] * D2R)
      sad = np.sin(a[d] * D2R)
      cad = np.cos(a[d] * D2R)
      sab = np.sin(a[b] * D2R)
      cab = np.cos(a[b] * D2R)
      sam = np.sin(a[m] * D2R)
      cam = np.cos(a[m] * D2R)
  
      for i in range(0, 360):
          fir = i * D2R
          s2alphan = (2. + 2. * iso) / \
              float(3. + (1. - 2. * f) * np.cos(2. * fir))
          if s2alphan > 1.:
              big_iso += 1
          else:
              alphan = np.arcsin(np.sqrt(s2alphan))
              sfi = np.sin(fir)
              cfi = np.cos(fir)
              san = np.sin(alphan)
              can = np.cos(alphan)
  
              xz = can * spd + san * sfi * spb + san * cfi * spm
              xn = can * cpd * cad + san * sfi * cpb * cab + \
                  san * cfi * cpm * cam
              xe = can * cpd * sad + san * sfi * cpb * sab + \
                  san * cfi * cpm * sam
  
              if np.fabs(xn) < EPSILON and np.fabs(xe) < EPSILON:
                  takeoff = 0.
                  az = 0.
              else:
                  az = np.arctan2(xe, xn)
                  if az < 0.:
                      az += np.pi * 2.
                  takeoff = np.arccos(xz / float(np.sqrt(xz * xz + xn * xn +
                                                         xe * xe)))
              if takeoff > np.pi / 2.:
                  takeoff = np.pi - takeoff
                  az += np.pi
                  if az > np.pi * 2.:
                      az -= np.pi * 2.
              r = np.sqrt(2) * np.sin(takeoff / 2.)
              si = np.sin(az)
              co = np.cos(az)
              if i == 0:
                  azi[i][0] = az
                  x[i] = x0 + radius_size * r * si
                  y[i] = y0 + radius_size * r * co
                  azp = az
              else:
                  if np.fabs(np.fabs(az - azp) - np.pi) < D2R * 10.:
                      azi[n][1] = azp
                      n += 1
                      azi[n][0] = az
                  if np.fabs(np.fabs(az - azp) - np.pi * 2.) < D2R * 2.:
                      if azp < az:
                          azi[n][0] += np.pi * 2.
                      else:
                          azi[n][0] -= np.pi * 2.
                  if n == 0:
                      x[j] = x0 + radius_size * r * si
                      y[j] = y0 + radius_size * r * co
                      j += 1
                  elif n == 1:
                      x2[j2] = x0 + radius_size * r * si
                      y2[j2] = y0 + radius_size * r * co
                      j2 += 1
                  elif n == 2:
                      x3[j3] = x0 + radius_size * r * si
                      y3[j3] = y0 + radius_size * r * co
                      j3 += 1
                  azp = az
      azi[n][1] = az
  
      if v[1] < 0.:
          rgb1 = 'b'
          rgb2 = 'w'
      else:
          rgb1 = 'w'
          rgb2 = 'b'
  
      cir = patches.Ellipse(xy, width=width[0], height=width[1])
      collect.append(cir)
      colors.append(rgb2)
      if n == 0:
          collect.append(xy2patch(x[0:360], y[0:360], res, xy))
          colors.append(rgb1)
          return colors, collect
      elif n == 1:
          for i in range(0, j):
              xp1[i] = x[i]
              yp1[i] = y[i]
          if azi[0][0] - azi[0][1] > np.pi:
              azi[0][0] -= np.pi * 2.
          elif azi[0][1] - azi[0][0] > np.pi:
              azi[0][0] += np.pi * 2.
          if azi[0][0] < azi[0][1]:
              az = azi[0][1] - D2R
              while az > azi[0][0]:
                  si = np.sin(az)
                  co = np.cos(az)
                  xp1[i] = x0 + radius_size * si
                  yp1[i] = y0 + radius_size * co
                  i += 1
                  az -= D2R
          else:
              az = azi[0][1] + D2R
              while az < azi[0][0]:
                  si = np.sin(az)
                  co = np.cos(az)
                  xp1[i] = x0 + radius_size * si
                  yp1[i] = y0 + radius_size * co
                  i += 1
                  az += D2R
          collect.append(xy2patch(xp1[0:i], yp1[0:i], res, xy))
          colors.append(rgb1)
          for i in range(0, j2):
              xp2[i] = x2[i]
              yp2[i] = y2[i]
          if azi[1][0] - azi[1][1] > np.pi:
              azi[1][0] -= np.pi * 2.
          elif azi[1][1] - azi[1][0] > np.pi:
              azi[1][0] += np.pi * 2.
          if azi[1][0] < azi[1][1]:
              az = azi[1][1] - D2R
              while az > azi[1][0]:
                  si = np.sin(az)
                  co = np.cos(az)
                  xp2[i] = x0 + radius_size * si
                  i += 1
                  yp2[i] = y0 + radius_size * co
                  az -= D2R
          else:
              az = azi[1][1] + D2R
              while az < azi[1][0]:
                  si = np.sin(az)
                  co = np.cos(az)
                  xp2[i] = x0 + radius_size * si
                  i += 1
                  yp2[i] = y0 + radius_size * co
                  az += D2R
          collect.append(xy2patch(xp2[0:i], yp2[0:i], res, xy))
          colors.append(rgb1)
          return colors, collect
      elif n == 2:
          for i in range(0, j3):
              xp1[i] = x3[i]
              yp1[i] = y3[i]
          for ii in range(0, j):
              xp1[i] = x[ii]
              i += 1
              yp1[i] = y[ii]
          if big_iso:
              ii = j2 - 1
              while ii >= 0:
                  xp1[i] = x2[ii]
                  i += 1
                  yp1[i] = y2[ii]
                  ii -= 1
              collect.append(xy2patch(xp1[0:i], yp1[0:i], res, xy))
              colors.append(rgb1)
              return colors, collect
  
          if azi[2][0] - azi[0][1] > np.pi:
              azi[2][0] -= np.pi * 2.
          elif azi[0][1] - azi[2][0] > np.pi:
              azi[2][0] += np.pi * 2.
          if azi[2][0] < azi[0][1]:
              az = azi[0][1] - D2R
              while az > azi[2][0]:
                  si = np.sin(az)
                  co = np.cos(az)
                  xp1[i] = x0 + radius_size * si
                  i += 1
                  yp1[i] = y0 + radius_size * co
                  az -= D2R
          else:
              az = azi[0][1] + D2R
              while az < azi[2][0]:
                  si = np.sin(az)
                  co = np.cos(az)
                  xp1[i] = x0 + radius_size * si
                  i += 1
                  yp1[i] = y0 + radius_size * co
                  az += D2R
          collect.append(xy2patch(xp1[0:i], yp1[0:i], res, xy))
          colors.append(rgb1)
  
          for i in range(0, j2):
              xp2[i] = x2[i]
              yp2[i] = y2[i]
          if azi[1][0] - azi[1][1] > np.pi:
              azi[1][0] -= np.pi * 2.
          elif azi[1][1] - azi[1][0] > np.pi:
              azi[1][0] += np.pi * 2.
          if azi[1][0] < azi[1][1]:
              az = azi[1][1] - D2R
              while az > azi[1][0]:
                  si = np.sin(az)
                  co = np.cos(az)
                  xp2[i] = x0 + radius_size * si
                  i += 1
                  yp2[i] = y0 + radius_size * co
                  az -= D2R
          else:
              az = azi[1][1] + D2R
              while az < azi[1][0]:
                  si = np.sin(az)
                  co = np.cos(az)
                  xp2[i] = x0 + radius_size * si
                  i += 1
                  yp2[i] = y0 + radius_size * co
                  az += D2R
          collect.append(xy2patch(xp2[0:i], yp2[0:i], res, xy))
          colors.append(rgb1)
          return colors, collect