def test_differential_operators():

    xyz_coords = (x, y, z) = symbols('x y z', real=True)
    (o3d, ex, ey, ez) = Ga.build('e', g=[1, 1, 1], coords=xyz_coords)
    f = o3d.mv('f', 'scalar', f=True)
    lap = o3d.grad*o3d.grad

    assert str(lap) == 'D{x}^2 + D{y}^2 + D{z}^2'
    assert str(lap * f) == 'D{x}^2f + D{y}^2f + D{z}^2f'

    sph_coords = (r, th, phi) = symbols('r theta phi', real=True)
    (sp3d, er, eth, ephi) = Ga.build('e', g=[1, r**2, r**2 * sin(th)**2], coords=sph_coords, norm=True)
    f = sp3d.mv('f', 'scalar', f=True)
    lap = sp3d.grad*sp3d.grad
    assert str(lap) == '2/r*D{r} + cos(theta)/(r**2*sin(theta))*D{theta} + D{r}^2 + r**(-2)*D{theta}^2 + 1/(r**2*sin(theta)**2)*D{phi}^2'
    assert str(lap * f) == 'D{r}^2f + 2*D{r}f/r + D{theta}^2f/r**2 + cos(theta)*D{theta}f/(r**2*sin(theta)) + D{phi}^2f/(r**2*sin(theta)**2)'

    A = o3d.mv('A','vector')
    xs = o3d.mv(x)

    assert o3d.grad*A == 0
    assert str(A*o3d.grad) == 'A__x*D{x} + A__y*D{y} + A__z*D{z} + e_x^e_y*(-A__y*D{x} + A__x*D{y}) + e_x^e_z*(-A__z*D{x} + A__x*D{z}) + e_y^e_z*(-A__z*D{y} + A__y*D{z})'
    assert o3d.grad*xs == ex
    assert str(xs*o3d.grad) == 'e_x*x*D{x} + e_y*x*D{y} + e_z*x*D{z}'
    assert str(o3d.grad*(o3d.grad+xs)) == 'D{x}^2 + D{y}^2 + D{z}^2 + e_x*D{}'
    assert str((o3d.grad+xs)*o3d.grad) == 'D{x}^2 + D{y}^2 + D{z}^2 + e_x*x*D{x} + e_y*x*D{y} + e_z*x*D{z}'

    return

def basic_multivector_operations_3D():
    Print_Function()

    g3d = Ga('e*x|y|z')
    (ex,ey,ez) = g3d.mv()

    A = g3d.mv('A','mv')

    A.Fmt(1,'A')
    A.Fmt(2,'A')
    A.Fmt(3,'A')

    A.even().Fmt(1,'%A_{+}')
    A.odd().Fmt(1,'%A_{-}')

    X = g3d.mv('X','vector')
    Y = g3d.mv('Y','vector')

    print 'g_{ij} = ',g3d.g

    X.Fmt(1,'X')
    Y.Fmt(1,'Y')

    (X*Y).Fmt(2,'X*Y')
    (X^Y).Fmt(2,'X^Y')
    (X|Y).Fmt(2,'X|Y')
    return

def properties_of_geometric_objects():
    Print_Function()
    global n, nbar

    g = '# # # 0 0,'+ \
        '# # # 0 0,'+ \
        '# # # 0 0,'+ \
        '0 0 0 0 2,'+ \
        '0 0 0 2 0'

    c3d = Ga('p1 p2 p3 n nbar',g=g)

    (p1,p2,p3,n,nbar) = c3d.mv()

    print 'g_{ij} =\n',c3d.g

    P1 = F(p1)
    P2 = F(p2)
    P3 = F(p3)

    print 'Extracting direction of line from L = P1^P2^n'

    L = P1^P2^n
    delta = (L|n)|nbar
    print '(L|n)|nbar =',delta

    print 'Extracting 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)

