# Source code for sympy.galgebra.GA

#!/usr/bin/python

"""
The module symbolicGA implements symbolic Geometric Algebra in python.
The relevant references for this module are:

1. "Geometric Algebra for Physicists" by C. Doran and A. Lazenby,
Cambridge University Press, 2003.

2. "Geometric Algebra for Computer Science" by Leo Dorst,
Daniel Fontijne, and Stephen Mann, Morgan Kaufmann Publishers, 2007.

3. SymPy Tutorial, http://docs.sympy.org/
"""
import sys
import numpy, sympy
import re as regrep
import sympy.galgebra.latex_ex
from sympy.core.decorators import deprecated

NUMPAT = regrep.compile( '([\-0-9])|([\-0-9]/[0-9])')
"""Re pattern for rational number"""

ZERO = sympy.Rational(0)
ONE  = sympy.Rational(1)
TWO  = sympy.Rational(2)
HALF = sympy.Rational(1,2)

from sympy.core import Pow as pow_type
from sympy import Abs as abs_type
from sympy.core import Mul as mul_type

@sympy.vectorize(0)
def substitute_array(array,*args):
return(array.subs(*args))

[docs]def is_quasi_unit_numpy_array(array):
"""
Determine if a array is square and diagonal with
entries of +1 or -1.
"""
shape = numpy.shape(array)
if len(shape) == 2 and (shape[0] == shape[1]):
n = shape[0]
ix = 0
while ix < n:
iy = 0
while iy < ix:
if array[ix][iy] != ZERO:
return(False)
iy += 1
if sympy.Abs(array[ix][ix]) != ONE:
return(False)
ix += 1
return(True)
else:
return(False)

@deprecated(value="This function is no longer needed.", issue=3379,
deprecated_since_version="0.7.2")
def set_main(main_program):
pass

def plist(lst):
if type(lst) == list:
for x in lst:
plist(x)
else:
sys.stderr.write(lst+'\n')
return

[docs]def numeric(num_str):
"""
Returns rational numbers compatible with symbols.
Input is a string representing a fraction or integer
or a simple integer.
"""
if type(num_str) == int:
a = num_str
b = 1
else:
tmp = num_str.split('/')
if len(tmp) == 1:
a = int(tmp[0])
b = 1
else:
a = int(tmp[0])
b = int(tmp[1])
return(sympy.Rational(a,b))

[docs]def collect(expr,lst):
"""
Wrapper for sympy.collect.
"""
lst = MV.scalar_to_symbol(lst)
return(sympy.collect(expr,lst))

def sqrt(expr):
return(sympy.sqrt(expr))

[docs]def isint(a):
"""
Test for integer.
"""
return(type(a) == int)

[docs]def make_null_array(n):
"""
Return list of n empty lists.
"""
a = []
for i in range(n):
a.append([])
return(a)

[docs]def test_int_flgs(lst):
"""
Test if all elements in list are 0.
"""
for i in lst:
if i:
return(1)
return(0)

[docs]def comb(N,P):
"""
Calculates the combinations of the integers [0,N-1] taken P at a time.
The output is a list of lists of integers where the inner lists are
the different combinations. Each combination is sorted in ascending
order.
"""
def rloop(n,p,combs,comb):
if p:
for i in range(n):
newcomb = comb+[i]
np = p-1
rloop(i,np,combs,newcomb)
else:
combs.append(comb)
combs = []
rloop(N,P,combs,[])
for comb in combs:
comb.sort()
return(combs)

[docs]def diagpq(p, q=0):
"""
Return string equivalent metric tensor for signature (p, q).
"""
n = p + q
D = []
rn = list(range(n))
for i in rn:
for j in rn:
if i ==j:
if i < p:
D.append('1 ')
else:
D.append('-1 ')
else:
D.append('0 ')
return ','.join(D)

[docs]def make_scalars(symnamelst):
"""
make_scalars takes a string of symbol names separated by
blanks and converts them to MV scalars and returns a list
of the symbols.
"""
symlst = sympy.symbols(symnamelst)
scalar_lst = []
for s in symlst:
tmp = MV(s,'scalar')
scalar_lst.append(tmp)
return(scalar_lst)

deprecated_since_version="0.7.2")
def make_symbols(symnamelst):
return sympy.symbols(symnamelst)

[docs]def israt(numstr):
"""
Test if string represents a rational number.
"""
global NUMPAT
if NUMPAT.match(numstr):
return(1)
return(0)

[docs]def dualsort(lst1, lst2):
"""
Inplace dual sort of lst1 and lst2 keyed on sorted lst1.
"""
_indices = list(range(len(lst1)))
_indices.sort(key=lst2.__getitem__)
lst1[:] = list(map(lst1.__getitem__, _indices))
lst2[:] = list(map(lst2.__getitem__, _indices))
return

[docs]def cp(A,B):
"""
Calculates the commutator product (A*B-B*A)/2 for
the objects A and B.
"""
return(HALF*(A*B-B*A))

[docs]def reduce_base(k,base):
"""
If base is a list of sorted integers [i_1,...,i_R] then reduce_base
sorts the list [k,i_1,...,i_R] and calculates whether an odd or even
number of permutations is required to sort the list. The sorted list
is returned and +1 for even permutations or -1 for odd permutations.
"""
if k in base:
return(0,base)
if k < base[0]:
return(1,[k,base[0]])
else:
return(-1,[base[0],k])
ilo = 0
if k < base[0]:
return(1,[k]+base)
if k > base[ihi]:
return(1,base+[k])
else:
return(-1,base+[k])
imid = ihi+ilo
return(-1,[base[0],k,base[1]])
while True:
if ihi-ilo == 1:
break
if base[imid] > k:
ihi = imid
else:
ilo = imid
imid = (ilo+ihi)/2
if ilo%2 == 1:
return(1,base[:ihi]+[k]+base[ihi:])
else:
return(-1,base[:ihi]+[k]+base[ihi:])

[docs]def sub_base(k,base):
"""
If base is a list of sorted integers [i_1,...,i_R] then sub_base returns
a list with the k^th element removed. Note that k=0 removes the first
element.  There is no test to see if k is in the range of the list.
"""
n = len(base)
if n == 1:
return([])
if n == 2:
if k == base[0]:
return([base[1]])
else:
return([base[0]])
return(base[:k]+base[k+1:])

[docs]def magnitude(vector):
"""
Calculate magnitude of vector containing trig expressions
and simplify.  This is a hack because of way he sign of
magsq is determined and because of the way that absolute
values are removed.
"""
magsq = sympy.expand((vector|vector)())
magsq = sympy.trigsimp(magsq,deep=True,recursive=True)
#print magsq
magsq_str = sympy.galgebra.latex_ex.LatexPrinter()._print(magsq)
if magsq_str[0] == '-':
magsq = -magsq
mag = unabs(sqrt(magsq))
#print mag
return(mag)

[docs]def LaTeX_lst(lst,title=''):
"""
Output a list in LaTeX format.
"""
if title != '':
sympy.galgebra.latex_ex.LaTeX(title)
for x in lst:
sympy.galgebra.latex_ex.LaTeX(x)
return

