Python scripts demonstrating half vector isocontour

halfvector.py

import numpy as np
import matplotlib as mpl
from mpl_toolkits.mplot3d import Axes3D
import matplotlib.pyplot as plt

def reflect(v, n):
    return 2 * np.dot(v, n) * n - v

def transform(beta, theta):
    l = np.array([-np.cos(beta), np.sin(beta), 0.0])
    n = np.array([0.0, 1.0, 0.0])
    r = reflect(l, n)
    
    d = np.sin(theta)
    proj = np.cos(theta) * n
    def f(t):
        p = np.array([np.cos(t), 0.0, np.sin(t)]) * d
        h = p + proj
        return reflect(l, h)
    
    ts = np.linspace(0.0, 2 * np.pi, 64)

    hs = map(f, ts)

    q = np.pi / 2 - beta
    m = np.matrix([
        [np.cos(q), -np.sin(q), 0],
        [np.sin(q), np.cos(q), 0],
        [0, 0, 1]
    ])
    
    def topview(c):
        tc = m * (np.matrix(c).transpose())
        return np.array([tc.item(0), tc.item(1), tc.item(2)])
        # return c

    ps = map(topview, hs)

    return np.array(list(ps))

def plot2(ps, ax):
    # ax.set_xlabel("X")
    # ax.set_ylabel("Z")
    ax.set_xlim(-0.5, 0.5)
    ax.set_ylim(-0.5, 0.5)
    x = ps[:, 0]
    z = ps[:, 2]
    ax.plot(x, z)

def plot3(ps):
    fig = plt.figure()
    ax = fig.gca(projection='3d')
    ax.set_xlabel("X")
    ax.set_ylabel("Y")
    ax.set_zlabel("Z")
    ax.set_xlim(-1.5, 1.5)
    ax.set_ylim(-1.5, 1.5)
    ax.set_zlim(-1.5, 1.5)
    ax.view_init(0, 90)
    x = ps[:, 0]
    y = ps[:, 1]
    z = ps[:, 2]
    ax.plot(x, y, z)
    plt.show()

def demo():
    plt.figure(1)
    
    ax = plt.subplot(221, aspect='equal')
    plot2(transform(np.pi / 3, np.pi / 16), ax)
    ax.set_title('pi / 3')

    ax = plt.subplot(222, aspect='equal')
    plot2(transform(np.pi / 4, np.pi / 16), ax)
    ax.set_title('pi / 4')

    ax = plt.subplot(223, aspect='equal')
    plot2(transform(np.pi / 6, np.pi / 16), ax)
    ax.set_title('pi / 6')

    ax = plt.subplot(224, aspect='equal')
    plot2(transform(np.pi / 12, np.pi / 16), ax)
    ax.set_title('pi / 12')

    plt.show()