************************************** Geometric Algebra Module for Sympy ************************************** :Author: Alan Bromborsky .. |release| replace:: 0.10 .. % Complete documentation on the extended LaTeX markup used for Python .. % documentation is available in ``Documenting Python'', which is part .. % of the standard documentation for Python. It may be found online .. % at: .. % .. % http://www.python.org/doc/current/doc/doc.html .. % \lstset{language=Python} .. % \input{macros} .. % This is a template for short or medium-size Python-related documents, .. % mostly notably the series of HOWTOs, but it can be used for any .. % document you like. .. % The title should be descriptive enough for people to be able to find .. % the relevant document. .. % Increment the release number whenever significant changes are made. .. % The author and/or editor can define 'significant' however they like. .. % At minimum, give your name and an email address. You can include a .. % snail-mail address if you like. .. % This makes the Abstract go on a separate page in the HTML version; .. % if a copyright notice is used, it should go immediately after this. .. % .. % \ifhtml .. % \chapter*{Front Matter\label{front}} .. % \fi .. % Copyright statement should go here, if needed. .. % ... .. % The abstract should be a paragraph or two long, and describe the .. % scope of the document. .. topic:: Abstract This document describes the implementation of a geometric algebra module in python that utilizes the :mod:`sympy` symbolic algebra library. The python module :mod:`GA` has been developed for coordinate free calculations using the operations (geometric, outer, and inner products etc.) of geometric algebra. The operations can be defined using a completely arbitrary metric defined by the inner products of a set of arbitrary vectors or the metric can be restricted to enforce orthogonality and signature constraints on the set of vectors. In addition the module includes the geometric, outer (curl) and inner (div) derivatives and the ability to define a curvilinear coordinate system. The module requires the numpy and the sympy modules. What is Geometric Algebra? ========================== Geometric algebra is the Clifford algebra of a real finite dimensional vector space or the algebra that results when a real finite dimensional vector space is extended with a product of vectors (geometric product) that is associative, left and right distributive, and yields a real number for the square (geometric product) of any vector [Hestenes,Lasenby]. The elements of the geometric algebra are called multivectors and consist of the linear combination of scalars, vectros, and the geometric product of two or more vectors. The additional axioms for the geometric algebra are that for any vectors :math:`a`, :math:`b`, and :math:`c` in the base vector space: .. math:: :nowrap: \begin{equation*} \begin{array}{c} a\lp bc \rp = \lp ab \rp c \\ a\lp b+c \rp = ab+ac \\ \lp a + b \rp c = ac+bc \\ aa = a^{2} \in \Re \end{array} \end{equation*} By induction these also apply to any multivectors. Several software packages for numerical geometric algebra calculations are available from Doran-Lazenby group and the Dorst group. Symbolic packages for Clifford algebra using orthongonal bases such as :math:`e_{i}e_{j}+e_{j}e_{i} = 2\eta_{ij}`, where :math:`\eta_{ij}` is a numeric array are available in Maple and Mathematica. The symbolic algebra module, :mod:`GA`, developed for python does not depend on an orthogonal basis representation, but rather is generated from a set of :math:`n` arbitrary symbolic vectors, :math:`a_{1},a_{2},\dots,a_{n}` and a symbolic metric tensor :math:`g_{ij} = a_{i}\cdot a_{j}`. In order not to reinvent the wheel all scalar symbolic algebra is handled by the python module :mod:`sympy`. The basic geometic algebra operations will be implemented in python by defining a multivector class, MV, and overloading the python operators in Table :ref:`1 ` where ``A`` and ``B`` are any two multivectors (In the case of ``+``, ``-``, ``*``, ``^``, ``|``, ``<<``, and ``>>`` the operation is also defined if ``A`` or ``B`` is a sympy symbol or a sympy real number). .. _table1: .. csv-table:: :header: " Operation ", " Result " :widths: 10, 40 " ``A+B`` ", " sum of multivectors " " ``A-B`` ", " difference of multivectors " " ``A*B`` ", " geometric product " " ``A^B`` ", " outer product of multivectors " " ``A|B`` ", " inner product of multivectors " " ``AB`` or ``A>>B`` ", " right contraction of multivectors " Table :ref:`1 `. Multivector operations for symbolicGA The option to use ``<`` or ``<<`` for left contraction and ``>`` or ``>>`` is given since the ``<`` and ``>`` operators do not have r-forms (there are no *__rlt__()* and *__rlt__()* functions to overload) while ``<<`` and ``>>`` do have r-forms so that *x << A* and *x >> A* are allowed where *x* is a scalar (symbol or integer) and *A* is a multivector. With ``<`` and ``>`` we can only have mixed modes (scalars and multivectors) if the multivector is the first operand. .. note:: Except for ``<`` and ``>`` all the multivector operators have r-forms so that as long as one of the operands, left or right, is a multivector the other can be a multivector or a scalar (sympy symbol or integer). .. warning:: Note that the operator order precedence is determined by python and is not neccessarily that used by geometric algebra. It is **absolutely essential** to use parenthesis in multivector expressions containing ``^``, ``|``, ``<``, ``>``, ``<<`` and/or ``>>``. As an example let ``A`` and ``B`` be any two multivectors. Then ``A + A*B = A +(A*B)``, but ``A+A^B = (2*A)^B`` since in python the ``^`` operator has a lower precedence than the '+' operator. In geometric algebra the outer and inner products and the left and right contractions have a higher precedence than the geometric product and the geometric product has a higher precedence than addition and subtraction. In python the ``^``, ``|``, ``<``, ``>``, ``<<`` and ``>>`` all have a lower precedence than ``+`` and ``-`` while ``*`` has a higher precedence than ``+`` and ``-``. .. _vbm: Vector Basis and Metric ======================= The two structures that define the :class:`MV` (multivector) class are the symbolic basis vectors and the symbolic metric. The symbolic basis vectors are input as a string with the symbol name separated by spaces. For example if we are calculating the geometric algebra of a system with three vectors that we wish to denote as ``a0``, ``a1``, and ``a2`` we would define the string variable: ``basis = 'a0 a1 a2'`` that would be input into the multivector setup function. The next step would be to define the symbolic metric for the geometric algebra of the basis we have defined. The default metric is the most general and is the matrix of the following symbols .. _eq1: .. math:: :label: 1 :nowrap: \begin{equation*} g = \lbrk \begin{array}{ccc} a0**2 & (a0.a1) & (a0.a2) \\ (a0.a1) & a1**2 & (a1.a2) \\ (a0.a2) & (a1.a2) & a2**2 \\ \end{array} \rbrk \end{equation*} where each of the :math:`g_{ij}` is a symbol representing all of the dot products of the basis vectors. Note that the symbols are named so that :math:`g_{ij} = g_{ji}` since for the symbol function :math:`(a0.a1) \ne (a1.a0)`. Note that the strings shown in equation :ref:`1 ` are only used when the values of :math:`g_{ij}` are output (printed). In the :mod:`GA` module (library) the :math:`g_{ij}` symbols are stored in a static member list of the multivector class :class:`MV` as the double list *MV.metric* (:math:`g_{ij}` = *MV.metric[i][j]*). The default definition of :math:`g` can be overwritten by specifying a string that will define :math:`g`. As an example consider a symbolic representation for conformal geometry. Define for a basis ``basis = 'a0 a1 a2 n nbar'`` and for a metric ``metric = '# # # 0 0, # # # 0 0, # # # 0 0, 0 0 0 0 2, 0 0 0 2 0'`` then calling `MV.setup(basis,metric)` would initialize .. math:: :nowrap: \begin{equation*} g = \lbrk \begin{array}{ccccc} a0**2 & (a0.a1) & (a0.a2) & 0 & 0\\ (a0.a1) & a1**2 & (a1.a2) & 0 & 0\\ (a0.a2) & (a1.a2) & a2**2 & 0 & 0 \\ 0 & 0 & 0 & 0 & 2 \\ 0 & 0 & 0 & 2 & 0 \end{array} \rbrk \end{equation*} Here we have specified that :math:`n` and :math:`nbar` are orthonal to all the :math:`a`'s, :math:`n**2 = nbar**2 = 0`, and :math:`(n.nbar) = 2`. Using :math:`\#` in the metric definition string just tells the program to use the default symbol for that value. When ``MV.setup`` is called multivector representations of the basis local to the program are instantiated. For our first example that means that the symbolic vectors named ``a0``, ``a1``, and ``a2`` are created and made available to the programmer for future calculations. In addition to the basis vectors the :math:`g_{ij}` are also made available to the programer with the following convention. If ``a0`` and ``a1`` are basis vectors, then their dot products are denoted by ``a0sq``, ``a2sq``, and ``a0dota1`` for use as python program varibles. If you print ``a0sq`` the output would be ``a0**2`` and the output for ``a0dota1`` would be ``(a0.a1)`` as shown in equation :ref:`1 `. If the default value are overridden the new values are output by print. For examle if :math:`g_{00} = 0` then "``print a0sq``" would output "0." More generally, if ``metric`` is not a string, but a list of lists or a two dimension numpy array, it is assumed that each element of ``metric`` is symbolic variable so that the :math:`g_{ij}` could be defined as symbolic functions as well as variables. For example instead of letting :math:`g_{01} = (a0.a1)` we could have :math:`g_{01} = cos(theta)` where we use a symbolic :math:`\cos` function. .. note:: Additionally ``MV.setup`` has an option for an othogonal basis where the signature of the metric space is defined by a string. For example if the signature of the vector space is :math:`(1,1,1)` (Euclidian 3-space) set ``metric = '[1,1,1]'`` Likewise if the signature is that of spacetime, :math:`(1,-1,-1,-1)` then define ``metric = '[1,-1,-1,-1]'``. Representation and Reduction of Multivector Bases ================================================= In our symbolic geometric algebra we assume that all multivectors of interest to us can be obtained from the symbolic basis vectors we have input, via the different operations available to geometric algebra. The first problem we have is representing the general multivector in terms terms of the basis vectors. To do this we form the ordered geometric products of the basis vectors and develop an internal representation of these products in terms of python classes. The ordered geometric products are all multivectors of the form :math:`a_{i_{1}}a_{i_{2}}\dots a_{i_{r}}` where :math:`i_{1}` .. note:: The empty list, ``[]``, represents the scalar 1. .. _table2: :: MV.basislabel = ['1', ['a0', 'a1', 'a2'], ['a0a1', 'a0a2', 'a1a2'], ['a0a1a2']] MV.basis = [[], [[0], [1], [2]], [[0, 1], [0, 2], [1, 2]], [[0, 1, 2]]] Table :ref:`2 `. Multivector basis labels and internal basis representation. Since there are :math:`2^{n}` bases and the number of bases with equal list lengths is the same as for the grade decomposition of a dimension :math:`n` geometric algebra we will call the collections of bases of equal length ``psuedogrades``. The critical operation in setting up the geometric algebra module is reducing the geomertric product of any two bases to a linear combination of bases so that we can calculate a multiplication table for the bases. First we represent the product as the concatenation of two base lists. For example ``a1a2*a0a1`` is represented by the list ``[1,2]+[0,1] = [1,2,0,1]``. The representation of the product is reduced via two operations, contraction and revision. The state of the reduction is saved in two lists of equal length. The first list contains symbolic scale factors (symbol or numeric types) for the corresponding interger list representing the product of bases. If we wish to reduce :math:`\lbrk i_{1},\dots,i_{r}\rbrk` the starting point is the coefficient list :math:`C = \lbrk 1 \rbrk` and the bases list :math:`B = \lbrk \lbrk i_{1},\dots,i_{r}\rbrk \rbrk`. We now operate on each element of the lists as follows: * ``contraction``: Consider a basis list :math:`B` with element :math:`B[j] = \lbrk i_{1},\dots,i_{l},i_{l+1},\dots,i_{s}\rbrk` where :math:`i_{l} = i_{l+1}`. Then the product of the :math:`l` and :math:`l+1` terms result in a scalar and :math:`B[j]` is replaced by the new list representation :math:`\lbrk i_{1},\dots,i_{l-1},i_{l+2},\dots,i_{r}\rbrk` which is of psuedo grade :math:`r-2` and :math:`C[j]` is replaced by the symbol :math:`g_{i_{l}i_{l}}C[j]`. * ``revision``: Consider a basis list :math:`B` with element :math:`B[j] = \lbrk i_{1},\dots,i_{l},i_{l+1},\dots,i_{s}\rbrk` where :math:`i_{l} > i_{l+1}`. Then the :math:`l` and :math:`l+1` elements must be reversed to be put in normal order, but we have :math:`a_{i_{l}}a_{i_{l+1}} = 2g_{i_{l}i_{l+1}}-a_{i_{l+1}}a_{i_{l}}` (From the geometric algebra definition of the dot product of two vectors). Thus we append the list representing the reduced element, :math:`\lbrk i_{1},\dots,i_{l-1},i_{l+2},\dots,i_{s}\rbrk`, to the pseudo bases list, :math:`B`, and append :math:`2g_{i{l}i_{l+1}}C[j]` to the coefficients list, then we replace :math:`B[j]` with :math:`\lbrk i_{1},\dots,i_{l+1},i_{l},\dots,i_{s}\rbrk` and :math:`C[j]` with :math:`-C[j]`. Both lists are increased by one element if :math:`g_{i_{l}i_{l+1}} \ne 0`. These processes are repeated untill every basis list in :math:`B` is in normal (ascending) order with no repeated elements. Then the coefficents of equivalent bases are summed and the bases sorted according to psuedograde and ascending order. We now have a way of calculating the geometric product of any two bases as a symbolic linear combination of all the bases with the coefficients determined by :math:`g`. The base multiplication table for our simple example of three vectors is given by (the coefficient of each psuedo base is enclosed with {} for clarity): .. include:: multable.dat :literal: Base Representation of Multivectors =================================== In terms of the bases defined an arbitrary multivector can be represented as a list of arrays (we use the numpy python module to implement arrays). If we have :math:`n` basis vectors we initialize the list ``self.mv = [0,0,...,0]`` with :math:`n+1` integers all zero. Each zero is a placeholder for an array of python objects (in this case the objects will be sympy symbol objects). If ``self.mv[r] = numpy.array([list of symbol objects])`` each entry in the ``numpy.array`` will be a coefficient of the corresponding psuedo base. ``self.mv[r] = 0`` indicates that the coefficients of every base of psuedo grade :math:`r` are 0. The length of the array ``self.mv[r]`` is :math:`n \choose r` the binomial coefficient. For example the psuedo basis vector ``a1`` would be represented as a multivector by the list: ``a1.mv = [0,numpy.array([numeric(0),numeric(1),numeric(0)]),0,0]`` and ``a0a1a2`` by: ``a0a1a2.mv = [0,0,0,numpy.array([numeric(1)])]`` The array is stuffed with sympy numeric objects instead of python integers so that we can perform symbolically manipulate sympy expressions that consist of scalar algebraic symbols and exact rational numbers which sympy can also represent. The ``numpy.array`` is used because operations of addition, substraction, and multiplication by an object are defined for the array if they are defined for the objects making up the array, which they are by sympy. We call this representation a base type because the ``r`` index is not a grade index since the bases we are using are not blades. In a blade representation the structure would be identical, but the bases would be replaced by blades and ``self.mv[r]`` would represent the ``r`` grade components of the multivector. The first use of the base representation is to store the results of the multiplication tabel for the bases in the class variable ``MV.mtabel``. This variable is a group of nested lists so that the geometric product of the ``igrade`` and ``ibase`` with the ``jgrade`` and ``jbase`` is ``MV.mtabel[igrade][ibase][jgrade][jbase]``. We can then use this table to calculate the geometric product of any two multivectors. Blade Representation of Multivectors ==================================== Since we can now calculate the symbolic geometric product of any two multivectors we can also calculate the blades corresponding to the product of the symbolic basis vectors using the formula .. math:: :nowrap: \begin{equation*} A_{r}\W b = \half\lp A_{r}b-\lp -1 \rp^{r}bA_{r} \rp, \end{equation*} where :math:`A_{r}` is a multivector of grade :math:`r` and :math:`b` is a vector. For our example basis the result is shown in Table :ref:`3 `. .. _table3: :: 1 = 1 a0 = a0 a1 = a1 a2 = a2 a0^a1 = {-(a0.a1)}1+a0a1 a0^a2 = {-(a0.a2)}1+a0a2 a1^a2 = {-(a1.a2)}1+a1a2 a0^a1^a2 = {-(a1.a2)}a0+{(a0.a2)}a1+{-(a0.a1)}a2+a0a1a2 Table :ref:`3 `. Bases blades in terms of bases. The important thing to notice about Table :ref:`3 ` is that it is a triagonal (lower triangular) system of equations so that using a simple back substitution algorithym we can solve for the psuedo bases in terms of the blades giving Table :ref:`4 `. .. _table4: :: 1 = 1 a0 = a0 a1 = a1 a2 = a2 a0a1 = {(a0.a1)}1+a0^a1 a0a2 = {(a0.a2)}1+a0^a2 a1a2 = {(a1.a2)}1+a1^a2 a0a1a2 = {(a1.a2)}a0+{-(a0.a2)}a1+{(a0.a1)}a2+a0^a1^a2 Table :ref:`4 `. Bases in terms of basis blades. Using Table :ref:`4 ` and simple substitution we can convert from a base multivector representation to a blade representation. Likewise, using Table :ref:`3 ` we can convert from blades to bases. Using the blade representation it becomes simple to program functions that will calculate the grade projection, reverse, even, and odd multivector functions. Note that in the multivector class ``MV`` there is a class variable for each instantiation, ``self.bladeflg``, that is set to zero for a base representation and 1 for a blade representation. One needs to keep track of which representation is in use since various multivector operations require conversion from one representation to the other. .. warning:: When the geometric product of two multivectors is calculated the module looks to see if either multivector is in blade representation. If either is the result of the geometric product is converted to a blade representation. One result of this is that if either of the multivectors is a simple vector (which is automatically a blade) the result will be in a blade representation. If ``a`` and ``b`` are vectors then the result ``a*b`` will be ``(a.b)+a^b`` or simply ``a^b`` if ``(a.b) = 0``. Outer and Inner Products, Left and Right Contractions ===================================================== In geometric algebra any general multivector :math:`A` can be decomposed into pure grade multivectors (a linear combination of blades of all the same order) so that in a :math:`n`-dimensional vector space .. math:: :nowrap: \begin{equation*} A = \sum_{r = 0}^{n}A_{r} \end{equation*} The geometric product of two pure grade multivectors :math:`A_{r}` and :math:`B_{s}` has the form .. math:: :nowrap: \begin{equation*} A_{r}B_{s} = \proj{A_{r}B_{s}}{\abs{r-s}}+\proj{A_{r}B_{s}}{\abs{r-s}+2}+\cdots+\proj{A_{r}B_{s}}{r+s} \end{equation*} where :math:`\proj{}{t}` projects the :math:`t` grade components of the multivector argument. The inner and outer products of :math:`A_{r}` and :math:`B_{s}` are then defined to be .. math:: :nowrap: \begin{equation*} A_{r}\cdot B_{s} = \proj{A_{r}B_{s}}{\abs{r-s}} \end{equation*} .. math:: :nowrap: \begin{equation*} A_{r}\wedge B_{s} = \proj{A_{r}B_{s}}{r+s} \end{equation*} and .. math:: :nowrap: \begin{equation*} A\cdot B = \sum_{r,s}A_{r}\cdot B_{s} \end{equation*} .. math:: :nowrap: \begin{equation*} A\wedge B = \sum_{r,s}A_{r}\wedge B_{s} \end{equation*} Likewise the right (:math:`\lfloor`) and left (:math:`\rfloor`) contractions are defined as .. math:: :nowrap: \begin{equation*} A_{r}\lfloor B_{s} = \left \{ \begin{array}{cc} \proj{A_{r}B_{s}}{r-s} & r \ge s \\ 0 & r < s \end{array} \right \} \end{equation*} .. math:: :nowrap: \begin{equation*} A_{r}\rfloor B_{s} = \left \{ \begin{array}{cc} \proj{A_{r}B_{s}}{s-r} & s \ge r \\ 0 & s < r \end{array} \right \} \end{equation*} and .. math:: :nowrap: \begin{equation*} A\lfloor B = \sum_{r,s}A_{r}\lfloor B_{s} \end{equation*} .. math:: :nowrap: \begin{equation*} A\rfloor B = \sum_{r,s}A_{r}\rfloor B_{s} \end{equation*} The ``MV`` class function for the outer product of the multivectors ``mv1`` and ``mv2`` is :: @staticmethod def outer_product(mv1,mv2): product = MV() product.bladeflg = 1 mv1.convert_to_blades() mv2.convert_to_blades() for igrade1 in MV.n1rg: if not isint(mv1.mv[igrade1]): pg1 = mv1.project(igrade1) for igrade2 in MV.n1rg: igrade = igrade1+igrade2 if igrade <= MV.n: if not isint(mv2.mv[igrade2]): pg2 = mv2.project(igrade2) pg1pg2 = pg1*pg2 product.add_in_place(pg1pg2.project(igrade)) return(product) The steps for calculating the outer product are: #. Convert ``mv1`` and ``mv2`` to blade representation if they are not already in that form. #. Project and loop through each grade ``mv1.mv[i1]`` and ``mv2.mv[i2]``. #. Calculate the geometric product ``pg1*pg2``. #. Project the ``i1+i2`` grade from ``pg1*pg2``. #. Accumulate the results for each pair of grades in the input multivectors. .. warning:: In the ``MV`` class we have overloaded the ``^`` operator to represent the outer product so that instead of calling the outer product function we can write ``mv1^ mv2``. Due to the precedence rules for python it is **absolutely essential** to enclose outer products in parenthesis. For the inner product of the multivectors ``mv1`` and ``mv2`` the ``MV`` class function is :: @staticmethod def inner_product(mv1,mv2,mode='s'): """ MV.inner_product(mv1,mv2) calculates the inner mode = 's' - symmetic (Doran & Lasenby) mode = 'l' - left contraction (Dorst) mode = 'r' - right contraction (Dorst) """ if isinstance(mv1,MV) and isinstance(mv2,MV): product = MV() product.bladeflg = 1 mv1.convert_to_blades() mv2.convert_to_blades() for igrade1 in range(MV.n1): if isinstance(mv1.mv[igrade1],numpy.ndarray): pg1 = mv1.project(igrade1) for igrade2 in range(MV.n1): igrade = igrade1-igrade2 if mode == 's': igrade = igrade.__abs__() else: if mode == 'l': igrade = -igrade if igrade >= 0: if isinstance(mv2.mv[igrade2],numpy.ndarray): pg2 = mv2.project(igrade2) pg1pg2 = pg1*pg2 product.add_in_place(pg1pg2.project(igrade)) return(product) else: if mode == 's': if isinstance(mv1,MV): product = mv1.scalar_mul(mv2) if isinstance(mv2,MV): product = mv2.scalar_mul(mv1) else: product = None return(product) The inner product is calculated the same way as the outer product except that in step 4, ``i1+i2`` is replaced by ``abs(i1-i2)`` or ``i1-i2`` for the right contraction or ``i2-i1`` for the left contraction. If ``i1-i2`` is less than zero there is no contribution to the right contraction. If ``i2-i1`` is less than zero there is no contribution to the left contraction. .. warning:: In the ``MV`` class we have overloaded the ``|`` operator for the inner product, ``>`` operator for the right contraction, and ``<`` operator for the left contraction. Instead of calling the inner product function we can write ``mv1|mv2``, ``mv1>mv2``, or ``mv1` holds for the geometric product of orthogonal vectors. The reverse is important in the theory of rotations in :math:`n`-dimensions. If :math:`R` is the product of an even number of vectors and :math:`RR^{\R} = 1` then :math:`RaR^{\R}` is a composition of rotations of the vector :math:`a`. If :math:`R` is the product of two vectors then the plane that :math:`R` defines is the plane of the rotation. That is to say that :math:`RaR^{\R}` rotates the component of :math:`a` that is projected into the plane defined by :math:`a` and :math:`b` where :math:`R=ab`. :math:`R` may be written :math:`R = e^{\frac{\theta}{2}U}`, where :math:`\theta` is the angle of rotation and :math:`u` is a unit blade :math:`\lp u^{2} = \pm 1\rp` that defines the plane of rotation. .. _recframe: Reciprocal Frames ================= If we have :math:`M` linearly independent vectors (a frame), :math:`a_{1},\dots,a_{M}`, then the reciprocal frame is :math:`a^{1},\dots,a^{M}` where :math:`a_{i}\cdot a^{j} = \delta_{i}^{j}`, :math:`\delta_{i}^{j}` is the Kronecker delta (zero if :math:`i \ne j` and one if :math:`i = j`). The reciprocal frame is constructed as follows: .. math:: :nowrap: \begin{equation*} E_{M} = a_{1}\W\dots\W a_{M} \end{equation*} .. math:: :nowrap: \begin{equation*} E_{M}^{-1} = \bfrac{E_{M}}{E_{M}^{2}} \end{equation*} Then .. math:: :nowrap: \begin{equation*} a^{i} = \lp -1\rp^{i-1}\lp a_{1}\W\dots\W \breve{a}_{i} \W\dots\W a_{M}\rp E_{M}^{-1} \end{equation*} where :math:`\breve{a}_{i}` indicates that :math:`a_{i}` is to be deleted from the product. In the standard notation if a vector is denoted with a subscript the reciprocal vector is denoted with a superscript. The multivector setup function ``MV.setup(basis,metric,rframe)`` has the argument ``rframe`` with a default value of ``False``. If it is set to ``True`` the reciprocal frame of the basis vectors is calculated. Additionaly there is the function ``reciprocal_frame(vlst,names='')`` external to the ``MV`` class that will calculate the reciprocal frame of a list, ``vlst``, of vectors. If the argument ``names`` is set to a space delimited string of names for the vectors the reciprocal vectors will be given these names. .. _deriv: Geometric Derivative ==================== If :math:`F` is a multivector field that is a function of a vector :math:`x = x^{i}\bm{\gamma}_{i}` (we are using the summation convention that pairs of subscripts and superscripts are summed over the dimension of the vector space) then the geometric derivative :math:`\nabla F` is given by (in this section the summation convention is used): .. math:: :nowrap: \begin{equation*} \nabla F = \bm{\gamma}^{i}\bfrac{\partial F}{\partial x^{i}} \end{equation*} If :math:`F_{R}` is a grade-:math:`R` multivector and :math:`F_{R} = F_{R}^{i_{1}\dots i_{R}}\bm{\gamma}_{i_{1}}\W\dots\W \bm{\gamma}_{i_{R}}` then .. math:: :nowrap: \begin{equation*} \nabla F_{R} = \bfrac{\partial F_{R}^{i_{1}\dots i_{R}}}{\partial x^{j}}\bm{\gamma}^{j}\lp\bm{\gamma}_{i_{1}}\W \dots\W \bm{\gamma}_{i_{R}} \rp \end{equation*} Note that :math:`\bm{\gamma}^{j}\lp\bm{\gamma}_{i_{1}}\W\dots\W \bm{\gamma}_{i_{R}} \rp` can only contain grades :math:`R-1` and :math:`R+1` so that :math:`\nabla F_{R}` also can only contain those grades. For a grade-:math:`R` multivector :math:`F_{R}` the inner (div) and outer (curl) derivatives are defined as .. math:: :nowrap: \begin{equation*} \nabla\cdot F_{R} = \left < \nabla F_{R}\right >_{R-1} \end{equation*} and .. math:: :nowrap: \begin{equation*} \nabla\W F_{R} = \left < \nabla F_{R}\right >_{R+1} \end{equation*} For a general multivector function :math:`F` the inner and outer derivatives are just the sum of the inner and outer dervatives of each grade of the multivector function. Curvilinear coordinates are derived from a vector function :math:`x(\bm{\theta})` where :math:`\bm{\theta} = \lp\theta_{1},\dots,\theta_{N}\rp` where the number of coordinates is equal to the dimension of the vector space. In the case of 3-dimensional spherical coordinates :math:`\bm{\theta} = \lp r,\theta,\phi \rp` and the coordinate generating function :math:`x(\bm{\theta})` is .. math:: :nowrap: \begin{equation*} x = r \cos\left({\phi}\right) \sin\left({\theta}\right){\bm{{\gamma}_{x}}}+ r \sin\left({\phi}\right) \sin\left({\theta}\right){\bm{{\gamma}_{y}}}+ r \cos\left({\theta}\right){\bm{{\gamma}_{z}}} \end{equation*} A coordinate frame is derived from :math:`x` by :math:`\bm{e}_{i} = \pdiff{x}{\theta^{i}}`. The following show the frame for spherical coordinates. .. math:: :nowrap: \begin{equation*} \bm{e}_{r} = \cos\left({\phi}\right) \sin\left({\theta}\right){\bm{{\gamma}_{x}}}+\sin\left({\phi}\right) \sin\left({\theta}\right){\bm{{\gamma}_{y}}}+\cos\left({\theta}\right){\bm{{\gamma}_{z}}} \end{equation*} .. math:: :nowrap: \begin{equation*} \bm{e}_{{\theta}} = \cos\left({\phi}\right) \cos\left({\theta}\right){\bm{{\gamma}_{x}}}+r \cos\left({\theta}\right) \sin\left({\phi}\right){\bm{{\gamma}_{y}}} - r \sin\left({\theta}\right){\bm{{\gamma}_{z}}} \end{equation*} .. math:: :nowrap: \begin{equation*} \bm{e}_{{\phi}} = - r \sin\left({\phi}\right) \sin\left({\theta}\right){\bm{{\gamma}_{x}}}+r \cos\left({\phi}\right) \sin\left({\theta}\right){\bm{{\gamma}_{y}}} \end{equation*} The coordinate frame generated in this manner is not necessarily normalized so define a normalized frame by .. math:: :nowrap: \begin{equation*} \bm{\hat{e}}_{i} = \bfrac{\bm{e}_{i}}{\sqrt{\abs{\bm{e}_{i}^{2}}}} = \bfrac{\bm{e}_{i}}{\abs{\bm{e}_{i}}} \end{equation*} This works for all :math:`\bm{e}_{i}^{2} \neq 0` since we have defined :math:`\abs{\bm{e}_{i}} = \sqrt{\abs{\bm{e}_{i}^{2}}}`. For spherical coordinates the normalized frame vectors are .. math:: :nowrap: \begin{equation*} \bm{\hat{e}}_{r} = \cos\left({\phi}\right) \sin\left({\theta}\right){\bm{{\gamma}_{x}}}+\sin\left({\phi}\right) \sin\left({\theta}\right){\bm{{\gamma}_{y}}}+\cos\left({\theta}\right){\bm{{\gamma}_{z}}} \end{equation*} .. math:: :nowrap: \begin{equation*} \bm{\hat{e}}_{{\theta}} = \cos\left({\phi}\right) \cos\left({\theta}\right){\bm{{\gamma}_{x}}}+\cos\left({\theta}\right) \sin\left({\phi}\right){\bm{{\gamma}_{y}}}- \sin\left({\theta}\right){\bm{{\gamma}_{z}}} \end{equation*} .. math:: :nowrap: \begin{equation*} \bm{\hat{e}}_{{\phi}} = - \sin\left({\phi}\right){\bm{{\gamma}_{x}}}+\cos\left({\phi}\right){\bm{{\gamma}_{y}}} \end{equation*} The geometric derivative in curvilinear coordinates is given by .. math:: :nowrap: \begin{align*} \nabla F_{R} & = \bm{\gamma}^{i}\pdiff{}{x^{i}}\lp F_{R}^{i_{1}\dots i_{R}} \bm{\hat{e}}_{i_{1}}\W\dots\W\bm{\hat{e}}_{i_{R}}\rp \\ & = \bm{e^{j}}\pdiff{}{\theta^{j}}\lp F_{R}^{i_{1}\dots i_{R}} \bm{\hat{e}}_{i_{1}}\W\dots\W\bm{\hat{e}}_{i_{R}}\rp \\ & = \lp\pdiff{}{\theta^{j}} F_{R}^{i_{1}\dots i_{R}}\rp \bm{e^{j}}\lp\bm{\hat{e}}_{i_{1}}\W\dots\W\bm{\hat{e}}_{i_{R}}\rp+ F_{R}^{i_{1}\dots i_{R}}\bm{e^{j}} \pdiff{}{\theta^{j}}\lp\bm{\hat{e}}_{i_{1}}\W\dots\W\bm{\hat{e}}_{i_{R}}\rp \\ & = \lp\pdiff{}{\theta^{j}} F_{R}^{i_{1}\dots i_{R}}\rp \bm{e^{j}}\lp\bm{\hat{e}}_{i_{1}}\W\dots\W\bm{\hat{e}}_{i_{R}}\rp+ F_{R}^{i_{1}\dots i_{R}}C\lbrc \bm{\hat{e}}_{i_{1}}\W\dots\W\bm{\hat{e}}_{i_{R}}\rbrc \end{align*} where .. math:: :nowrap: \begin{equation*} C\lbrc \bm{\hat{e}}_{i_{1}}\W\dots\W\bm{\hat{e}}_{i_{R}}\rbrc = \bm{e^{j}}\pdiff{}{\theta^{j}} \lp\bm{\hat{e}}_{i_{1}}\W\dots\W\bm{\hat{e}}_{i_{R}}\rp \end{equation*} are the connection multivectors for the curvilinear coordinate system. For a spherical coordinate system they are .. math:: :nowrap: \begin{equation*} C\lbrc\bm{\hat{e}}_{r}\rbrc = \frac{2}{r} \end{equation*} .. math:: :nowrap: \begin{equation*} C\lbrc\bm{\hat{e}}_{\theta}\rbrc = \frac{\cos\left({\theta}\right)}{r \sin\left({\theta}\right)} +\frac{1}{r}\bm{\hat{e}}_{r}\W\bm{\hat{e}}_{\theta} \end{equation*} .. math:: :nowrap: \begin{equation*} C\lbrc\bm{\hat{e}}_{\phi}\rbrc = \frac{1}{r}{\bm{\bm{\hat{e}}_{r}}}\W\bm{\hat{e}}_{{\phi}}+ \frac{\cos\left({\theta}\right)}{r \sin\left({\theta}\right)}\bm{\hat{e}}_{{\theta}}\W\bm{\hat{e}}_{{\phi}} \end{equation*} .. math:: :nowrap: \begin{equation*} C\lbrc\hat{e}_{r}\W\hat{e}_{\theta}\rbrc = - \frac{\cos\left({\theta}\right)}{r \sin\left({\theta}\right)} \bm{\hat{e}}_{r}+\frac{1}{r}\bm{\hat{e}}_{{\theta}} \end{equation*} .. math:: :nowrap: \begin{equation*} C\lbrc\bm{\hat{e}}_{r}\W\bm{\hat{e}}_{\phi}\rbrc = \frac{1}{r}\bm{\hat{e}}_{{\phi}} - \frac{\cos\left({\theta}\right)}{r \sin\left({\theta}\right)}\bm{\hat{e}}_{r}\W\bm{\hat{e}}_{{\theta}}\W\bm{\hat{e}}_{{\phi}} \end{equation*} .. math:: :nowrap: \begin{equation*} C\lbrc\bm{\hat{e}}_{\theta}\W\bm{\hat{e}}_{\phi}\rbrc = \frac{2}{r}\bm{\hat{e}}_{r}\W \bm{\hat{e}}_{\theta}\W\bm{\hat{e}}_{\phi} \end{equation*} .. math:: :nowrap: \begin{equation*} C\lbrc\bm{\hat{e}}_r\W\bm{\hat{e}}_{\theta}\W\bm{\hat{e}}_{\phi}\rbrc = 0 \end{equation*} Module Components ================= Initializing Multivector Class ------------------------------ The multivector class is initialized with: .. function:: MV.setup(basis,metric='',rframe=False,coords=None,debug=False,offset=0) The *basis* and *metric* parameters were described in section :ref:`vbm`. If *rframe=True* the reciprocal frame of the symbolic bases vectors is calculated. If *debug=True* the data structure required to initialize the :class:`MV` class are printer out. *coords* is a list of :class:`sympy` symbols equal in length to the number of basis vectors. These symbols are used as the arguments of a multivector field as a function of position and for calculating the derivatives of a multivector field (if *coords* is defined then *rframe* is automatically set equal to *True*). *offset* is an integer that is added to the multivector coefficient index. For example if one wishes to start labeling vector coefficient indexes at one instead of zero then set *offset=1*. Additionally, :func:`MV.setup` calculates the pseudo scalar, :math:`I` and its inverse, :math:`I^{-1}` and makes them available to the programmer as *MV.I* and *MV.Iinv*. After :func:`MV.setup` is run one can reinialize the :class:`MV` class with curvilinear coordinates using: .. function:: MV.rebase(x,coords,base_name,debug=False,debug_level=0) A typical usage of ``MV.rebase`` for generating spherical curvilinear coordinate is:: metric = '1 0 0,0 1 0,0 0 1' MV.setup('gamma_x gamma_y gamma_z',metric,True) coords = make_symbols('r theta phi') x = r*(sympy.cos(theta)*gamma_z+sympy.sin(theta)*\ (sympy.cos(phi)*gamma_x+sympy.sin(phi)*gamma_y)) x.set_name('x') MV.rebase(x,coords,'e',True) The input parameters for ``MV.rebase`` are * ``x``: Vector function of coordinates (derivatives define curvilinear basis) * ``coords``: List of sympy symbols for curvilinear coordinates * ``debug``: If ``True`` printout (LaTeX) all quantities required for derivative calculation * ``debug_level``: Set to 0,1,2, or 3 to stop curvilinear calculation before all quatities are calculated. This is done when debugging new curvilinear coordinate systems since simplification of expressions is not sufficiently automated to insure success of process of any coordinate system defined by vector function ``x`` To date ``MV.rebase`` works for cylindrical and spherical coordinate systems in any number of dimensions (until the execution time becomes too long). To make it work for these systems required creating some hacks for expression simplification since both trigsimp and simplify were not general enough to perform the required simplification. .. function:: MV.set_str_format(str_mode=0) If ``str_mode=0`` the string representation of the multivector contains no newline characters (prints on one line). If ``str_mode=1`` the string representation of the multivector places a newline after each grade of the multivector (prints one grade per line). If ``str_mode=2`` the string representation of the multivector places a newline after each base of the multivector (prints one base per line). In both cases bases with zero coefficients are not printed. .. note:: This function directly affects the way multivectors are printed with the print command since it interacts with the :func:`__str__` function for the multivector class which is used by the ``print`` command. .. csv-table:: :header: " ``str_mode`` ", " Effect on print " :widths: 10, 50 " 0 ","One multivector per line" " 1 ","One grade per line" " 2 ","One base per line" Instantiating a Multivector --------------------------- Now that grades and bases have been described we can show all the ways that a multivector can be instantiated. As an example assume that the multivector space is initialized with ``MV.setup('e1 e2 e3')``. Then the vectors ``e1``, ``e2``, and ``e3`` are made available (broadcast) for use in the program . .. warning:: This is only true if the statement ``set_main(sys.modules[__name__])`` appears immediately after the ``from sympy.galgebra.GA import *`` statement. So that multivectors could be instantiated with statements such as (``a1``, ``a2``, and ``a3`` are ``sympy`` symbols):: x = a1*e1+a2*e2+a3*e3 y = x*e1*e2 z = x|y w = x^y or with the multivector class constructor: .. class:: MV(value='',mvtype='',mvname='',fct=False) *mvname* is a string that defines the name of the multivector for output purposes. *value* and *type* are defined by the following table and *fct* is a switch that will convert the symbolic coefficients of a multivector to functions if coordinate variables have been defined when :func:`MV.setup` is called: .. csv-table:: :header: " mvtype ", " value ", " result " :widths: 10, 30, 30 " default ", " default ", " Zero multivector " " 'basisvector' ", " int i ", " :math:`\mbox{i}^{th}` basis vector " " 'basisbivector' ", " int i ", " :math:`\mbox{i}^{th}` basis bivector " " 'scalar' ", " symbol x ", " symbolic scalar of value x" " ", " string s ", " symbolic scalar of value Symbol(s) " " ", " int i ", " sympy integer of value i" " 'grade' ", " [int r, 1-D symbol array A] ", " X.grade(r) = A " " ", " [int r, string s] ", " symbolic grade r multivector " " 'vector' ", " 1-D symbol array A ", " X.grade(1) = A " " ", " string s ", " symbolic vector " " 'grade2' ", " 1-D symbol array A ", " X.grade(2) = A " " 'pseudo' ", " symbol x ", " X.grade(n) = x " " 'spinor' ", " string s ", " symbolic even multivector " " default "," string s ", " symbolic general multivector " If the *value* argument has the option of being a string s then a general symbolic multivector will be constructed constrained by the value of *mvtype*. The string s will be the base name of the multivector symbolic coefficients. If *coords* is not defined in :func:`MV.setup` the indices of the multivector bases are appended to the base name with a double underscore (superscript notation). If *coords* is defined the coordinate names will replace the indices in the coefficient names. For example if the base string is *A* and the coordinates *(x,y,z)* then the coefficients of a spinor in 3d space would be *A*, *A__xy*, *A__xz*, and *A__yz*. If the :mod:`latex_ex` is used to print the multivector the coefficients would print as :math:`A`, :math:`A^{xy}`, :math:`A^{xz}`, and :math:`A^{yz}`. If the *fct* argrument of :func:`MV` is set to *True* and the *coords* argument in :func:`MV.setup` is defined the symbolic coefficients of the multivector are functions of the coordinates. Basic Multivector Class Functions --------------------------------- .. function:: __call__(self,igrade=0,ibase=0) ``__call__`` returns the the ``igrade``, ``ibase`` coefficient of the multivector. The defaults return the scalar component of the multivector. .. function:: convert_to_blades(self) Convert multivector from the base representation to the blade representation. If multivector is already in blade representation nothing is done. .. function:: convert_from_blades(self) Convert multivector from the blade representation to the base representation. If multivector is already in base representation nothing is done. .. function:: project(self,r) If r is a integer return the grade-:math:`r` components of the multivector. If r is a multivector return the grades of the multivector that correspond to the non-zero grades of r. For example if one is projecting a general multivector and r is a spinor, ``A.project(r)`` will return only the even grades of the multivector A since a spinor only has even grades that are non-zero. .. function:: even(self) Return the even grade components of the multivector. .. function:: odd(self) Return the odd grade components of the multivector. .. function:: rev(self) Return the reverse of the multivector. See section :ref:`reverse`. .. function:: is_pure(self) Return true if multivector is pure grade (all grades but one are zero). .. function:: diff(self,x) Return the partial derivative of the multivector function with respect to variable :math:`x`. .. function:: grad(self) Return the geometric derivative of the multivector function. .. function:: grad_ext(self) Return the outer (curl) derivative of the multivector function. Equivalent to :func:`curl`. .. function:: grad_int(self) Return the inner (div) derivative of the multivector function. Equivalent to :func:`div`. .. warning:: If *A* is a vector field in three dimensions :math:`\nabla\cdot {\bf A}` = A.grad_int() = A.div(), but :math:`\nabla\times {\bf A}` = -MV.I*A.grad_ext() = -MV.I*A.curl(). Note that grad_int() lowers the grade of all blades by one grade and grad_ext() raises the grade of all blades by one. .. function:: set_coef(self,grade,base,value) Set the multivector coefficient of index *(grade,base)* to *value*. Sympy Functions Applied Inplace to Multivector Coefficients ----------------------------------------------------------- All the following fuctions belong to the :class:`MV` class and apply the corresponding :mod:`sympy` function to each component of a multivector. All these functions perform the operations inplace (``None`` is returned) on each coefficient. For example if you wished to simplify all the components of the multivector ``A`` you would invoke ``A.simplify()``. The argument list for each function is the same as for the corresponding :mod:`sympy` function. The only function that differs in its action from the :mod:`sympy` version is :func:`trigsimp` in its case the function ``TrigSimp is applied`` (see documentation on :func:`TrigSimp`). .. function:: collect(self,lst) .. function:: sqrfree(self,lst) .. function:: subs(self,*args) .. function:: simplify(self) .. function:: trigsimp(self) .. function:: cancel(self) .. function:: trim(self) .. function:: expand(self) Helper Functions ---------------- These are functions in :mod:`GA`, but not in the multivector (:class:`MV`) class. .. function:: set_main(main_program) :func:`set_main` passes the argument *main_program* from the main program to the :mod:`GA` module. The argument must be *sys.modules[__name__]* and the call should be placed immediately after :mod:`sys` and :mod:`GA` are imported. The purpose of this call is to allow :mod:`GA` to broadcast to the main program :mod:`sympy` variables and multivectors created by calls to :mod:`GA`. It is used by :func:`MV.setup` and :func:`make_symbols`. .. function:: make_symbols(symnamelst) :func:`make_symbols` creates a list of :mod:`sympy` symbols with names defined by the space delimited string *symnamelst*. In addition to returning the symbol list the function broadcasts the named symbols to the main program. For example if you make the call:: syms = make_symbols('x y ab') Not only will *syms* contain the symbols, but you can also directly use *x*, *y*, and *ab* as symbols in your program. .. warning:: You can only directly use *x*, *y*, and *ab* as symbols in your program if the statement ``set_main(sys.modules[__name__])`` appears immediately after the ``from sympy.galgebra.GA import *`` statement. .. function:: set_names(var_lst,var_str) :func:`set_names` allows one to name a list, *var_lst*, of multivectors enmass. The names are in *var_str*, a blank separated string of names. An error is generated if the number of name is not equal to the length of *var_lst*. .. function:: reciprocal_frame(vlst,names='') :func:`reciprocal_frame` implements the proceedure described in section :ref:`recframe`. *vlst* is a list of independent vectors that you wish the reciprocal frame calculated for. *names* is a blank separated string of names for the reciprocal vectors if names are required by you application. The function returns a list containing the reciprocal vectors. .. function:: TrigSimp(f) In general :func:`sympy.trigsimp` will not catch all the trigonometric simplifications in an :mod:`sympy` expression. Neither will :func:`TrigSimp`, but it will catch a lot more of them. :func:`TrigSimp` is so simple it is show below in its entirety. All it does is apply :func:`sympy.trigsimp` to the expressions generated by :func:`sympy.cse`. :: def TrigSimp(f): (w,g) = sympy.cse(f) g = sympy.trigsimp(g[0]) for sub in reversed(w): g = g.subs(sub[0],sub[1]) g = sympy.trigsimp(g) return(g) .. function:: S(x) :func:`S` instanciates a scaler multivector of value *x*, where *x* can be a :mod:`sympy` variable or integer. This is just a shorthand method for constructing scalar multivectors and can be used when there is any ambiguity in a multivector expression as to whether a symbol or constant should be treated as a scalar multivector or not. Examples ======== Algebra ------- The following examples of geometric algebra (not calculus) are all in the file :program:`testsymbolicGA.py` which is included in the sympy distribution examples under the galbebra directory. The section of code in the program for each example is show with the respective output following the code section. Example Header ^^^^^^^^^^^^^^ This is the header of :program:`testsymbolicGA.py` that allows access to the required modules and also allow variables and certain multivectors to be broadcast from the :mod:`GA` module to the main program. .. include:: headerGAtest.py :literal: Basic Geometric Algebra Operations ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ Example of basic geometric algebra operation of geometric, outer, and inner products. .. include:: BasicGAtest.py :literal: Examples of Conformal Geometry ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ Examples of comformal geometry [Lasenby,Chapter 10]. The examples show that basic geometric entities (lines, circles, planes, and spheres) in three dimensions can be represented by blades in a five dimensional (conformal) space. .. include:: conformalgeometryGAtest.py :literal: Calculation of Reciprocal Frame ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ Example shows the calculation of the reciprocal frame for three arbitrary vectors and verifies that the calculated reciprocal vectors have the correct properties. .. include:: reciprocalframeGAtest.py :literal: Hyperbolic Geometry ^^^^^^^^^^^^^^^^^^^ Examples of calculation of distance in hyperbolic geometry [Lasenby,pp373-375]. This is a good example of the utility of not restricting the basis vector to be orthogonal. Note that most of the calculation is simplifying a scalar expression. .. include:: hyperbolicGAtest.py :literal: Calculus -------- The calculus examples all use the extened LaTeXoutput module, ``latex_ex``, for clarity. Maxwell's Equations ^^^^^^^^^^^^^^^^^^^ In geometric calculus the equivalent of the electromagnetic tensor is :math:`F = E\gamma_{0}+IB\gamma_{0}` where a spacetime vector is given by :math:`x = x^{0}\gamma_{0}+x^{1}\gamma_{1}+x^{2}\gamma_{2}+x^{3}\gamma_{3}` where :math:`x^{0} = ct`, the pseudoscalar :math:`I = \gamma_{0}\gamma_{1}\gamma_{2}\gamma_{3}`, :math:`E` and :math:`B` are four vectors where the time component is zero and the spatial components equal to the electric and magnetic field components. Then Maxwell's equations can be all written as :math:`\nabla F = J` with :math:`J` the four current. This example shows that this equations generates all of Maxwell's equations correctly (in our units:math:`c=1`) [Lasenby,pp229-231]. ``Begin Program Maxwell.py`` .. include:: Maxwell.py :literal: ``End Program Maxwell.py`` ``Begin Program Output`` :math:`I` Pseudo-Scalar .. math:: :nowrap: \begin{equation*} I = {\gamma}_{t}{\gamma}_{x}{\gamma}_{y}{\gamma}_{z} \end{equation*} :math:`B` Magnetic Field Bi-Vector .. math:: :nowrap: \begin{equation*} B = - {B^{x}}{\gamma}_{t}{\gamma}_{x}- {B^{y}}{\gamma}_{t}{\gamma}_{y}- {B^{z}}{\gamma}_{t}{\gamma}_{z} \end{equation*} :math:`F` Electric Field Bi-Vector .. math:: :nowrap: \begin{equation*} E = - {E^{x}}{\gamma}_{t}{\gamma}_{x}- {E^{y}}{\gamma}_{t}{\gamma}_{y}- {E^{z}}{\gamma}_{t}{\gamma}_{z} \end{equation*} :math:`E+IB` Electo-Magnetic Field Bi-Vector .. math:: :nowrap: \begin{equation*} F = - {E^{x}}{\gamma}_{t}{\gamma}_{x}- {E^{y}}{\gamma}_{t}{\gamma}_{y}- {B^{z}}{\gamma}_{x}{\gamma}_{y}- {E^{z}}{\gamma}_{t}{\gamma}_{z}+ {B^{y}}{\gamma}_{x}{\gamma}_{z}- {B^{x}}{\gamma}_{y}{\gamma}_{z} \end{equation*} :math:`J` Four Current .. math:: :nowrap: \begin{equation*} J = {J^{t}}{\gamma}_{t}+ {J^{x}}{\gamma}_{x}+ {J^{y}}{\gamma}_{y}+ {J^{z}}{\gamma}_{z} \end{equation*} Geometric Derivative of Electo-Magnetic Field Bi-Vector .. math:: :nowrap: \begin{align*} \nabla F & = \left(\partial_{z} {E^{z}} + \partial_{y} {E^{y}} + \partial_{x} {E^{x}}\right){\gamma}_{t} \\ & + \left(-\partial_{t} {E^{x}} + \partial_{y} {B^{z}} -\partial_{z} {B^{y}}\right){\gamma}_{x} \\ & + \left(\partial_{z} {B^{x}} -\partial_{t} {E^{y}} -\partial_{x} {B^{z}}\right){\gamma}_{y} \\ & + \left(-\partial_{y} {B^{x}} -\partial_{t} {E^{z}} + \partial_{x} {B^{y}}\right){\gamma}_{z} \\ & + \left(-\partial_{x} {E^{y}} -\partial_{t} {B^{z}} + \partial_{y} {E^{x}}\right){\gamma}_{t}{\gamma}_{x}{\gamma}_{y} \\ & + \left(-\partial_{x} {E^{z}} + \partial_{t} {B^{y}} + \partial_{z} {E^{x}}\right){\gamma}_{t}{\gamma}_{x}{\gamma}_{z} \\ & + \left(-\partial_{t} {B^{x}} -\partial_{y} {E^{z}} + \partial_{z} {E^{y}}\right){\gamma}_{t}{\gamma}_{y}{\gamma}_{z} \\ & + \left(\partial_{y} {B^{y}} + \partial_{z} {B^{z}} + \partial_{x} {B^{x}}\right){\gamma}_{x}{\gamma}_{y}{\gamma}_{z}\end{align*} All Maxwell Equations are .. math:: :nowrap: \begin{equation*} \nabla F = J \end{equation*} Div :math:`E` and Curl :math:`H` Equations .. math:: :nowrap: \begin{align*} <\nabla F>_1 -J & = \left(-{J^{t}} + \partial_{z} {E^{z}} + \partial_{y} {E^{y}} + \partial_{x} {E^{x}}\right){\gamma}_{t} \\ & + \left(-{J^{x}} -\partial_{t} {E^{x}} + \partial_{y} {B^{z}} -\partial_{z} {B^{y}}\right){\gamma}_{x} \\ & + \left(\partial_{z} {B^{x}} -\partial_{t} {E^{y}} -{J^{y}} -\partial_{x} {B^{z}}\right){\gamma}_{y} \\ & + \left(-\partial_{y} {B^{x}} -\partial_{t} {E^{z}} -{J^{z}} + \partial_{x} {B^{y}}\right){\gamma}_{z}\end{align*} .. math:: :nowrap: \begin{equation*} = 0 \end{equation*} Curl :math:`E` and Div :math:`B` equations .. math:: :nowrap: \begin{align*} <\nabla F>_3 & = \left(-\partial_{x} {E^{y}} -\partial_{t} {B^{z}} + \partial_{y} {E^{x}}\right){\gamma}_{t}\W {\gamma}_{x}\W {\gamma}_{y} \\ & + \left(-\partial_{x} {E^{z}} + \partial_{t} {B^{y}} + \partial_{z} {E^{x}}\right){\gamma}_{t}\W {\gamma}_{x}\W {\gamma}_{z} \\ & + \left(-\partial_{t} {B^{x}} -\partial_{y} {E^{z}} + \partial_{z} {E^{y}}\right){\gamma}_{t}\W {\gamma}_{y}\W {\gamma}_{z} \\ & + \left(\partial_{y} {B^{y}} + \partial_{z} {B^{z}} + \partial_{x} {B^{x}}\right){\gamma}_{x}\W {\gamma}_{y}\W {\gamma}_{z}\end{align*} .. math:: :nowrap: \begin{equation*} = 0 \end{equation*} ``End Program Output`` Dirac's Equation ^^^^^^^^^^^^^^^^ The geometric algebra/calculus allows one to formulate the Dirac equation in real terms (no :math:`\sqrt{-1}`). Spinors :math:`\lp\Psi\rp` are even multivectors in space time (Minkowski space with signature (1,-1,-1,-1)) and the Dirac equation becomes :math:`\nabla \Psi I \sigma_{z}-eA\Psi = m\Psi\gamma_{t}`. All the terms in the real equation are explined in Doran and Lasenby [Lasenby,pp281-283]. ``Begin Program Dirac.py`` .. include:: Dirac.py :literal: ``End Program Dirac.py`` ``Begin Program Output`` :math:`A` is 4-vector potential .. math:: :nowrap: \begin{equation*} A = {A^{t}}{\gamma}_{t}+ {A^{x}}{\gamma}_{x}+ {A^{y}}{\gamma}_{y}+ {A^{z}}{\gamma}_{z} \end{equation*} :math:`\bm{\psi}` is 8-component real spinor (even multi-vector) .. math:: :nowrap: \begin{equation*} {}\bm{{\psi}} = {{\psi}}+ {{\psi}^{tx}}{\gamma}_{t}{\gamma}_{x}+ {{\psi}^{ty}}{\gamma}_{t}{\gamma}_{y}+ {{\psi}^{xy}}{\gamma}_{x}{\gamma}_{y}+ {{\psi}^{tz}}{\gamma}_{t}{\gamma}_{z}+ {{\psi}^{xz}}{\gamma}_{x}{\gamma}_{z}+ {{\psi}^{yz}}{\gamma}_{y}{\gamma}_{z}+ {{\psi}^{txyz}}{\gamma}_{t}{\gamma}_{x}{\gamma}_{y}{\gamma}_{z} \end{equation*} Dirac equation in terms of real geometric algebra/calculus :math:`\lp\nabla \bm{\psi} I \sigma_{z}-eA\bm{\psi} = m\bm{\psi}\gamma_{t}\rp` Spin measured with respect to :math:`z` axis .. math:: :nowrap: \begin{align*} \nabla \bm{\psi} I \sigma_{z}-eA\bm{\psi}-m\bm{\psi}\gamma_{t} & = \left(-e {A^{y}} {{\psi}^{ty}} -m {{\psi}} -\partial_{z} {{\psi}^{txyz}} -\partial_{y} {{\psi}^{tx}} -e {A^{z}} {{\psi}^{tz}} -e {A^{x}} {{\psi}^{tx}} + \partial_{x} {{\psi}^{ty}} -e {A^{t}} {{\psi}} + \partial_{t} {{\psi}^{xy}}\right){\gamma}_{t} \\ & + \left(-e {A^{y}} {{\psi}^{xy}} + \partial_{z} {{\psi}^{yz}} -e {A^{z}} {{\psi}^{xz}} -e {A^{t}} {{\psi}^{tx}} + m {{\psi}^{tx}} -\partial_{x} {{\psi}^{xy}} + \partial_{y} {{\psi}} -e {A^{x}} {{\psi}} -\partial_{t} {{\psi}^{ty}}\right){\gamma}_{x} \\ & + \left(-e {A^{y}} {{\psi}} + e {A^{x}} {{\psi}^{xy}} -e {A^{t}} {{\psi}^{ty}} -\partial_{x} {{\psi}} + \partial_{t} {{\psi}^{tx}} -\partial_{z} {{\psi}^{xz}} -e {A^{z}} {{\psi}^{yz}} + m {{\psi}^{ty}} -\partial_{y} {{\psi}^{xy}}\right){\gamma}_{y} \\ & + \left(-\partial_{z} {{\psi}^{xy}} + e {A^{x}} {{\psi}^{xz}} -e {A^{t}} {{\psi}^{tz}} -\partial_{x} {{\psi}^{yz}} + \partial_{y} {{\psi}^{xz}} + e {A^{y}} {{\psi}^{yz}} + \partial_{t} {{\psi}^{txyz}} + m {{\psi}^{tz}} -e {A^{z}} {{\psi}}\right){\gamma}_{z} \\ & + \left(\partial_{y} {{\psi}^{ty}} -e {A^{y}} {{\psi}^{tx}} + e {A^{x}} {{\psi}^{ty}} + \partial_{z} {{\psi}^{tz}} -e {A^{z}} {{\psi}^{txyz}} + \partial_{x} {{\psi}^{tx}} -e {A^{t}} {{\psi}^{xy}} -\partial_{t} {{\psi}} -m {{\psi}^{xy}}\right){\gamma}_{t}{\gamma}_{x}{\gamma}_{y} \\ & + \left(-\partial_{y} {{\psi}^{tz}} + \partial_{x} {{\psi}^{txyz}} -e {A^{t}} {{\psi}^{xz}} + e {A^{x}} {{\psi}^{tz}} + e {A^{y}} {{\psi}^{txyz}} + \partial_{z} {{\psi}^{ty}} -m {{\psi}^{xz}} -e {A^{z}} {{\psi}^{tx}} -\partial_{t} {{\psi}^{yz}}\right){\gamma}_{t}{\gamma}_{x}{\gamma}_{z} \\ & + \left(-e {A^{t}} {{\psi}^{yz}} -\partial_{z} {{\psi}^{tx}} + \partial_{x} {{\psi}^{tz}} -m {{\psi}^{yz}} + \partial_{y} {{\psi}^{txyz}} + \partial_{t} {{\psi}^{xz}} -e {A^{z}} {{\psi}^{ty}} -e {A^{x}} {{\psi}^{txyz}} + e {A^{y}} {{\psi}^{tz}}\right){\gamma}_{t}{\gamma}_{y}{\gamma}_{z} \\ & + \left(\partial_{z} {{\psi}} + e {A^{y}} {{\psi}^{xz}} + m {{\psi}^{txyz}} -\partial_{t} {{\psi}^{tz}} -\partial_{x} {{\psi}^{xz}} -e {A^{x}} {{\psi}^{yz}} -e {A^{t}} {{\psi}^{txyz}} -\partial_{y} {{\psi}^{yz}} -e {A^{z}} {{\psi}^{xy}}\right){\gamma}_{x}{\gamma}_{y}{\gamma}_{z}\end{align*} .. math:: :nowrap: \begin{equation*} = 0 \end{equation*} ``End Program Output`` Spherical Coordinates ^^^^^^^^^^^^^^^^^^^^^ Curvilinear coodinates are implemented as shown in section :ref:`deriv`. The gradient of a scalar function and the divergence and curl of a vector function (:math:`-I\lp\nabla\W A\rp` is the curl in three dimensions in the notation of geometric algebra) to demonstrate the formulas derived in section :ref:`deriv`. ``Begin Program coords.py`` .. include:: coords.py :literal: ``End Program coords.py`` ``Begin Program Output`` Gradient of Scalar Function :math:`\psi` .. math:: :nowrap: \begin{equation*} \nabla\psi = \partial_{r} {{\psi}}{}\bm{e}_{r}+\frac{1 \partial_{{\theta}} {{\psi}}}{r}{}\bm{e}_{{\theta}}+\frac{1 \partial_{{\phi}} {{\psi}}}{r \operatorname{sin}\left({\theta}\right)}{}\bm{e}_{{\phi}} \end{equation*} Div and Curl of Vector Function :math:`A` .. math:: :nowrap: \begin{equation*} A = {A^{r}}{}\bm{e}_{r}+ {A^{{\theta}}}{}\bm{e}_{{\theta}}+ {A^{{\phi}}}{}\bm{e}_{{\phi}} \end{equation*} .. math:: :nowrap: \begin{equation*} \nabla \cdot A = \left(\frac{{A^{{\theta}}} \operatorname{cos}\left({\theta}\right)}{r \operatorname{sin}\left({\theta}\right)} + 2 \frac{{A^{r}}}{r} + \partial_{r} {A^{r}} + \frac{\partial_{{\theta}} {A^{{\theta}}}}{r} + \frac{\partial_{{\phi}} {A^{{\phi}}}}{r \operatorname{sin}\left({\theta}\right)}\right) \end{equation*} .. math:: :nowrap: \begin{align*} -I\lp\nabla \W A\rp & = \left(\frac{\partial_{{\theta}} {A^{{\phi}}}}{r} + \frac{{A^{{\phi}}} \operatorname{cos}\left({\theta}\right)}{r \operatorname{sin}\left({\theta}\right)} -\frac{\partial_{{\phi}} {A^{{\theta}}}}{r \operatorname{sin}\left({\theta}\right)}\right){}\bm{e}_{r} \\ & - \left(\frac{{A^{{\phi}}}}{r} + \partial_{r} {A^{{\phi}}} -\frac{\partial_{{\phi}} {A^{r}}}{r \operatorname{sin}\left({\theta}\right)}\right){}\bm{e}_{{\theta}} \\ & + \left(\frac{{A^{{\theta}}}}{r} + \partial_{r} {A^{{\theta}}} -\frac{\partial_{{\theta}} {A^{r}}}{r}\right){}\bm{e}_{{\phi}}\end{align*} ``End Program Output`` .. seealso:: `Hestenes `_ ``Clifford Algebra to Geometric Calculus`` by D.Hestenes and G. Sobczyk, Kluwer Academic Publishers, 1984. `Lasenby `_ ``Geometric Algebra for Physicists`` by C. Doran and A. Lasenby, Cambridge University Press, 2003.