[docs]def unabs(x):
"""
Remove absolute values from expressions so a = sqrt(a**2).
This is a hack.
"""
if type(x) == mul_type:
y = unabs(x.args[0])
for yi in x.args[1:]:
y *= unabs(yi)
return(y)
if type(x) == pow_type:
if x.args[1] == HALF and type(x.args[0]) == add_type:
return(x)
y = 1/unabs(x.args[0])
return(y)
if len(x.args) == 0:
return(x)
if type(x) == abs_type:
return(x.args[0])
return(x)

def function_lst(fstr,xtuple):
sys.stderr.write(fstr+'\n')
fct_lst = []
for xstr in fstr.split():
f = sympy.Function(xstr)(*xtuple)
fct_lst.append(f)
return(fct_lst)

[docs]def vector_fct(Fstr,x):
"""
Create a list of functions of arguments x.  One function is
created for each variable in x.  Fstr is a string that is
the base name of each function while each function in the
list is given the name Fstr+'__'+str(x[ix]) so that if
Fstr = 'f' and str(x[1]) = 'theta' then the LaTeX output
of the second element in the output list would be 'f^{\\theta}'.
"""
nx = len(x)
Fvec = []
for ix in range(nx):
ftmp = sympy.Function(Fstr+'__'+sympy.galgebra.latex_ex.LatexPrinter.str_basic(x[ix]))(*tuple(x))
Fvec.append(ftmp)
return(Fvec)

def print_lst(lst):
for x in lst:
print(x)
return

[docs]def normalize(elst,nname_lst):
"""
Normalize a list of vectors and rename the normalized vectors. 'elist' is the list
(or array) of vectors to be normalized and nname_lst is a list of the names for the
normalized vectors.  The function returns the numpy arrays enlst and mags containing
the normalized vectors (enlst) and the magnitudes of the original vectors (mags).
"""
i = 0
mags = numpy.array(MV.n*[ZERO],dtype=numpy.object)
enlst= numpy.array(MV.n*[ZERO],dtype=numpy.object)
for (e,nname) in zip(elst,nname_lst):
emag = magnitude(e)
emaginv = 1/emag
mags[i] = emag
enorm = emaginv*e
enorm.name = nname
enlst[i] = enorm
i += 1
return(enlst,mags)

def build_base(base_index,base_vectors,reverse=False):
base = base_vectors[base_index[0]]
if len(base_index) > 1:
for i in base_index[1:]:
base = base^base_vectors[i]
if reverse:
base = base.rev()
return(base)

class MV(object):

is_setup = False
basislabel_lst = 0
curvilinear_flg = False
coords = None

@staticmethod
"""
Pad list with zeros to length n. If length is > n
"""
nvalue = len(value)
if nvalue < n:
value = value+(n-nvalue)*[ZERO]
if nvalue > n:
value = value[:n]
return(value)

@staticmethod
def define_basis(basis):
"""
Calculates all the MV static variables needed for
basis operations.  See reference 5 section 2.
"""
MV.vbasis     = basis
MV.vsyms      = sympy.symbols(MV.vbasis)
MV.n          = len(MV.vbasis)
MV.nrg        = list(range(MV.n))
MV.n1         = MV.n+1
MV.n1rg       = list(range(MV.n1))
MV.npow       = 2**MV.n
MV.index      = list(range(MV.n))
MV.gabasis    = [[]]
MV.basis      = (MV.n+1)*[0]
MV.basislabel = (MV.n+1)*[0]
MV.basis[0]   = []
MV.basislabel[0] = '1'
MV.basislabel_lst = [['1']]
MV.nbasis = numpy.array((MV.n+1)*[1],dtype=numpy.object)
MV.gabasis += [tmp]
ntmp = len(tmp)
for i in range(ntmp):
tmp_lst = []
bstr = ''
for j in tmp[i]:
bstr += MV.vbasis[j]
tmp_lst.append(MV.vbasis[j])
MV.basis_map = [{'':0}]
tmpdict = {}
ibases = 0
for base in bases:
tmpdict[str(base)] = ibases
ibases += 1
MV.basis_map.append(tmpdict)

if MV.debug:
print('basis strings =',MV.vbasis)
print('basis symbols =',MV.vsyms)
print('basis labels  =',MV.basislabel)
print('basis         =',MV.basis)
print('index         =',MV.index)
return

@staticmethod
def define_metric(metric):
"""
Calculates all the MV static variables needed for
metric operations.  See reference 5 section 2.
"""
if MV.metric_str:
MV.g = []
MV.metric = numpy.array(MV.n*[MV.n*[ZERO]],dtype=numpy.object)
if metric == '':
metric = numpy.array(MV.n*[MV.n*['#']],dtype=numpy.object)
for i in MV.index:
for j in MV.index:
gij = metric[i][j]
if israt(gij):
MV.metric[i][j] = numeric(gij)
else:
if gij == '#':
if i == j:
gij = '('+MV.vbasis[j]+'**2)'
else:
gij = '('+MV.vbasis[min(i,j)]+'.'+MV.vbasis[max(i,j)]+')'
tmp = sympy.Symbol(gij)
MV.metric[i][j] = tmp
if i <= j:
MV.g.append(tmp)
else:
MV.metric = metric
MV.g = []
for row in metric:
g_row = []
for col in metric:
g_row.append(col)
MV.g.append(g_row)
if MV.debug:
print('metric =',MV.metric)
return

@staticmethod
def define_reciprocal_frame():
"""
Calculates unscaled reciprocal vectors (MV.brecp) and scale
factor (MV.Esq). The ith scaled reciprocal vector is
(1/MV.Esq)*MV.brecp[i].  The pseudoscalar for the set of
basis vectors is MV.E.
"""
if MV.tables_flg:
MV.E = MV.bvec[0]
MV.brecp = []
for i in range(1,MV.n):
MV.E = MV.E^MV.bvec[i]
for i in range(MV.n):
tmp = ONE
if i%2 != 0:
tmp = -ONE
for j in range(MV.n):
if i != j:
tmp = tmp^MV.bvec[j]
tmp = tmp*MV.E
MV.brecp.append(tmp)
MV.Esq = (MV.E*MV.E)()
MV.Esq_inv = ONE/MV.Esq
for i in range(MV.n):
MV.brecp[i] = MV.brecp[i]*MV.Esq_inv
return

@staticmethod
def reduce_basis_loop(blst):
"""
Makes one pass through basis product representation for
reduction of representation to normal form.
See reference 5 section 3.
"""
nblst = len(blst)
if nblst <= 1:
return(1)
jstep = 1
while jstep <nblst:
istep = jstep-1
if blst[istep] == blst[jstep]:
i = blst[istep]
if len(blst) >2:
blst = blst[:istep]+blst[jstep+1:]
else:
blst = []
if len(blst) <= 1 or jstep == nblst-1:
blst_flg = 0
else:
blst_flg = 1
return(MV.metric[i][i],blst,blst_flg)
if blst[istep] > blst[jstep]:
blst1 = blst[:istep]+blst[jstep+1:]
a1 = TWO*MV.metric[blst[jstep]][blst[istep]]
blst = blst[:istep]+[blst[jstep]]+[blst[istep]]+blst[jstep+1:]
if len(blst1) <= 1:
blst1_flg = 0
else:
blst1_flg = 1
return(a1,blst1,blst1_flg,blst)
jstep +=1
return(1)

