def dotem(a, b):
        return (a*b).sum(axis=0)

    def get_B2D(x):

        r    = (x - x0)/R0
        rabs = np.sqrt((r**2).sum(axis=0)) # [..., None]
        rhat = r / rabs

        B = B0 * (3*dotem(p, rhat) * rhat - p) / rabs**3

        return B

    import numpy as np
    import matplotlib.pyplot as plt

    halfpi, pi, twopi, fourpi = [f*np.pi for f in (0.5, 1, 2, 4)]
    degs, rads                = 180/pi, pi/180

    mu0  = fourpi * 1E-07
    R0   = 6371.     # km
    B0   = 3.12E-5   # T
    inc  = rads * 11.0

    p    = np.array([f(inc) for f in (np.sin, np.zeros_like, np.cos)], dtype=float)[:, None, None]
    x0   = np.array([-500, 0, 0], dtype=float)[:, None, None]  # offset (manifests as South Atlantic Anomaly)

    pp   = p[:, 0, 0]
    xx0  = x0[:, 0, 0]

    R    = 6778.

    lat0, lat1 = rads*(90-50), rads*(90+50)
    lat      = np.linspace( lat0,  lat1, 301) # 0 to 180 degrees (spherical coordinates)
    lon      = np.linspace(-pi, pi, 601)
    LAT, LON = np.meshgrid(lat, lon, indexing='ij')
    (clat, clon), (slat, slon) = [[f(x) for x in (LAT, LON)] for f in (np.cos, np.sin)]

    XYZ = np.stack([R * x for x in (slat*clon, slat*slon, clat)])

    B = get_B2D(XYZ)
    Bmag = np.sqrt((B**2).sum(axis=0))

    if True:
        plt.figure()
        plt.imshow(Bmag)
        plt.colorbar()
        plt.show()