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