@staticmethod
def reduce_basis(blst):
"""
Repetitively applies reduce_basis_loop to basis
product representation until normal form is
realized.  See reference 5 section 3.
"""
if blst == []:
blst_coef   = []
blst_expand = []
for i in MV.n1rg:
blst_coef.append([])
blst_expand.append([])
blst_expand[0].append([])
blst_coef[0].append(ONE)
return(blst_coef,blst_expand)
blst_expand = [blst]
blst_coef      = [ONE]
blst_flg         = [1]
while test_int_flgs(blst_flg):
for i in range(len(blst_flg)):
if blst_flg[i]:
tmp = MV.reduce_basis_loop(blst_expand[i])
if tmp ==1:
blst_flg[i] = 0
else:
if len(tmp) == 3:
blst_coef[i] = tmp[0]*blst_coef[i]
blst_expand[i] = tmp[1]
blst_flg[i] = tmp[2]
else:
blst_coef[i]      = -blst_coef[i]
blst_flg[i]        = 1
blst_expand[i] = tmp[3]
blst_coef.append(-blst_coef[i]*tmp[0])
blst_expand.append(tmp[1])
blst_flg.append(tmp[2])
(blst_coef,blst_expand) = MV.combine_common_factors(blst_coef,blst_expand)
return(blst_coef,blst_expand)

@staticmethod
def combine_common_factors(blst_coef,blst_expand):
new_blst_coef = []
new_blst_expand = []
for i in range(MV.n1):
new_blst_coef.append([])
new_blst_expand.append([])
nfac = len(blst_coef)
for ifac in range(nfac):
blen = len(blst_expand[ifac])
new_blst_coef[blen].append(blst_coef[ifac])
new_blst_expand[blen].append(blst_expand[ifac])
for i in range(MV.n1):
if len(new_blst_coef[i]) > 1:
MV.contract(new_blst_coef[i],new_blst_expand[i])
return(new_blst_coef,new_blst_expand)

@staticmethod
def contract(coefs,bases):
dualsort(coefs,bases)
n = len(bases)-1
i = 0
while i < n:
j = i+1
if bases[i] == bases[j]:
coefs[i] += coefs[j]
bases.pop(j)
coefs.pop(j)
n -= 1
else:
i += 1
n = len(coefs)
i = 0
while i < n:
if coefs[i] == ZERO:
coefs.pop(i)
bases.pop(i)
n -= 1
else:
i +=1
return

@staticmethod
def convert(coefs,bases):
mv = MV()
if len(coef) > 0:
nbaserg = list(range(len(base)))
for ibase in nbaserg:
else:
return(mv)

@staticmethod
def set_str_format(str_mode=0):
MV.str_mode = str_mode
return

@staticmethod
def str_rep(mv):
"""
Converts internal representation of a multivector to a string
for outputting.  If lst_mode = 1, str_rep outputs a list of
strings where each string contains one multivector coefficient
concatenated with the corresponding base or blade symbol.

MV.str_mode     Effect
0           Print entire multivector on one line (default)
1           Print each grade on a single line
2           Print each base on a single line
"""

else:
labels = MV.basislabel
else:
mv.compact()
if isinstance(mv.mv[0],int):
value = ''
else:
value = (mv.mv[0][0]).__str__()
value = value.replace(' ','')
dummy = sympy.Dummy('dummy')
j = 0
if x != ZERO:
xstr = (x*dummy).__str__()
xstr = xstr.replace(' ','')
if xstr[0] != '-' and len(value) > 0:
xstr = '+'+xstr
if xstr.find('_dummy') < 2 and xstr[-5:] != '_dummy':
else:
if MV.str_mode == 2:
xstr += '\n'
value += xstr
j += 1
if MV.str_mode == 1:
value += '\n'
if value == '':
value = '0'
#value = value.replace(' ','')
value = value.replace('dummy','1')
return(value)

@staticmethod
def xstr_rep(mv,lst_mode=0):
"""
Converts internal representation of a multivector to a string
for outputing.  If lst_mode = 1, str_rep outputs a list of
strings where each string contains one multivector coefficient
concatenated with the corresponding base or blade symbol.
"""
if lst_mode:
outlst = []
else:
labels = MV.basislabel
else:
value = ''
tmp = []
j = 0
if x != ZERO:
xstr = x.__str__()
if xstr == '+1' or xstr == '1' or xstr == '-1':
if xstr == '+1' or xstr == '1':
xstr = '+'
else:
xstr = '-'
else:
if xstr[0] != '+':
xstr = '+('+xstr+')'
else:
xstr = '+('+xstr[1:]+')'
if MV.str_mode and not lst_mode:
value += value+'\n'
if lst_mode:
tmp.append(value)
j += 1
if lst_mode:
if len(tmp) > 0:
outlst.append(tmp)
value = ''
if not lst_mode:
if len(value) > 1 and value[0] == '+':
value = value[1:]
if len(value) == 0:
value = '0'
else:
value = outlst
return(value)

@staticmethod
def setup(basis,metric='',rframe=False,coords=None,debug=False,offset=0):
"""
MV.setup initializes the MV class by calculating the static
multivector tables required for geometric algebra operations
on multivectors.  See reference 5 section 2 for details on
basis and metric arguments.
"""
MV.is_setup = True
MV.metric_str = False
MV.debug = debug
MV.tables_flg = 0
MV.str_mode  = 0
MV.basisroot = ''
MV.index_offset = offset
if coords == None:
MV.coords = None
else:
MV.coords = tuple(coords)
rframe= True
if type(basis) == str:
basislst = basis.split()
if len(basislst) == 1:
MV.basisroot = basislst[0]
basislst = []
for coord in coords:
basislst.append(MV.basisroot+'_'+str(coord))
MV.define_basis(basislst)
if type(metric) == str:
MV.metric_str = True
if len(metric) > 0:
if metric[0] == '[' and metric[-1] == ']':
tmps = metric[1:-1].split(',')
N = len(tmps)
metric = []
itmp = 0
for tmp in tmps:
xtmp = N*['0']
xtmp[itmp] = tmp
itmp += 1
metric.append(xtmp)
else:
tmps = metric.split(',')
metric = []
for tmp in tmps:
xlst = tmp.split()
xtmp = []
for x in xlst:
xtmp.append(x)
metric.append(xtmp)

MV.define_metric(metric)
MV.multiplication_table()
MV.tables_flg = 1
isym = 0
MV.bvec = []
for name in MV.vbasis:
bvar = MV(value=isym,mvtype='basisvector',mvname=name)
MV.bvec.append(bvar)
isym += 1
if rframe:
MV.define_reciprocal_frame()
MV.I = MV(ONE,'pseudo','I')
MV.ZERO = MV()
Isq = (MV.I*MV.I)()
MV.Iinv = (1/Isq)*MV.I
return MV.bvec

@staticmethod
def set_coords(coords):
MV.coords = coords
return

@staticmethod
def scalar_fct(fct_name):
"""
Create multivector scalar function with name fct_name (string) and
independent variables coords (list of variable).  Default variables are
those associated with each dimension of vector space.
"""
phi = sympy.Function(fct_name)(*MV.coords)
Phi = MV(phi,'scalar')
Phi.name = fct_name
return(Phi)

