import numpy as np
    import matplotlib.pyplot as plt

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

    normals = np.vstack((np.eye(3), -np.eye(3)))

    theta = np.linspace(0, pi, 181)
    phi   = np.linspace(-pi, pi, 361)

    Phi, Theta = np.meshgrid(phi, theta)

    (cp, sp), (ct, st) = [[f(thing) for f in (np.cos, np.sin)] for thing in (Phi, Theta)]

    directions = np.stack([cp*st, sp*st, ct])

    cos_thetas = [np.maximum((directions * normal[:, None, None]).sum(axis=0), 0) for normal in normals]

    if True:
        plt.figure()
        for i, (cth, normal) in enumerate(zip(cos_thetas, normals)):
            plt.subplot(4, 3, i+1)
            plt.imshow(cth)
            # plt.colorbar()
            plt.title(str(normal.astype(int)), fontsize=14)
        cth = sum(cos_thetas)
        plt.subplot(2, 1, 2)
        plt.imshow(cth/cth.max(), vmin=0)
        plt.colorbar()
        plt.xlabel('phi (longitude)', fontsize=14)
        plt.ylabel('theta (latitude)', fontsize=14)
        plt.show()