def extracting_vectors_from_conformal_2_blade():
    Print_Function()

    g = '0 -1 #,'+ \
        '-1 0 #,'+ \
        '# # #'

    e2b = Ga('P1 P2 a',g=g)

    (P1,P2,a) = e2b.mv()

    print 'g_{ij} =\n',e2b.g

    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
    return

def test_reciprocal_frame():
    """
    Test of formula for general reciprocal frame of three vectors.
    Let three independent vectors be e1, e2, and e3. The reciprocal
    vectors E1, E2, and E3 obey the relations:

    e_i.E_j = delta_ij*(e1^e2^e3)**2
    """
    g = '1 # #,'+ \
        '# 1 #,'+ \
        '# # 1'

    g3dn = Ga('e1 e2 e3',g=g)

    (e1,e2,e3) = g3dn.mv()

    E = e1^e2^e3
    Esq = (E*E).scalar()
    Esq_inv = 1 / Esq

    E1 = (e2^e3)*E
    E2 = (-1)*(e1^e3)*E
    E3 = (e1^e2)*E

    w = (E1|e2)
    w = w.expand()
    assert w.scalar() == 0

    w = (E1|e3)
    w = w.expand()
    assert w.scalar() == 0

    w = (E2|e1)
    w = w.expand()
    assert w.scalar() == 0

    w = (E2|e3)
    w = w.expand()
    assert w.scalar() == 0

    w = (E3|e1)
    w = w.expand()
    assert w.scalar() == 0

    w = (E3|e2)
    w = w.expand()
    assert w.scalar() == 0

    w = (E1|e1)
    w = (w.expand()).scalar()
    Esq = expand(Esq)
    assert simplify(w/Esq) == 1

    w = (E2|e2)
    w = (w.expand()).scalar()
    assert simplify(w/Esq) == 1

    w = (E3|e3)
    w = (w.expand()).scalar()
    assert simplify(w/Esq) == 1

def derivatives_in_rectangular_coordinates():
    Print_Function()

    X = (x, y, z) = symbols('x y z')
    o3d = Ga('e_x e_y e_z', g=[1, 1, 1], coords=X)
    (ex, ey, ez) = o3d.mv()
    grad = o3d.grad

    f = o3d.mv('f', 'scalar', f=True)
    A = o3d.mv('A', 'vector', f=True)
    B = o3d.mv('B', 'bivector', f=True)
    C = o3d.mv('C', 'mv', f=True)
    print 'f =', f
    print 'A =', A
    print 'B =', B
    print 'C =', C

    print 'grad*f =', grad * f
    print 'grad|A =', grad | A
    print 'grad*A =', grad * A

    print '-I*(grad^A) =', -o3d.I() * (grad ^ A)
    print 'grad*B =', grad * B
    print 'grad^B =', grad ^ B
    print 'grad|B =', grad | B

    print 'grad<A =', grad < A
    print 'grad>A =', grad > A
    print 'grad<B =', grad < B
    print 'grad>B =', grad > B
    print 'grad<C =', grad < C
    print 'grad>C =', grad > C

    return

def rounding_numerical_components():
    Print_Function()
    o3d = Ga('e_x e_y e_z',g=[1,1,1])
    (ex,ey,ez) = o3d.mv()

    X = 1.2*ex+2.34*ey+0.555*ez
    Y = 0.333*ex+4*ey+5.3*ez

    print 'X =',X
    print 'Nga(X,2) =',Nga(X,2)
    print 'X*Y =',X*Y
    print 'Nga(X*Y,2) =',Nga(X*Y,2)
    return