@staticmethod
def vector_fct(fct_name,vars=''):
"""
Create multivector vector function with name fct_name (string) and
independent variables coords (list of variable).  Default variables are
those associated with each dimension of vector space.
"""
if isinstance(vars,str):
Acoefs = vector_fct(fct_name,MV.coords)
else:
Acoefs =numpy.array(MV.n*[ZERO],dtype=numpy.object)
x = MV.coords
if isinstance(vars,sympy.core.symbol.Symbol):
for icoef in MV.nrg:
Acoefs[icoef] = sympy.Function(fct_name+'__'+\
sympy.galgebra.latex_ex.LatexPrinter.str_basic(x[icoef]))(vars)
else:
for icoef in MV.nrg:
Acoefs[icoef] = sympy.Function(fct_name+'__'+\
sympy.galgebra.latex_ex.LatexPrinter.str_basic(x[icoef]))(*tuple(vars))
A = MV(Acoefs,'vector',fct_name)
return(A)

@staticmethod
def rebase(x,coords,base_name='',debug=False,debug_level=0):
"""
Define curvilinear coordinates for previously defined vector (multivector) space (MV.setup has been run)
with position vector, x, that is a vector function of the independent coordinates, coords (list of
sympy variables equal in length to dimension of vector space), and calculate:

1. Frame (basis) vectors
2. Normalized frame (basis) vectors.
3. Metric tensor
4. Reciprocal frame vectors
5. Reciprocal metric tensor
6. Connection multivectors

The basis vectors are named with the base_name (string) and a subscript derived from the name of each
coordinate.  So that if the base name is 'e' and the coordinated are [r,theta,z] the variable names
of the frame vectors would be e_r, e_theta, and e_z.  For LaTeX output the names of the frame vectors
would be e_{r}, e_{\theta}, and e_{z}.  Everything needed to compute the geometric, outer, and inner
derivatives of multivector functions in curvilinear coordinates is calculated.

If debug is True all the quantities in the above list are output in LaTeX format.

Currently rebase works with cylindrical and spherical coordinates in any dimension.  The limitation is the
ability to automatically simplify complex sympy expressions generated while calculating the quantities in
the above list.  This is why the debug option is included.  The debug_level can equal 0,1,2, or 3 and
determines how far in the list to calculate (input 0 to do the entire list) while debugging.
"""
#Form root names for basis, reciprocal basis, normalized basis, and normalized reciprocal basis

if base_name == '':
base_name = MV.basisroot+'prm'

LaTeX_base = sympy.galgebra.latex_ex.LatexPrinter.extended_symbol(base_name)
bm = '\\bm{'+LaTeX_base+'}'
bmhat = '\\hat{'+bm+'}'
bstr = bmhat+'_{[i_{1},\dots, i_{R}]}'
base_name += 'bm'
base_name_hat = base_name+'hat'

base_name_lst   = []
nbase_name_lst  = []
rbase_name_lst  = []
rnbase_name_lst = []
coords_lst = []

for coord in coords:
coord_str = sympy.galgebra.latex_ex.LatexPrinter.str_basic(coord)
coords_lst.append(coord_str)
base_name_lst.append(base_name+'_'+coord_str)
rbase_name_lst.append(base_name+'__'+coord_str)
nbase_name_lst.append(base_name_hat+'_'+coord_str)
rnbase_name_lst.append(base_name_hat+'__'+coord_str)

if not (MV.n == len(coords) == len(base_name_lst)):
print('rebaseMV inputs not congruent:')
print('MV.n =',MV.n)
print('coords =',coords)
print('bases =',base_name)
sys.exit(1)

if isinstance(x,MV):

#Calculate basis vectors from derivatives of position vector x

bases = numpy.array(MV.n*[ZERO],dtype=numpy.object)
i = 0
for coord in coords:
ei = x.diff(coords[i])
ei.set_name(base_name_lst[i])
bases[i] = ei
i += 1

#Calculate normalizee basis vectors and basis vector magnitudes

if debug:
print('Coordinate Generating Vector')
print(x)
print('Basis Vectors')
for base in bases:
print(base)

else:

#Input basis vectors as N vector fields

bases = x

for (base,name) in zip(bases,base_name_lst):
base.set_name(name)

if debug:
print('Basis Vectors')
for base in bases:
print(base)
if debug_level == 1:
return

if debug_level == 1:
return

#Calculate normalized basis vectors and magnitudes of
#unormalized basis vectors

(nbases,mags) = normalize(bases,nbase_name_lst)

if debug:
print('Magnitudes')
print('\\abs{'+LaTeX_base+'_{i}} = ',mags)
print('Normalized Basis Vectors')
for nbase in nbases:
print(nbase)

g =  numpy.array(MV.n*[MV.n*[ZERO]],dtype=numpy.object)

for irow in MV.nrg:
for icol in MV.nrg:
magsq = sympy.expand((nbases[irow]|nbases[icol])())
g[irow][icol]  = sympy.simplify(sympy.trigsimp(magsq,deep=True,recursive=True))

if debug:
print('Metric $\\hat{g}_{ij} = \\hat{'+LaTeX_base+'}_{i}\\cdot \\hat{'+\ LaTeX_base+'}_{j}$')
print(r'\hat{g}_{ij} =',sympy.galgebra.latex_ex.LaTeX(g))

if debug_level == 2:
return

#Calculate reciprocal normalized basis vectors

rnbases = []

if is_quasi_unit_numpy_array(g):
ibasis = 0
while ibasis < MV.n:
base = g[ibasis][ibasis]*nbases[ibasis]
base.set_name(rnbase_name_lst[ibasis])
base.simplify()
base.trigsimp()
rnbases.append(base)
ibasis += 1
else:
rnbases = reciprocal_frame(nbases,rnbase_name_lst)
ibase = 0
for base in rnbases:
base.simplify()
base.trigsimp()
rnbases[ibase] = base
ibase += 1

if debug:
if debug_level != 0:
sympy.galgebra.latex_ex.MV_format(1)
print('Reciprocal Normalized Basis Vectors')
for rnbase in rnbases:
print(rnbase)

if debug_level == 3:
return

#Calculate components of inverse vectors

Acoef = []

for ibasis in MV.nrg:
evec = numpy.array(MV.n*[ZERO],dtype=numpy.object)
for jbasis in MV.nrg:
evec[jbasis] = (MV.bvec[ibasis]|rnbases[jbasis])()
Acoef.append(evec)

#Calculat metric tensors

gr = numpy.array(MV.n*[MV.n*[ZERO]],dtype=numpy.object)

for irow in MV.nrg:
for icol in MV.nrg:
magsq = sympy.expand((rnbases[irow]|rnbases[icol])())
gr[irow][icol] = sympy.simplify(sympy.trigsimp(magsq,deep=True,recursive=True))

if debug:
print('Metric $\\hat{g}^{ij} = \\hat{'+LaTeX_base+'}^{i}\\cdot \\hat{'+\ LaTeX_base+'}^{j}$')
print(r'\hat{g}^{ij} =',sympy.galgebra.latex_ex.LaTeX(gr))

