-
Notifications
You must be signed in to change notification settings - Fork 1
/
Copy pathsphere-11.py
115 lines (104 loc) · 4.94 KB
/
sphere-11.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
from sympy import Eq, Matrix, cancel, expand, fraction, gcd_list, lcm_list, poly, solve, symbols
def sphere(P1, P2, P3, P4):
# return F(x, y, z) such that F(x, y, z) = 0 is the sphere's equation
g, h, j, k, x, y, z = symbols('g, h, j, k, x, y, z')
sphere_eq = Eq(x**2 + y**2 + z**2 + g*x + h*y + j*z + k, 0)
sphere_eqs = []
sphere_eqs.append(sphere_eq.subs(x, P1[0]).subs(y, P1[1]).subs(z, P1[2]))
sphere_eqs.append(sphere_eq.subs(x, P2[0]).subs(y, P2[1]).subs(z, P2[2]))
sphere_eqs.append(sphere_eq.subs(x, P3[0]).subs(y, P3[1]).subs(z, P3[2]))
sphere_eqs.append(sphere_eq.subs(x, P4[0]).subs(y, P4[1]).subs(z, P4[2]))
s = solve(sphere_eqs, (g, h, j, k))
return fraction(cancel(x**2 + y**2 + z**2 + s[g]*x + s[h]*y + s[j]*z + s[k]))[0]
def coplanar(P1, P2, P3, P4):
mat = []
mat.append([P1[0], P1[1], P1[2], 1])
mat.append([P2[0], P2[1], P2[2], 1])
mat.append([P3[0], P3[1], P3[2], 1])
mat.append([P4[0], P4[1], P4[2], 1])
return Eq(fraction(cancel(Matrix(mat).det()))[0], 0)
def plane(P1, P2, P3):
x, y, z = symbols('x, y, z')
return coplanar((x, y, z), P1, P2, P3)
def intersect(P1, P2, P3, slope, sphere_l, ai0, ai1, ai2):
axes = symbols('x, y, z')
axis0, axis1, axis2 = axes[ai0], axes[ai1], axes[ai2]
# return the intersect point P of:
# `axis2 = plane_p(axis0, axis1)` passing through P1, P2 and P3
# `axis1 = plane_q(axis0)` with `slope = axis1/axis0` passing through P1
# `sphere_l(x, y, z) = 0`
plane_p = solve(plane(P1, P2, P3), axis2)[0]
plane_q = slope*(axis0 - P1[ai0]) + P1[ai1]
sphere_coeffs = poly(sphere_l.subs(axis2, plane_p).subs(axis1, plane_q), axis0).all_coeffs()
P = {}
P[axis0] = -sphere_coeffs[1]/sphere_coeffs[0] - P1[ai0]
P[axis1] = plane_q.subs(axis0, P[axis0])
P[axis2] = plane_p.subs(axis0, P[axis0]).subs(axis1, P[axis1])
return P[axes[0]], P[axes[1]], P[axes[2]]
def multiplied(x, y, z, w):
x1, y1, z1, w1 = fraction(cancel(x)), fraction(cancel(y)), fraction(cancel(z)), fraction(cancel(w))
lcd = lcm_list([x1[1], y1[1], z1[1], w1[1]])
return x1[0]*cancel(lcd/x1[1]), y1[0]*cancel(lcd/y1[1]), z1[0]*cancel(lcd/z1[1]), w1[0]*cancel(lcd/w1[1])
def to_homogeneous(P):
return multiplied(P[0], P[1], P[2], 1)
def reduced(x, y, z, w):
gcd = gcd_list([x, y, z, w])
return cancel(x/gcd), cancel(y/gcd), cancel(z/gcd), cancel(w/gcd)
def cross(P1, P2, P3):
x10, x11, x12, x13 = P1
x20, x21, x22, x23 = P2
x30, x31, x32, x33 = P3
# generated by cross-3d.py
x = -x11*x22*x33 + x11*x23*x32 + x12*x21*x33 - x12*x23*x31 - x13*x21*x32 + x13*x22*x31
y = x10*x22*x33 - x10*x23*x32 - x12*x20*x33 + x12*x23*x30 + x13*x20*x32 - x13*x22*x30
z = -x10*x21*x33 + x10*x23*x31 + x11*x20*x33 - x11*x23*x30 - x13*x20*x31 + x13*x21*x30
w = x10*x21*x32 - x10*x22*x31 - x11*x20*x32 + x11*x22*x30 + x12*x20*x31 - x12*x21*x30
return reduced(x, y, z, w)
def on_sphere(sphere_h, P):
x, y, z, w = symbols('x, y, z, w')
return expand(sphere_h.subs(x, P[0]).subs(y, P[1]).subs(z, P[2]).subs(w, P[3])) == 0
def main():
# https://imomath.com/index.cgi?page=inversion (Problem 11)
a, b, c, d, e, f, g, h, j, x, y, z, w = symbols('a, b, c, d, e, f, g, h, j, x, y, z, w')
# quick test a special case
# b, d, e = 0, 0, 0
O, A, B, C = (0, 0, 0), (a, 0, 0), (b, c, 0), (d, e, f)
sphere_l = sphere(O, A, B, C)
print('Sphere Equation:', sphere_l, '= 0')
# If we choose plane $y=kx$ through point O, we'll get G with too many terms:
# G[0] has 12319 terms
# G[1] has 9959 terms
# G[2] has 6337 terms
# G[3] has 8760 terms
# D = intersect(O, B, C, g, sphere_l, 1, 2, 0)
# E = intersect(O, C, A, h, sphere_l, 2, 0, 1)
# F = intersect(O, A, B, j, sphere_l, 0, 1, 2)
D = intersect(C, B, O, g, sphere_l, 1, 2, 0)
E = intersect(A, C, O, h, sphere_l, 2, 0, 1)
F = intersect(B, A, O, j, sphere_l, 0, 1, 2)
sphere_h = poly(sphere_l, (x, y, z)).homogenize(w).expr
print('Sphere Equation in Homogeneous:', sphere_h, '= 0')
A, B, C = to_homogeneous(A), to_homogeneous(B), to_homogeneous(C)
D, E, F = to_homogeneous(D), to_homogeneous(E), to_homogeneous(F)
print('D:', D)
print('E:', E)
print('F:', F)
print('Is A on Sphere?', on_sphere(sphere_h, A))
print('Is B on Sphere?', on_sphere(sphere_h, B))
print('Is C on Sphere?', on_sphere(sphere_h, C))
print('Is D on Sphere?', on_sphere(sphere_h, D))
print('Is E on Sphere?', on_sphere(sphere_h, E))
print('Is F on Sphere?', on_sphere(sphere_h, F))
CDE, BFD, AEF = cross(C, D, E), cross(B, F, D), cross(A, E, F)
print('CDE:', CDE)
print('BFD:', BFD)
print('AEF:', AEF)
G = cross(CDE, BFD, AEF)
print('G:', G)
print('G[0] has', len(G[0].args), 'terms')
print('G[1] has', len(G[1].args), 'terms')
print('G[2] has', len(G[2].args), 'terms')
print('G[3] has', len(G[3].args), 'terms')
print('Is G on Sphere?', on_sphere(sphere_h, G))
if __name__ == '__main__':
main()