I am trying to write a program that, given one vertex of the tetrahedron circumscribed by the unit sphere, finds the other 3 vertices, their fibers in S^3 and the stereographic projections of these fibers onto R^3. So each of these fibers and their stereographic projections are functions of theta (in fact they are circles), and I want to plot these circles.
Do I have to explicitly write out each function? I am trying to use a list of functions using the 'lambda' keyword, but it's not working as I hoped. I'm new to python, so any help would be great.
My code:
import numpy as np
import matplotlib as mpl
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
from collections import defaultdict
from _future_ import division
# Tetrahedron
# given the first vertex of the tetrahedron, calculate the
# coordinates of the other 3 vertices
v1 = np.array([np.sqrt(2.)/3., -1/3., np.sqrt(2./3.)], dtype = float)
# find two unit vectors orthogonal to v_1
x=np.array([1.,0.,0.])
b1=np.cross(v1,x)
u1=b1/np.linalg.norm(b1)
b2=np.cross(v1,b1)
u2=b2/np.linalg.norm(b2)
# find other vertices
v2=-1/3.*v1+np.sqrt(8)/3.*u1
v3=-1/3.*v1+np.sqrt(2)/3.*(-1.*u1+np.sqrt(3.)*u2)
v4=-1/3.*v1+np.sqrt(2)/3.*(-1.*u1-np.sqrt(3.)*u2)
V=np.array([v1,v2,v3,v4])
# Find fibers of each vertex (S^3)
k=[]
for i in range(4):
a.append(lambda the, z=i:k*(1.+V[z][2])*np.cos(the))
b=[] # 2nd coord.
for i in range(4):
b.append(lambda the, z=i:k*(V[z][0]*np.sin(the)-\
V[z][1]*np.cos(the)))
c=[] # 3rd coord.
for i in range(4):
c.append(lambda the, z=i:k*(V[z][0]*np.cos(the)+\
V[z][1]*np.sin(the)))
d=[] # 4th coord.
for i in range(4):
d.append(lambda the, z=i:k*(np.sin(the)+V[z][2]*np.sin(the)))
# Find stereographic projection of fibers (R^3)
q=[]
for i in range(4):
q.append(lambda the, z=i:b[z]/(1.-a[z]))
r=[]
for i in range(4):
r.append(lambda the, z=i:c[z]/(1.-a[z]))
s=[]
for i in range(4):
s.append(lambda the, z=i:d[z]/(1.-a[z]))
Now, I would like to create a 3d plot. I wrote
fig = plt.figure()
ax = fig.gca(projection='3d')
the = np.linspace(0, 2*np.pi, 100)
q=b[0]/(1-a[0])
r=c[0]/(1-a[0])
s=d[0]/(1-a[0])
ax.plot(q,r,s)
plt.show()
but I get "TypeError: unsupported operand type(s) for -: 'int' and 'function'"
Update: Changed my approach. Here is my new code. I'm still wondering how I can loop over the vertices for plotting so I don't have so much code reuse.
import numpy as np
import matplotlib as mpl
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
from collections import defaultdict
# Tetrahedron
# given the first vertex of the tetrahedron, calculate the
# coordinates of the other 3 vertices
v1 = np.array([np.sqrt(2.)/3., -1/3., np.sqrt(2./3.)], dtype = float)
# find two unit vectors orthogonal to v_1
x=np.array([1.,0.,0.])
b1=np.cross(v1,x)
u1=b1/np.linalg.norm(b1)
b2=np.cross(v1,b1)
u2=b2/np.linalg.norm(b2)
# find other vertices
v2=-1/3.*v1+np.sqrt(8)/3.*u1
v3=-1/3.*v1+np.sqrt(2)/3.*(-1.*u1+np.sqrt(3.)*u2)
v4=-1/3.*v1+np.sqrt(2)/3.*(-1.*u1-np.sqrt(3.)*u2)
V=np.array([v1,v2,v3,v4])
k=[]
for i in range(4):
k.append(1./(2.+2.*V[i][2]))
fig = plt.figure()
ax = fig.gca(projection='3d')
the = np.linspace(0, 2*np.pi, 100)
q0=k[0]*(V[0][0]*np.sin(the)-V[0][1]*np.cos(the))/\
(1.-k[0]*(1.+V[0][2])*np.cos(the))
r0=k[0]*(V[0][0]*np.cos(the)+V[0][1]*np.sin(the))/\
(1.-k[0]*(1.+V[0][2])*np.cos(the))
s0=k[0]*(1.+V[0][2])*np.sin(the)/\
(1.-k[0]*(1.+V[0][2])*np.cos(the))
ax.plot(q0,r0,s0,'b')
q1=k[1]*(V[1][0]*np.sin(the)-V[1][1]*np.cos(the))/\
(1.-k[1]*(1.+V[1][2])*np.cos(the))
r1=k[1]*(V[1][0]*np.cos(the)+V[1][1]*np.sin(the))/\
(1.-k[1]*(1.+V[1][2])*np.cos(the))
s1=k[1]*(1.+V[1][2])*np.sin(the)/\
(1.-k[1]*(1.+V[1][2])*np.cos(the))
ax.plot(q1,r1,s1,'r')
q2=k[2]*(V[2][0]*np.sin(the)-V[2][1]*np.cos(the))/\
(1.-k[2]*(1.+V[2][2])*np.cos(the))
r2=k[2]*(V[2][0]*np.cos(the)+V[2][1]*np.sin(the))/\
(1.-k[2]*(1.+V[2][2])*np.cos(the))
s2=k[2]*(1.+V[2][2])*np.sin(the)/\
(1.-k[2]*(1.+V[2][2])*np.cos(the))
ax.plot(q2,r2,s2,'g')
q3=k[3]*(V[3][0]*np.sin(the)-V[3][1]*np.cos(the))/\
(1.-k[3]*(1.+V[3][2])*np.cos(the))
r3=k[3]*(V[3][0]*np.cos(the)+V[3][1]*np.sin(the))/\
(1.-k[3]*(1.+V[3][2])*np.cos(the))
s3=k[3]*(1.+V[3][2])*np.sin(the)/\
(1.-k[3]*(1.+V[3][2])*np.cos(the))
ax.plot(q3,r3,s3,'m')
plt.show()
You can loop over the vertices by zipping the k and V lists. If you want a legend or a different color cycle, those would have to be added to the call to zip.
for k0, V0 in zip(k, V):
q=k0*(V0[0]*np.sin(the)-V0[1]*np.cos(the))/\
(1.-k0*(1.+V0[2])*np.cos(the))
r=k0*(V0[0]*np.cos(the)+V0[1]*np.sin(the))/\
(1.-k0*(1.+V0[2])*np.cos(the))
s=k0*(1.+V0[2])*np.sin(the)/\
(1.-k0*(1.+V0[2])*np.cos(the))
ax.plot(q,r,s,)