#!/usr/bin/python
import sys
import os,sympy,time
from sympy.galgebra.GA import set_main, make_symbols, types, MV, ZERO, ONE, HALF
from sympy import collect
set_main(sys.modules[__name__])
def F(x):
"""
Conformal Mapping Function
"""
Fx = HALF*((x*x)*n+2*x-nbar)
return(Fx)
def make_vector(a,n = 3):
if type(a) == types.StringType:
sym_str = ''
for i in range(n):
sym_str += a+str(i)+' '
sym_lst = make_symbols(sym_str)
sym_lst.append(ZERO)
sym_lst.append(ZERO)
a = MV(sym_lst,'vector')
return(F(a))
if __name__ == '__main__':
ti = time.time()
MV.setup('a b c d e')
MV.set_str_format(1)
print 'e|(a^b) =',e|(a^b)
print 'e|(a^b^c) =',e|(a^b^c)
print 'a*(b^c)-b*(a^c)+c*(a^b) =',a*(b^c)-b*(a^c)+c*(a^b)
print 'e|(a^b^c^d) =',e|(a^b^c^d)
print -d*(a^b^c)+c*(a^b^d)-b*(a^c^d)+a*(b^c^d)
print (a^b)|(c^d)
# FIXME: currently broken
"""
print 'Example: non-euclidian distance calculation'
metric = '0 # #,# 0 #,# # 1'
MV.setup('X Y e',metric)
MV.set_str_format(1)
L = X^Y^e
B = L*e
Bsq = (B*B)()
print 'L = X^Y^e is a non-euclidian line'
print 'B = L*e =',B
BeBr =B*e*B.rev()
print 'B*e*B.rev() =',BeBr
print 'B^2 =',Bsq
print 'L^2 =',(L*L)()
make_symbols('s c Binv M S C alpha')
Bhat = Binv*B # Normalize translation generator
R = c+s*Bhat # Rotor R = exp(alpha*Bhat/2)
print 's = sinh(alpha/2) and c = cosh(alpha/2)'
print 'R = exp(alpha*B/(2*|B|)) =',R
Z = R*X*R.rev()
Z.expand()
Z.collect([Binv,s,c,XdotY])
print 'R*X*R.rev() =',Z
W = Z|Y
W.expand()
W.collect([s*Binv])
print '(R*X*rev(R)).Y =',W
M = 1/Bsq
W.subs(Binv**2,M)
W.simplify()
Bmag = sympy.sqrt(XdotY**2-2*XdotY*Xdote*Ydote)
W.collect([Binv*c*s,XdotY])
W.subs(2*XdotY**2-4*XdotY*Xdote*Ydote,2/(Binv**2))
W.subs(2*c*s,S)
W.subs(c**2,(C+1)/2)
W.subs(s**2,(C-1)/2)
W.simplify()
W.subs(1/Binv,Bmag)
W = W().expand()
print '(R*X*R.rev()).Y =',W
nl = '\n'
Wd = collect(W,[C,S],exact=True,evaluate=False)
print 'Wd =',Wd
Wd_1 = Wd[ONE]
Wd_C = Wd[C]
Wd_S = Wd[S]
print '|B| =',Bmag
Wd_1 = Wd_1.subs(Bmag,1/Binv)
Wd_C = Wd_C.subs(Bmag,1/Binv)
Wd_S = Wd_S.subs(Bmag,1/Binv)
print 'Wd[ONE] =',Wd_1
print 'Wd[C] =',Wd_C
print 'Wd[S] =',Wd_S
lhs = Wd_1+Wd_C*C
rhs = -Wd_S*S
lhs = lhs**2
rhs = rhs**2
W = (lhs-rhs).expand()
W = (W.subs(1/Binv**2,Bmag**2)).expand()
print 'W =',W
W = (W.subs(S**2,C**2-1)).expand()
print 'W =',W
W = collect(W,[C,C**2],evaluate=False)
print 'W =',W
a = W[C**2]
b = W[C]
c = W[ONE]
print 'a =',a
print 'b =',b
print 'c =',c
D = (b**2-4*a*c).expand()
print 'Setting to 0 and solving for C gives:'
print 'Discriminant D = b^2-4*a*c =',D
C = (-b/(2*a)).expand()
print 'C = cosh(alpha) = -b/(2*a) =',C
"""
print '\nExample: Conformal representations of circles, lines, spheres, and planes'
metric = '1 0 0 0 0,0 1 0 0 0,0 0 1 0 0,0 0 0 0 2,0 0 0 2 0'
MV.setup('e0 e1 e2 n nbar',metric,debug=0)
MV.set_str_format(1)
e = n+nbar
#conformal representation of points
A = make_vector(e0) # point a = (1,0,0) A = F(a)
B = make_vector(e1) # point b = (0,1,0) B = F(b)
C = make_vector(-1*e0) # point c = (-1,0,0) C = F(c)
D = make_vector(e2) # point d = (0,0,1) D = F(d)
X = make_vector('x',3)
print 'a = e0, b = e1, c = -e0, and d = e2'
print 'A = F(a) = 1/2*(a*a*n+2*a-nbar), etc.'
print 'Circle through a, b, and c'
print 'Circle: A^B^C^X = 0 =',(A^B^C^X)
print 'Line through a and b'
print 'Line : A^B^n^X = 0 =',(A^B^n^X)
print 'Sphere through a, b, c, and d'
print 'Sphere: A^B^C^D^X = 0 =',(A^B^C^D^X)
print 'Plane through a, b, and d'
print 'Plane : A^B^n^D^X = 0 =',(A^B^n^D^X)
L = (A^B^e)^X
print 'Hyperbolic Circle: (A^B^e)^X = 0 =',L
#MV.LaTeX()
metric = '# # # 0 0,'+ \
'# # # 0 0,'+ \
'# # # 0 0,'+ \
'0 0 0 0 2,'+ \
'0 0 0 2 0'
MV.setup('p1 p2 p3 n nbar',metric,debug=0)
MV.set_str_format(1)
P1 = F(p1)
P2 = F(p2)
P3 = F(p3)
print '\nExtracting direction of line from L = P1^P2^n'
L = P1^P2^n
delta = (L|n)|nbar
print '(L.n).nbar=',delta
print '\nExtracting plane of circle from C = P1^P2^P3'
C = P1^P2^P3
delta = ((C^n)|n)|nbar
print '((C^n).n).nbar=',delta
print '(p2-p1)^(p3-p1)=',(p2-p1)^(p3-p1)
metric = '1 # #,'+ \
'# 1 #,'+ \
'# # 1,'
MV.setup('e1 e2 e3',metric)
print 'Example: Reciprocal Frames e1, e2, and e3 unit vectors.\n\n'
E = e1^e2^e3
Esq = (E*E)()
print 'E =',E
print 'E^2 =',Esq
Esq_inv = 1/Esq
E1 = (e2^e3)*E
E2 = (-1)*(e1^e3)*E
E3 = (e1^e2)*E
print 'E1 = (e2^e3)*E =',E1
print 'E2 =-(e1^e3)*E =',E2
print 'E3 = (e1^e2)*E =',E3
w = (E1|e2)
w.collect(MV.g)
w = w().expand()
print 'E1|e2 =',w
w = (E1|e3)
w.collect(MV.g)
w = w().expand()
print 'E1|e3 =',w
w = (E2|e1)
w.collect(MV.g)
w = w().expand()
print 'E2|e1 =',w
w = (E2|e3)
w.collect(MV.g)
w = w().expand()
print 'E2|e3 =',w
w = (E3|e1)
w.collect(MV.g)
w = w().expand()
print 'E3|e1 =',w
w = (E3|e2)
w.collect(MV.g)
w = w().expand()
print 'E3|e2 =',w
w = (E1|e1)
w = w().expand()
Esq = Esq.expand()
print '(E1|e1)/E^2 =',w/Esq
w = (E2|e2)
w = w().expand()
print '(E2|e2)/E^2 =',w/Esq
w = (E3|e3)
w = w().expand()
print '(E3|e3)/E^2 =',w/Esq
print '\nExtracting vectors from conformal 2 blade B = P1^P2'
metric = ' 0 -1 #,'+ \
'-1 0 #,'+ \
' # # #,'
MV.setup('P1 P2 a',metric)
B = P1^P2
Bsq = B*B
print 'B^2 =',Bsq
ap = a-(a^B)*B
print "a' = a-(a^B)*B =",ap
Ap = ap+ap*B
Am = ap-ap*B
print "A+ = a'+a'*B =",Ap
print "A- = a'-a'*B =",Am
print '(A+)^2 =',Ap*Ap
print '(A-)^2 =',Am*Am
aB = a|B
print 'a.B =',aB
tf = time.time()
print 1000.0*(tf-ti)
