Skip to content

Instantly share code, notes, and snippets.

Embed
What would you like to do?
Python scripts demonstrating half vector isocontour
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()
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
You can’t perform that action at this time.