def test_basic_multivector_operations():

    g3d, ex, ey, ez = Ga.build('e*x|y|z')

    A = g3d.mv('A', 'mv')

    assert str(A) == 'A + A__x*e_x + A__y*e_y + A__z*e_z + A__xy*e_x^e_y + A__xz*e_x^e_z + A__yz*e_y^e_z + A__xyz*e_x^e_y^e_z'

    X = g3d.mv('X', 'vector')
    Y = g3d.mv('Y', 'vector')

    assert str(X) == 'X__x*e_x + X__y*e_y + X__z*e_z'
    assert str(Y) == 'Y__x*e_x + Y__y*e_y + Y__z*e_z'

    assert str((X*Y)) == '(e_x.e_x)*X__x*Y__x + (e_x.e_y)*X__x*Y__y + (e_x.e_y)*X__y*Y__x + (e_x.e_z)*X__x*Y__z + (e_x.e_z)*X__z*Y__x + (e_y.e_y)*X__y*Y__y + (e_y.e_z)*X__y*Y__z + (e_y.e_z)*X__z*Y__y + (e_z.e_z)*X__z*Y__z + (X__x*Y__y - X__y*Y__x)*e_x^e_y + (X__x*Y__z - X__z*Y__x)*e_x^e_z + (X__y*Y__z - X__z*Y__y)*e_y^e_z'
    assert str((X ^ Y)) == '(X__x*Y__y - X__y*Y__x)*e_x^e_y + (X__x*Y__z - X__z*Y__x)*e_x^e_z + (X__y*Y__z - X__z*Y__y)*e_y^e_z'
    assert str((X | Y)) == '(e_x.e_x)*X__x*Y__x + (e_x.e_y)*X__x*Y__y + (e_x.e_y)*X__y*Y__x + (e_x.e_z)*X__x*Y__z + (e_x.e_z)*X__z*Y__x + (e_y.e_y)*X__y*Y__y + (e_y.e_z)*X__y*Y__z + (e_y.e_z)*X__z*Y__y + (e_z.e_z)*X__z*Y__z'

    g2d, ex, ey = Ga.build('e*x|y')

    X = g2d.mv('X', 'vector')
    A = g2d.mv('A', 'spinor')

    assert str(X) == 'X__x*e_x + X__y*e_y'
    assert str(A) == 'A + A__xy*e_x^e_y'

    assert str((X | A)) == '-A__xy*((e_x.e_y)*X__x + (e_y.e_y)*X__y)*e_x + A__xy*((e_x.e_x)*X__x + (e_x.e_y)*X__y)*e_y'
    assert str((X < A)) == '(-(e_x.e_y)*A__xy*X__x - (e_y.e_y)*A__xy*X__y + A*X__x)*e_x + ((e_x.e_x)*A__xy*X__x + (e_x.e_y)*A__xy*X__y + A*X__y)*e_y'
    assert str((A > X)) == '((e_x.e_y)*A__xy*X__x + (e_y.e_y)*A__xy*X__y + A*X__x)*e_x + (-(e_x.e_x)*A__xy*X__x - (e_x.e_y)*A__xy*X__y + A*X__y)*e_y'

    o2d, ex, ey = Ga.build('e*x|y', g=[1, 1])

    X = o2d.mv('X', 'vector')
    A = o2d.mv('A', 'spinor')

    assert str(X) == 'X__x*e_x + X__y*e_y'
    assert str(A) == 'A + A__xy*e_x^e_y'

    assert str((X*A)) == '(A*X__x - A__xy*X__y)*e_x + (A*X__y + A__xy*X__x)*e_y'
    assert str((X | A)) == '-A__xy*X__y*e_x + A__xy*X__x*e_y'
    assert str((X < A)) == '(A*X__x - A__xy*X__y)*e_x + (A*X__y + A__xy*X__x)*e_y'
    assert str((X > A)) == 'A*X__x*e_x + A*X__y*e_y'

    assert str((A*X)) == '(A*X__x + A__xy*X__y)*e_x + (A*X__y - A__xy*X__x)*e_y'
    assert str((A | X)) == 'A__xy*X__y*e_x - A__xy*X__x*e_y'
    assert str((A < X)) == 'A*X__x*e_x + A*X__y*e_y'
    assert str((A > X)) == '(A*X__x + A__xy*X__y)*e_x + (A*X__y - A__xy*X__x)*e_y'

    return

