#!/usr/bin/python
#test_symbolicGA.py

from GAsympy import *

set_main(sys.modules[__name__])

def F(x):
    """
    Conformal Mapping Function
    """
    Fx = HALF*((x*x)*n+2*x-nbar)
    #print 'F(x) =',Fx
    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__':

    """
    Note that is A and B are multivectors:

        A*B - Geometric Product
        A^B - Outer Product
        A|B - Inner Product

    Outer and Inner product terms must be enclosed in parenthesis
    since ^ and | have a low precedence than +, -, and *.

    """
    """
    Reference "Geometric Algebra for Physicists" pp 373-375
    """

    print 'Example: non-euclidian distance calculation'

    metric = '0 # #,'+ \
             '# 0 #,'+ \
             '# # 1,'

    MV.setup('X Y e',metric,debug=0)
    MV.set_str_format(1)

    """
    X and Y are conformal mappings of two vectors on the Poincare disk
    """

    L = X^Y^e

    """
    B is the bivector generator of translations on the circle defined by L
    """

    B = L*e

    print 'B =',B
    print 'B^2 =',(B*B)()
    print 'L^2 =',(L*L)()

    make_scalars('s c Binv')

    """
    s = sinh(alpha/2)
    c = cosh(alpha/2)
    Binv = 1/sqrt(B*B)  It can be shown that B*B > 0
    """

    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 =',R
    Z = R*X*R.rev()
    Z.expand()
    Z.collect([Binv,s,c,MV.metric[0][1]])
    print 'R*X*R.rev() =',Z
    W = Z^Y
    W.expand()
    W.collect([Binv])
    print '(R*X*rev(R))^Y =',W

    """
    Coefficients of W blades must be reduced via hyperbolic trig
    substitutions to yield solution for alpha.
    """
    print '(R*X*R.rev())^Y =',W
    print 'Setting blade coefficients to zero allows one to '+\
          'solve for alpha in terms of dot products!'

    print '\n\n\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)
    #MV.LaTeX()

    #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)

    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