if debug_level == 4:
return

#Calculate bases and reciprocal bases for curvilinear mulitvectors

MV_bases = [[ONE]]
MV_rbases = [[ONE]]
for index in base_index:
base = build_base(index,nbases)
base.simplify()
base.trigsimp()
rbase = build_base(index,rnbases,True)
rbase.simplify()
rbase.trigsimp()

#Calculate connection multivectors for geometric derivative

MV_connect = [[ZERO]]
ibase = 0
sum = MV()
itheta = 0
for (theta,etheta) in zip(coords,rnbases):
psum = (1/mags[itheta])*etheta*base.diff(theta)
psum.trigsimp()
sum += psum
itheta += 1
sum.simplify()
sum.trigsimp()
ibase += 1

if debug:
print('Curvilinear Bases: $'+bstr+' = '+bmhat+'_{i_{1}}\\W\\dots\\W'+bmhat+'_{i_{R}}$')
ibase = 0
sub_str = ''
for i in index:
sub_str += sympy.galgebra.latex_ex.LatexPrinter.extended_symbol(coords_lst[i])
base_str = bmhat+'_{['+sub_str+']} = '
print(base_str,base)
ibase += 1

if debug_level == 5:
return

#Calculate representation of connection multivectors in curvilinear system

MV_Connect = [[ZERO]]
ibase = 0

ibase = 0
while ibase < nbase:
Cm1 =  numpy.array(m1base*[ZERO],dtype=numpy.object)
Cp1 =  numpy.array(p1base*[ZERO],dtype=numpy.object)
X = C(0)
else:
X = MV()
jbase = 0
while jbase < m1base:
jbase += 1
jbase = 0
while jbase < p1base:
jbase += 1
X.simplify()
X.trigsimp()
ibase += 1
else:
ibase = 0
while ibase < nbase:
Cm1 =  numpy.array(m1base*[ZERO],dtype=numpy.object)
jbase = 0
while jbase < m1base:
jbase += 1
X = MV()
X.mv[MV.n-1] = Cm1
X.simplify()
X.trigsimp()
ibase += 1

base_str = ''
for coord in coords:
base_str += base_name+'_'+sympy.galgebra.latex_ex.LatexPrinter.str_basic(coord)+' '
base_str = base_str[:-1]

old_names = MV.vbasis

MV.setup(base_str,g,True,coords)

MV.curvilinear_flg = True
MV.Connect = MV_Connect
sympy.galgebra.latex_ex.LatexPrinter.latex_bases()
MV.Rframe = numpy.array(MV.n*[ZERO],dtype=numpy.object)
ibasis = 0
while ibasis < MV.n:
base = MV()
jbasis = 0
while jbasis < MV.n:
jbasis += 1
base.scalar_mul_inplace(1/mags[ibasis])
MV.Rframe[ibasis] = base
ibasis += 1

MV.org_basis = []
for ibasis in MV.nrg:
evec = MV(Acoef[ibasis],'vector',old_names[ibasis])
MV.org_basis.append(evec)

if MV.coords[0] == sympy.Symbol('t'):
MV.dedt = []
# I can't fix this no name error because I have no clue what the
# original writer was trying to do.
for coef in dedt_coef:
MV.dedt.append(MV(coef,'vector'))
else:
MV.dedt = None

if debug:
print('Representation of Original Basis Vectors')
for evec in MV.org_basis:
print(evec)

print('Renormalized Reciprocal Vectors '+\
'$\\bfrac{'+bmhat+'^{k}}{\\abs{\\bm{'+LaTeX_base+'}_{k}}}$')

ibasis = 0
while ibasis < MV.n:
c_str = sympy.galgebra.latex_ex.LatexPrinter.extended_symbol(coords_lst[ibasis])
print('\\bfrac{\\bm{\\hat{'+LaTeX_base+\
'}}^{'+c_str+'}}{\\abs{\\bm{'+LaTeX_base+\
'}_{'+c_str+'}}} =',MV.Rframe[ibasis])
ibasis += 1

title_str = 'Connection Multivectors: $C\\lbrc'+bstr+'\\rbrc = '+\ '\\bfrac{'+bmhat+'^{k}}{\\abs{'+bmhat+\ '_{k}}}\\pdiff{'+bstr+'}{\\theta^{k}}$'

print(title_str)
ibase = 0
sub_str = ''
for i in index:
sub_str += sympy.galgebra.latex_ex.LatexPrinter.extended_symbol(coords_lst[i])

base_str = 'C\\lbrc\\hat{'+LaTeX_base+'}_{['+sub_str+']}\\rbrc = '
print(base_str,base)
ibase += 1
return

@staticmethod
"""
Set multivector output to blade representation.
"""
return

@staticmethod
def print_bases():
"""
Set multivector output to base representation.
"""
return

@staticmethod
def multiplication_table():
"""
Calculate geometric product base multiplication table.
See reference 5 section 3 for details.
"""
MV.mtable = []
MV.mtable.append([])
base1 = []
else:
base2 = []
else:
base = base1+base2
(coefs,bases) = MV.reduce_basis(base)
product = MV.convert(coefs,bases)
if MV.debug:
print('Multiplication Table:')
for level1 in MV.mtable:
for level2 in level1:
for level3 in level2:
for mv in level3:
mv.printmv()
return

@staticmethod
def geometric_product(mv1,mv2):
"""
MV.geometric_product(mv1,mv2) calculates the geometric
product the multivectors mv1 and mv2 (mv1*mv2). See
reference 5 section 3.
"""
product = MV()
if isinstance(mv1,MV) and isinstance(mv2,MV):
if xi != ZERO:
if xj != ZERO:
else:
if isinstance(mv1,MV):
product = mv1.scalar_mul(mv2)
else:
product = mv2.scalar_mul(mv1)
return(product)

@staticmethod
"""
Calculate the outer product of a multivector blade
and a vector.  See reference 5 section 5 for details.
"""
w = w12-w21
else:
w = w12+w21
w.name = name
return(w*HALF)

@staticmethod
"""
Calculate basis blades in terms of bases. See reference 5
section 5 for details. Used to convert from blade to base
representation.
"""
MV.btable = []
basis_str = MV.basislabel[1]
tmp = [MV(value=ONE,mvtype='scalar',mvname='1')]
tmp = []
for ibase in range(MV.nbasis[1]):
tmp.append(MV(value=ibase,mvtype='basisvector',mvname=basis_str[ibase]))
tmp = []
name = ''
name += basis_str[i]+'^'
name = name[:-1]
tmp.append(b1Wv2)
MV.btable.append(tmp)
if MV.debug:
print(mv)
return

@staticmethod
"""
Calculate bases in terms of basis blades. See reference 5
section 5 for details. Used to convert from base to blade
representation.
"""
MV.ibtable = []
tmp = [MV(value=ONE,mvtype='scalar',mvname='1')]
tmp = []
for ibase in range(MV.nbasis[1]):
tmp.append(MV(value=ibase,mvtype='basisvector'))
tmp = []
MV.ibtable.append(tmp)
if MV.debug:
mv.printmv()
return