def test_extracting_vectors_from_conformal_2_blade():

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

    (cf1d, P1, P2, a) = Ga.build('P1 P2 a', g=metric)

    B = P1 ^ P2
    Bsq = B*B
    assert str(Bsq) == '1'
    ap = a - (a ^ B)*B
    assert str(ap) == '-(P2.a)*P1 - (P1.a)*P2'

    Ap = ap + ap*B
    Am = ap - ap*B

    assert str(Ap) == '-2*(P2.a)*P1'
    assert str(Am) == '-2*(P1.a)*P2'

    assert str(Ap*Ap) == '0'
    assert str(Am*Am) == '0'

    aB = a | B
    assert str(aB) == '-(P2.a)*P1 + (P1.a)*P2'

    return

def test_conformal_representations_of_circles_lines_spheres_and_planes():
    global n, nbar

    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'

    (cf3d, ex, ey, ez, n, nbar) = Ga.build('e_x e_y e_z n nbar', g=metric)

    x, y, z = symbols('x y z', real=True)

    e = n + nbar
    #conformal representation of points

    A = make_vector(ex, cf3d)
    B = make_vector(ey, cf3d)
    C = make_vector(-ex, cf3d)
    D = make_vector(ez, cf3d)
    X = make_vector(x*ex + y*ey +z*ez, cf3d)

    assert A == ex + (n - nbar)/S(2)
    assert B == ey + (n - nbar)/S(2)
    assert C == -ex + (n - nbar)/S(2)
    assert D == ez + (n - nbar)/S(2)
    assert X == x*ex + y*ey + z*ez + (x**2/2 + y**2/2 + z**2/2)*n - nbar/2

    assert A ^ B ^ C ^ X == -z*(ex^ey^ez^n) + z*(ex^ey^ez^nbar) + ((x**2 + y**2 + z**2 - S(1))/2)*(ex^ey^n^nbar)
    assert A ^ B ^ n ^ X == -z*(ex^ey^ez^n) + ((x + y - S(1))/2)*(ex^ey^n^nbar) + (z/2)*(ex^ez^n^nbar) - (z/2)*(ey^ez^n^nbar)
    assert A ^ B ^ C ^ D ^ X == ((-x**2 - y**2 - z**2 + S(1))/2)*ex^ey^ez^n^nbar
    assert A ^ B ^ n ^ D ^ X == ((-x - y - z + S(1))/2)*(ex^ey^ez^n^nbar)

    L = (A ^ B ^ e) ^ X

    assert L == -z*(ex^ey^ez^n) - z*(ex^ey^ez^nbar) + (-x**2/2 + x - y**2/2 + y - z**2/2 - S(1)/2)*(ex^ey^n^nbar) + z*(ex^ez^n^nbar) - z*(ey^ez^n^nbar)

    return

def test_derivatives_in_rectangular_coordinates():

    X = (x, y, z) = symbols('x y z')
    o3d, ex, ey, ez = Ga.build('e_x e_y e_z', g=[1,1,1], coords=X)
    grad = o3d.grad

    f = o3d.mv('f', 'scalar', f=True)
    A = o3d.mv('A', 'vector', f=True)
    B = o3d.mv('B', 'bivector', f=True)
    C = o3d.mv('C', 'mv', f=True)

    assert str(f) == 'f'
    assert str(A) == 'A__x*e_x + A__y*e_y + A__z*e_z'
    assert str(B) == 'B__xy*e_x^e_y + B__xz*e_x^e_z + B__yz*e_y^e_z'
    assert str(C) == 'C + C__x*e_x + C__y*e_y + C__z*e_z + C__xy*e_x^e_y + C__xz*e_x^e_z + C__yz*e_y^e_z + C__xyz*e_x^e_y^e_z'

    assert str(grad*f) == 'D{x}f*e_x + D{y}f*e_y + D{z}f*e_z'
    assert str(grad | A) == 'D{x}A__x + D{y}A__y + D{z}A__z'
    assert str(grad*A) == 'D{x}A__x + D{y}A__y + D{z}A__z + (-D{y}A__x + D{x}A__y)*e_x^e_y + (-D{z}A__x + D{x}A__z)*e_x^e_z + (-D{z}A__y + D{y}A__z)*e_y^e_z'

    assert str(-o3d.I()*(grad ^ A)) == '(-D{z}A__y + D{y}A__z)*e_x + (D{z}A__x - D{x}A__z)*e_y + (-D{y}A__x + D{x}A__y)*e_z'
    assert str(grad*B) == '(-D{y}B__xy - D{z}B__xz)*e_x + (D{x}B__xy - D{z}B__yz)*e_y + (D{x}B__xz + D{y}B__yz)*e_z + (D{z}B__xy - D{y}B__xz + D{x}B__yz)*e_x^e_y^e_z'
    assert str(grad ^ B) == '(D{z}B__xy - D{y}B__xz + D{x}B__yz)*e_x^e_y^e_z'
    assert str(grad | B) == '(-D{y}B__xy - D{z}B__xz)*e_x + (D{x}B__xy - D{z}B__yz)*e_y + (D{x}B__xz + D{y}B__yz)*e_z'

    assert str(grad < A) == 'D{x}A__x + D{y}A__y + D{z}A__z'
    assert str(grad > A) == 'D{x}A__x + D{y}A__y + D{z}A__z'
    assert str(grad < B) == '(-D{y}B__xy - D{z}B__xz)*e_x + (D{x}B__xy - D{z}B__yz)*e_y + (D{x}B__xz + D{y}B__yz)*e_z'
    assert str(grad > B) == '0'
    assert str(grad < C) == 'D{x}C__x + D{y}C__y + D{z}C__z + (D{x}C - D{y}C__xy - D{z}C__xz)*e_x + (D{y}C + D{x}C__xy - D{z}C__yz)*e_y + (D{z}C + D{x}C__xz + D{y}C__yz)*e_z + D{z}C__xyz*e_x^e_y - D{y}C__xyz*e_x^e_z + D{x}C__xyz*e_y^e_z'
    assert str(grad > C) == 'D{x}C__x + D{y}C__y + D{z}C__z + D{x}C*e_x + D{y}C*e_y + D{z}C*e_z'

    return

def Lorentz_Tranformation_in_Geog_Algebra():
    Print_Function()
    (alpha,beta,gamma) = symbols('alpha beta gamma')
    (x,t,xp,tp) = symbols("x t x' t'",real=True)
    (st2d,g0,g1) = Ga.build('gamma*t|x',g=[1,-1])

    from sympy import sinh,cosh

    R = cosh(alpha/2)+sinh(alpha/2)*(g0^g1)
    X = t*g0+x*g1
    Xp = tp*g0+xp*g1
    print 'R =',R

    print r"#%t\bm{\gamma_{t}}+x\bm{\gamma_{x}} = t'\bm{\gamma'_{t}}+x'\bm{\gamma'_{x}} = R\lp t'\bm{\gamma_{t}}+x'\bm{\gamma_{x}}\rp R^{\dagger}"

    Xpp = R*Xp*R.rev()
    Xpp = Xpp.collect()
    Xpp = Xpp.trigsimp()
    print r"%t\bm{\gamma_{t}}+x\bm{\gamma_{x}} =",Xpp
    Xpp = Xpp.subs({sinh(alpha):gamma*beta,cosh(alpha):gamma})

    print r'%\f{\sinh}{\alpha} = \gamma\beta'
    print r'%\f{\cosh}{\alpha} = \gamma'

    print r"%t\bm{\gamma_{t}}+x\bm{\gamma_{x}} =",Xpp.collect()
    return