@staticmethod
def outer_product(mv1,mv2):
"""
MV.outer_product(mv1,mv2) calculates the outer (exterior,wedge)
product of the multivectors mv1 and mv2 (mv1^mv2). See
reference 5 section 6.
"""
if type(mv1) == type(MV) and type(mv2) == type(MV):
if mv1.is_scalar() and mv2.is_scalar():
return(mv1*mv2)

if isinstance(mv1,MV) and isinstance(mv2,MV):
product = MV()
pg1pg2 = pg1*pg2
else:
if isinstance(mv1,MV):
product = mv1.scalar_mul(mv2)
if isinstance(mv2,MV):
product = mv2.scalar_mul(mv1)
return(product)

@staticmethod
def inner_product(mv1,mv2,mode='s'):
"""
MV.inner_product(mv1,mv2) calculates the inner

mode = 's' - symmetric (Doran & Lasenby)
mode = 'l' - left contraction (Dorst)
mode = 'r' - right contraction (Dorst)
"""
if type(mv1) == type(MV) and type(mv2) == type(MV):
if mv1.is_scalar() and mv2.is_scalar():
return(mv1*mv2)

if isinstance(mv1,MV) and isinstance(mv2,MV):
product = MV()
if mode == 's':
else:
if mode == 'l':
pg1pg2 = pg1*pg2
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)

@staticmethod
"""
of the multivectors mv1 and mv2 (mv1+mv2).
"""
sum = MV()
if isinstance(mv1,MV) and isinstance(mv2,MV):
for i in MV.n1rg:
if isinstance(mv1.mv[i],numpy.ndarray) and isinstance(mv2.mv[i],numpy.ndarray):
sum.mv[i] = mv1.mv[i]+mv2.mv[i]
else:
if isinstance(mv1.mv[i],numpy.ndarray) and not isinstance(mv2.mv[i],numpy.ndarray):
sum.mv[i] = +mv1.mv[i]
else:
if isinstance(mv2.mv[i],numpy.ndarray) and not isinstance(mv1.mv[i],numpy.ndarray):
sum.mv[i] = +mv2.mv[i]
return(sum)
else:
if isinstance(mv1,MV):
return(mv1+MV(mv2,'scalar'))
else:
return(MV(mv1,'scalar')+mv2)

@staticmethod
def subtraction(mv1,mv2):
"""
MV.subtraction(mv1,mv2) calculates the difference
of the multivectors mv1 and mv2 (mv1-mv2).
"""
diff = MV()
if isinstance(mv1,MV) and isinstance(mv2,MV):
for i in MV.n1rg:
if isinstance(mv1.mv[i],numpy.ndarray) and isinstance(mv2.mv[i],numpy.ndarray):
diff.mv[i] = mv1.mv[i]-mv2.mv[i]
else:
if isinstance(mv1.mv[i],numpy.ndarray) and not isinstance(mv2.mv[i],numpy.ndarray):
diff.mv[i] = +mv1.mv[i]
else:
if not isinstance(mv1.mv[i],numpy.ndarray) and isinstance(mv2.mv[i],numpy.ndarray):
diff.mv[i] = -mv2.mv[i]
return(diff)
else:
if isinstance(mv1,MV):
return(mv1-MV(mv2,'scalar'))
else:
return(MV(mv1,'scalar')-mv2)

@staticmethod
def vdiff(vec,x):
dvec = numpy.array(len(vec)*[ZERO])
ivec = 0
for veci in vec:
dvec[ivec] = sympy.diff(veci,x)
ivec += 1
return(dvec)

@staticmethod
def scalar_to_symbol(scalar):
if isinstance(scalar,MV):
return(scalar.mv[0][0])
if type(scalar) == list:
sym = []
for x in scalar:
sym.append(MV.scalar_to_symbol(x))
return(sym)
return(scalar)

def __init__(self,value='',mvtype='',mvname='',fct=False,vars=None):
"""
Initialization of multivector X. Inputs are as follows

mvtype           value                       result

default         default                  Zero multivector
'basisvector'   int i                    ith basis vector
'basisbivector' int i                    ith basis bivector
'scalar'        symbol x                 scalar of value x
string s
[int i, string s]
'vector'        symbol array A           X.grade(1) = A
string s
string s
'pseudo'        symbol x                 X.grade(n) = x
string s
'spinor'        string s                 spinor with coefficients
s__indices and name sbm

mvname is name of multivector.
If fct is 'True' and MV.coords is defined in MV.setup then a
multivector field of MV.coords is instantiated.
"""

self.name      = mvname
self.mv        = MV.n1*[0]
if mvtype == 'basisvector':
self.mv[1] = numpy.array(MV.nbasis[1]*[ZERO],dtype=numpy.object)
self.mv[1][value] = ONE
if mvtype == 'basisbivector':
self.mv[2] = numpy.array(MV.nbasis[2]*[ZERO],dtype=numpy.object)
self.mv[2][value] = ONE
if mvtype == 'scalar':
if isinstance(value,str):
value = sympy.Symbol(value)
if isinstance(value,int):
value = sympy.Rational(value)
self.mv[0] = numpy.array([value],dtype=numpy.object)
if mvtype == 'pseudo':
if isinstance(value,str):
value = sympy.Symbol(value)
self.mv[MV.n] = numpy.array([value],dtype=numpy.object)
if mvtype == 'vector':
if isinstance(value,str): #Most general vector
symbol_str = ''
for ibase in MV.nrg:
if MV.coords == None:
symbol = value+'__'+str(ibase+MV.index_offset)
symbol_str += symbol+' '
else:
symbol = value+'__'+(MV.coords[ibase]).name
symbol_str += symbol+' '
symbol_lst = sympy.symbols(symbol_str)
self.mv[1] = numpy.array(symbol_lst,dtype=numpy.object)
self.name = value
else:
self.mv[1] = numpy.array(value,dtype=numpy.object)
if isinstance(value,str): #Most general grade-2 multivector
if value != '':
symbol_str = ''
for base in MV.basis[2]:
symbol = value+MV.construct_index(base)
symbol_str += symbol+' '
symbol_lst = sympy.symbols(symbol_str)
self.mv[2] = numpy.array(symbol_lst,dtype=numpy.object)
else:
self.mv[2] = numpy.array(value,dtype=numpy.object)
coefs  = value[1]
if isinstance(coefs,str): #Most general pure grade multivector
base_symbol = coefs
coefs = []
self.mv[0] = numpy.array([sympy.Symbol(base_symbol)],dtype=numpy.object)
else:
for base in bases:
coef = base_symbol+MV.construct_index(base)
coef = sympy.Symbol(coef)
coefs.append(coef)
else:
if mvtype == 'base':
self.mv[value[0]] = numpy.array(MV.nbasis[value[0]]*[ZERO],dtype=numpy.object)
self.mv[value[0]][value[1]] = ONE
if mvtype == 'spinor':
if isinstance(value,str): #Most general spinor
symbol_str = ''
symbol_lst = []
symbol = sympy.Symbol(value+MV.construct_index(base))
symbol_lst.append(symbol)
else:
self.mv[0] = numpy.array([sympy.Symbol(value)],dtype=numpy.object)
self.name = value+'bm'
if isinstance(value,str) and mvtype == '': #Most general multivector
if value != '':
symbol_str = ''
symbol_lst = []
symbol = sympy.Symbol(value+MV.construct_index(base))
symbol_lst.append(symbol)
else:
self.mv[0] = numpy.array([sympy.Symbol(value)],dtype=numpy.object)
self.name = value+'bm'
if fct:
if vars is not None:
vars = tuple(vars)
coef = sympy.galgebra.latex_ex.LatexPrinter.str_basic(self.mv[0][0])
if vars == None and MV.coords is not None:
self.mv[0]= numpy.array([sympy.Function(coef)(*MV.coords)],dtype=numpy.object)
else:
self.mv[0]= numpy.array([sympy.Function(coef)(*vars)],dtype=numpy.object)
else:
if vars == None and MV.coords is not None:
else:

@staticmethod
def construct_index(base):
index_str = ''
if len(base) == 0:
return('')
if MV.coords == None:
for ix in base:
index_str += str(ix+MV.index_offset)
else:
for ix in base:
index_str += (MV.coords[ix]).name
return('__'+index_str)

def set_name(self,namestr):
self.name = namestr
return

"""
"""
for i in range(MV.n,-1,-1):
if isinstance(self.mv[i],numpy.ndarray):
return(i)
return(-1)

@staticmethod
def coord(xname,offset=0):
xi_str = ''
for i in MV.nrg:
xi_str += xname+str(i+offset)+' '
xi = sympy.symbols(xi_str)
x = MV(xi,'vector')
return(x)

def x(self,i):
if isint(self.mv[1]):
return(ZERO)
return(self.mv[1][i])

return

@staticmethod
def named(mvname,value='',mvtype=''):
name = mvname
tmp = MV(value=value,mvtype=mvtype,mvname=name)
setattr(sys.modules[__name__],name,tmp)
return

@staticmethod
def printnm(tpl):
for a in tpl:
print(a.name,' =',a.mv)
return

def __str__(self):
return(MV.str_rep(self))

def printmv(self,name=''):
title = ''
if name:
title += name+' = '
else:
if self.name:
title += self.name+' = '
print(title+MV.str_rep(self))
return

else:
return

"""
X.add_in_place(mv) increments multivector X by multivector
mv.
"""
for i in MV.n1rg:
if isinstance(self.mv[i],numpy.ndarray) and isinstance(mv.mv[i],numpy.ndarray):
self.mv[i] += mv.mv[i]
else:
if not isinstance(self.mv[i],numpy.ndarray) and isinstance(mv.mv[i],numpy.ndarray):
self.mv[i] = +mv.mv[i]
return

def sub_in_place(self,mv):
"""
X.sub_in_place(mv) decrements multivector X by multivector
mv.
"""
for i in MV.n1rg:
if isinstance(self.mv[i],numpy.ndarray) and isinstance(mv.mv[i],numpy.ndarray):
self.mv[i] -= mv.mv[i]
else:
if not isinstance(self.mv[i],numpy.ndarray) and isinstance(mv.mv[i],numpy.ndarray):
self.mv[i] = +mv.mv[i]
return

def __pos__(self):
p = MV()
for i in MV.n1rg:
if isinstance(self.mv[i],numpy.ndarray):
p.mv[i] = +self.mv[i]
return(p)

def __neg__(self):
n = MV()
for i in MV.n1rg:
if isinstance(self.mv[i],numpy.ndarray):
n.mv[i] = -self.mv[i]
return(n)

return

def __sub_ab__(self,mv):
self.sub_in_place(mv)
return

def __sub__(self,mv):
"""See MV.subtraction(self,mv)"""
return(MV.subtraction(self,mv))

def __rsub__(self,mv):
"""See MV.subtraction(mv,self)"""
return(MV.subtraction(mv,self))

def __xor__(self,mv):
"""See MV.outer_product(self,mv)"""
return(MV.outer_product(self,mv))

def __pow__(self,mv):
"""See MV.outer_product(self,mv)"""
return(MV.outer_product(self,mv))

def __rxor__(self,mv):
"""See MV.outer_product(mv,self)"""
return(MV.outer_product(mv,self))

def __or__(self,mv):
"""See MV.inner_product(self,mv)"""
return(MV.inner_product(self,mv))

def __ror__(self,mv):
"""See MV.inner_product(mv,self)"""
return(MV.inner_product(mv,self))

def __lt__(self,mv):
"""See MV.inner_product(self,mv)"""
return(MV.inner_product(self,mv,'l'))

def __lshift__(self,mv):
"""See MV.inner_product(self,mv)"""
return(MV.inner_product(self,mv,'l'))

def __rlshift__(self,mv):
"""See MV.inner_product(self,mv)"""
return(MV.inner_product(mv,self,'l'))

def lc(self,mv):
return(MV.inner_product(self,mv,'l'))

def __gt__(self,mv):
"""See MV.inner_product(self,mv)"""
return(MV.inner_product(self,mv,'r'))

def __rshift__(self,mv):
"""See MV.inner_product(self,mv)"""
return(MV.inner_product(self,mv,'r'))

def __rrshift__(self,mv):
"""See MV.inner_product(self,mv)"""
return(MV.inner_product(mv,self,'r'))

def rc(self,mv):
return(MV.inner_product(self,mv,'r'))

def scalar_mul(self,c):
"""
Y = X.scalar_mul(c), multiply multivector X by scalar c and return
result.
"""
mv = MV()
for i in MV.n1rg:
if isinstance(self.mv[i],numpy.ndarray):
#print self.mv[i]
#print c,type(c)
mv.mv[i] = self.mv[i]*c
return(mv)

def scalar_mul_inplace(self,c):
"""
X.scalar_mul_inplace(c), multiply multivector X by scalar c and save
result in X.
"""
for i in MV.n1rg:
if isinstance(self.mv[i],numpy.ndarray):
self.mv[i] = self.mv[i]*c
return

def __mul__(self,mv):
"""See MV.geometric_product(self,mv)"""
return(MV.geometric_product(self,mv))

def __rmul__(self,mv):
"""See MV.geometric_product(mv,self)"""
return(MV.geometric_product(mv,self))

def __div__(self,scalar):
div = MV()
for i in MV.n1rg:
if isinstance(self.mv[i],numpy.ndarray):
div.mv[i] = self.mv[i]/scalar
return(div)

def __div_ab__(self,scalar):
for i in MV.n1rg:
if isinstance(self.mv[i],numpy.ndarray):
self.mv[i] /= scalar
return

"""
X(i,j) returns symbol in ith grade and jth base or blade of
multivector X.
"""
return(ZERO)

@staticmethod
def equal(mv1,mv2):
mv1.compact()
if isinstance(mv2,MV):
mv2.compact()
if not isinstance(mv2,MV) and pure_grade != 0:
return(False)
if not isinstance(mv2,MV) and pure_grade == 0:
if isinstance(mv1.mv[0],int):
return(mv2 == 0)
else:
return(mv1.mv[0][0] == mv2)
for (mvi,mvj) in zip(mv1.mv,mv2.mv):
if isint(mvi) ^ isint(mvj):
return(False)
if isinstance(mvi,numpy.ndarray) and isinstance(mvj,numpy.ndarray):
for (x,y) in zip(mvi,mvj):
if x != y:
return(False)
return(True)