def main():
    Print_Function()

    (a, b, c) = abc = symbols('a,b,c',real=True)
    (o3d, ea, eb, ec) = Ga.build('e_a e_b e_c', g=[1, 1, 1], coords=abc)
    grad = o3d.grad

    x = symbols('x',real=True)
    A = o3d.lt([[x*a*c**2,x**2*a*b*c,x**2*a**3*b**5],\
                [x**3*a**2*b*c,x**4*a*b**2*c**5,5*x**4*a*b**2*c],\
                [x**4*a*b**2*c**4,4*x**4*a*b**2*c**2,4*x**4*a**5*b**2*c]])
    print 'A =',A

    v = a*ea+b*eb+c*ec

    print 'v =',v

    f = v|A(v)

    print r'%f = v\cdot \f{A}{v} =',f

    (grad * f).Fmt(3,r'%\nabla f')

    Av = A(v)

    print r'%\f{A}{v} =', Av

    (grad * Av).Fmt(3,r'%\nabla \f{A}{v}')

    return

def main():
    Eprint()

    X = (x,y,z) = symbols('x y z',real=True)
    (o3d,ex,ey,ez) = Ga.build('e_x e_y e_z',g=[1,1,1],coords=(x,y,z))

    A = x*(ey^ez) + y*(ez^ex) + z*(ex^ey)
    print 'A =', A
    print 'grad^A =',(o3d.grad^A).simplify()
    print

    f = o3d.mv(1/sqrt(x**2 + y**2 + z**2))
    print 'f =', f
    print 'grad*f =',(o3d.grad*f).simplify()
    print

    B = f*A
    print 'B =', B
    print

    Curl_B = o3d.grad^B

    print 'grad^B =', Curl_B.simplify()

    return

def derivatives_in_spherical_coordinates():
    Print_Function()
    coords = (r,th,phi) = symbols('r theta phi', real=True)
    (sp3d,er,eth,ephi) = Ga.build('e_r e_theta e_phi',g=[1,r**2,r**2*sin(th)**2],coords=coords,norm=True)
    grad = sp3d.grad

    f = sp3d.mv('f','scalar',f=True)
    A = sp3d.mv('A','vector',f=True)
    B = sp3d.mv('B','bivector',f=True)

    print 'f =',f
    print 'A =',A
    print 'B =',B

    print 'grad*f =',grad*f
    print 'grad|A =',grad|A
    print 'grad\\times A = -I*(grad^A) =',-sp3d.i*(grad^A)
    print '%\\nabla^{2}f =',grad|(grad*f)
    print 'grad^B =',grad^B

    """
    print '( \\nabla\\W\\nabla )\\bm{e}_{r} =',((grad^grad)*er).trigsimp()
    print '( \\nabla\\W\\nabla )\\bm{e}_{\\theta} =',((grad^grad)*eth).trigsimp()
    print '( \\nabla\\W\\nabla )\\bm{e}_{\\phi} =',((grad^grad)*ephi).trigsimp()
    """

    return

def test_extract_plane_and_line():
    """
    Show that conformal trivector encodes planes and lines. See D&L section
    10.4.2
    """
    metric = '# # # 0 0,' + \
             '# # # 0 0,' + \
             '# # # 0 0,' + \
             '0 0 0 0 2,' + \
             '0 0 0 2 0'

    cf3d, p1, p2, p3, n, nbar = Ga.build('p1 p2 p3 n nbar', g=metric)

    P1 = F(p1, n, nbar)
    P2 = F(p2, n, nbar)
    P3 = F(p3, n, nbar)

    #Line through p1 and p2
    L = P1 ^ P2 ^ n
    delta = (L | n) | nbar
    delta_test = 2*p1 - 2*p2
    diff = delta - delta_test
    assert diff == S.Zero

    #Plane through p1, p2, and p3
    C = P1 ^ P2 ^ P3
    delta = ((C ^ n) | n) | nbar
    delta_test = 2*(p1 ^ p2) - 2*(p1 ^ p3) + 2*(p2 ^ p3)
    diff = delta - delta_test
    assert diff == S.Zero