def __eq__(self,mv):
return(MV.equal(self,mv))

def copy(self,sub=0):
"""
Y = X.copy(), make a deep copy of multivector X in multivector
Y so that Y can be modified without affecting X.
"""
cpy = MV()
cpy.name      = self.name
for i in MV.n1rg:
if sub:
if isinstance(self.mv[i],numpy.ndarray):
cpy.mv[i] = -self.mv[i]
else:
if isinstance(self.mv[i],numpy.ndarray):
cpy.mv[i] = +self.mv[i]
return(cpy)

return
if isinstance(base,numpy.ndarray):
else:
ibase = base
if coef == ZERO:
return
return

"""
to blade representation.  See reference 5 section 5.
"""
return
if (not coef == ZERO):
return

"""
to base representation.  See reference 5 section 5.
"""
return
if (not coef == ZERO):
return

def project(self,r):
"""
Grade projection operator. For multivector X, X.project(r)
returns multivector of grade r components of X if r is an
integer. If r is a multivector X.project(r) returns a
multivector consisting of the grade of X for which r has non-
zero grades. For example if X is a general multivector and
r is a general spinor then X.project(r) will return the even
"""
if isinstance(r,int):
if r > MV.n:
if not isinstance(self.mv[r],numpy.ndarray):
if isinstance(r,MV):
proj = MV()
for i in MV.n1rg:
if not isinstance(r.mv[i],int):
proj.mv[i] = self.mv[i]
return(proj)
return(None)

def even(self):
"""
Even grade projection operator. For multivector X, X.even()
returns multivector of even grade components of X.
"""

def odd(self):
"""
Odd grade projection operator. For multivector X, X.odd()
returns multivector of odd grade components of X.
"""

def rev(self):
"""
Revision operator. For multivector X, X.rev()
returns reversed multivector of X.
"""
revmv = MV()
else:
return(revmv)

cse_lst = []
return(cse_lst)

div_lst = []
return(div_lst)

# don't know which one is correctly named

def div(self):

def collect(self,lst):
"""
Applies sympy collect function to each component
of multivector.
"""
return

def sqrfree(self,lst):
"""
Applies sympy sqrfree function to each component
of multivector.
"""
return

def flatten(self):
flst = []
else:
flst.append(coef)
return(flst)

def subs(self,*args):
X = MV()
return(X)

def sub_mv(self,mv1,mv2):
mv1_flat = mv1.flatten()
mv2_flat = mv2.flatten()
self.sub_scalar(mv1_flat,mv2_flat)
return

def sub_scalar(self,expr1,expr2):
if (isinstance(expr1,list) and isinstance(expr2,list)) or \
(isinstance(expr1,tuple) and isinstance(expr2,tuple)):
for (var1,var2) in zip(expr1,expr2):
self.sub_scalar(var1,var2)
else:
if expr1 != ZERO:
return

def simplify(self):
"""
Applies sympy simplify function
to each component of multivector.
"""
return

def trigsimp(self):
"""
Applies sympy trigsimp function
to each component of multivector
using the recursive option.
"""
return

def cancel(self):
"""
Applies sympy cancel function
to each component of multivector.
"""
return

def expand(self):
"""
Applies sympy expand function to each component
of multivector.
"""
return

def is_pure(self):
self.compact()
for i in MV.n1rg:
if isinstance(self.mv[i],numpy.ndarray):
for base in self.mv[i]:
if base != 0:
break
return(-1)
return(0)

def is_scalar(self):
if self.is_pure() == 0:
return(True)
return(False)

def compact(self):
"""
Convert zero numpy arrays to single integer zero place holder
in grade list for instantiated multivector. For example if
numpy array of grade one components is a zero array then replace
with single integer equal to zero.
"""
for i in MV.n1rg:
if not isint(self.mv[i]):
zero_flg = True
for x in self.mv[i]:
if x != 0:
zero_flg = False
break
if zero_flg:
self.mv[i] = 0
return

def diff(self,x):
"""
Calculate partial derivative of multivector with respect to
argument x.
"""
D = MV()
return(D)

def ddt(self):
if MV.coords[0] != sympy.Symbol('t'):
return(MV())
dxdt = self.diff(MV.coords[0])
for ibase in MV.nrg:
dxdt += self.mv[1][ibase]*MV.dedt[ibase]
#dxdt.simplify()
#dxdt.trigsimp()
return(dxdt)

"""
Calculate geometric (grad) derivative of multivector function
"""
D = []
dD = MV()
for theta in MV.coords:
D.append(self.diff(theta))
if MV.curvilinear_flg:
recp = MV.Rframe
else:
recp = MV.brecp
for (rbase,iD) in zip(recp,D):
if type(coefs) != int:
return(dD)

"""
Calculate outer (exterior,curl) derivative of multivector function.
"""
D = []
dD = MV()
for ix in MV.coords:
D.append(self.diff(ix))
if MV.curvilinear_flg:
recp = MV.Rframe
else:
recp = MV.brecp
for (irbase,iD) in zip(recp,D):
if type(coefs) != int:
return(dD)

def curl(self):

"""
Calculate inner (interior,div) derivative of multivector function.
"""
D = []
dD = MV()
for ix in MV.coords:
D.append(self.diff(ix))
if MV.curvilinear_flg:
recp = MV.Rframe
else:
recp = MV.brecp
for (irbase,iD) in zip(recp,D):
if type(coefs) != int:
return(dD)

def mag2(self):
"""
Calculate scalar component of square of multivector.
"""
return((self|self)())

def Name(self):
"""
Get LaTeX name of multivector.
"""
return(sympy.galgebra.latex_ex.LatexPrinter.extended_symbol(self.name))

[docs]def set_names(var_lst,var_str):
"""
Set the names of a list of multivectors (var_lst) for a space delimited
string (var_str) containing the names.
"""
var_str_lst = var_str.split()
if len(var_lst) == len(var_str_lst):
for (var,var_name) in zip(var_lst,var_str_lst):
var.set_name(var_name)
return
sys.stderr.write('Error in set_names. Lists incongruent!\n')
sys.exit()
return

[docs]def reciprocal_frame(vlst,names=''):
"""
Calculate reciprocal frame of list (vlst) of vectors.  If desired name each
vector in list of reciprocal vectors with names in space delimited string
(names).
"""
E = vlst[0]
recp = []
if type(names) != str:
name_lst = names
else:
if names != '':
name_lst = names.split()
for i in range(1,MV.n):
E = E^vlst[i]
for i in range(MV.n):
tmp = ONE
if i%2 != 0:
tmp = -ONE
for j in range(MV.n):
if i != j:
tmp = tmp^vlst[j]
tmp = tmp*E
recp.append(tmp)
Esq = sympy.trigsimp(E.mag2(),deep=True,recursive=True)
print(Esq)
print(sympy.simplify(Esq))
Esq_inv = ONE/Esq
i = 0
for i in range(MV.n):
recp[i].trigsimp()
recp[i] = recp[i]*Esq_inv
if names != '':
recp[i].set_name(name_lst[i])
i += 1
return(recp)

def S(value):
return(MV(value,'scalar'))