def Maxwells_Equations_in_Geom_Calculus():
    Print_Function()
    X = symbols('t x y z',real=True)
    (st4d,g0,g1,g2,g3) = Ga.build('gamma*t|x|y|z',g=[1,-1,-1,-1],coords=X)

    I = st4d.i

    B = st4d.mv('B','vector',f=True)
    E = st4d.mv('E','vector',f=True)
    B.set_coef(1,0,0)
    E.set_coef(1,0,0)
    B *= g0
    E *= g0
    J = st4d.mv('J','vector',f=True)
    F = E+I*B

    print r'\text{Pseudo Scalar\;\;}I =',I
    print '\\text{Magnetic Field Bi-Vector\\;\\;} B = \\bm{B\\gamma_{t}} =',B
    print '\\text{Electric Field Bi-Vector\\;\\;} E = \\bm{E\\gamma_{t}} =',E
    print '\\text{Electromagnetic Field Bi-Vector\\;\\;} F = E+IB =',F
    print '%\\text{Four Current Density\\;\\;} J =',J
    gradF = st4d.grad*F
    print '#Geom Derivative of Electomagnetic Field Bi-Vector'
    gradF.Fmt(3,'grad*F')

    print '#Maxwell Equations'
    print 'grad*F = J'
    print '#Div $E$ and Curl $H$ Equations'
    (gradF.get_grade(1)-J).Fmt(3,'%\\grade{\\nabla F}_{1} -J = 0')
    print '#Curl $E$ and Div $B$ equations'
    (gradF.get_grade(3)).Fmt(3,'%\\grade{\\nabla F}_{3} = 0')
    return

def test_vector_extraction():
    """
    Show that conformal bivector encodes two points. See D&L Section 10.4.1
    """
    metric = '0 -1 #,' + \
             '-1 0 #,' + \
             '# # #'

    cext, P1, P2, a = Ga.build('P1 P2 a', g=metric)
    """
    P1 and P2 are null vectors and hence encode points in conformal space.
    Show that P1 and P2 can be extracted from the bivector B = P1^P2. a is a
    third vector in the conformal space with a.B not 0.
    """
    B = P1 ^ P2
    Bsq = B*B
    ap = a - (a ^ B)*B
    Ap = ap + ap*B
    Am = ap - ap*B

    P1dota = cext.g[0,2]
    P2dota = cext.g[1,2]
    Ap_test = (-2*P2dota)*P1
    Am_test = (-2*P1dota)*P2
    assert Ap == Ap_test
    assert Am == Am_test
    Ap2 = Ap*Ap
    Am2 = Am*Am
    assert Ap2 == S.Zero
    assert Am2 == S.Zero

def test_substitution():

    o3d, e_x, e_y, e_z = Ga.build('e_x e_y e_z', g='1 0 0, 0 1 0, 0 0 1')
    x, y, z = symbols('x y z')

    X = x*e_x + y*e_y + z*e_z
    Y = X.subs([(x, 2), (y, 3), (z, 4)])
    assert Y == 2*e_x + 3*e_y + 4*e_z

def basic_multivector_operations_2D():
    Print_Function()
    g2d = Ga('e*x|y')
    (ex,ey) = g2d.mv()

    print 'g_{ij} =',g2d.g

    X = g2d.mv('X','vector')
    A = g2d.mv('A','spinor')

    X.Fmt(1,'X')
    A.Fmt(1,'A')

    (X|A).Fmt(2,'X|A')
    (X<A).Fmt(2,'X<A')
    (A>X).Fmt(2,'A>X')
    return