#!/usr/bin/python
"""
The module GA 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.
"""
import sys
import numpy
import 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
from sympy.core import Add as add_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):
"""
Returns string equivalent metric tensor for signature (p, q).
"""
n = p + q
D = []
for i in xrange(p):
D.append((i*'0 ' +'1 '+ (n-i-1)*'0 ')[:-1])
for i in xrange(p,n):
D.append((i*'0 ' +'-1 '+ (n-i-1)*'0 ')[:-1])
return ','.join(D)
[docs]def arbitrary_metric(n):
"""
Returns string equivalent metric tensor for arbitrary signature.
"""
return ','.join(n*[(n*'# ')[:-1]])
[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(useinstead="sympy.symbols()", issue=3379,
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 = range(len(lst1))
_indices.sort(key=lst2.__getitem__)
lst1[:] = map(lst1.__getitem__, _indices)
lst2[:] = 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)
grade = len(base)
if grade == 1:
if k < base[0]:
return(1, [k, base[0]])
else:
return(-1, [base[0], k])
ilo = 0
ihi = grade - 1
if k < base[0]:
return(1, [k] + base)
if k > base[ihi]:
if grade % 2 == 0:
return(1, base + [k])
else:
return(-1, base + [k])
imid = ihi + ilo
if grade == 2:
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
def pad_zeros(value, n):
"""
Pad list with zeros to length n. If length is > n
truncate list. Return padded list.
"""
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 = range(MV.n)
MV.n1 = MV.n + 1
MV.n1rg = range(MV.n1)
MV.npow = 2**MV.n
MV.index = 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)
for igrade in range(1, MV.n + 1):
tmp = comb(MV.n, igrade)
MV.gabasis += [tmp]
ntmp = len(tmp)
MV.nbasis[igrade] = ntmp
MV.basis[igrade] = tmp
gradelabels = []
gradelabel_lst = []
for i in range(ntmp):
tmp_lst = []
bstr = ''
for j in tmp[i]:
bstr += MV.vbasis[j]
tmp_lst.append(MV.vbasis[j])
gradelabel_lst.append(tmp_lst)
gradelabels.append(bstr)
MV.basislabel_lst.append(gradelabel_lst)
MV.basislabel[igrade] = gradelabels
MV.basis_map = [{'':0}]
igrade = 1
while igrade <= MV.n:
tmpdict = {}
bases = MV.gabasis[igrade]
ibases = 0
for base in bases:
tmpdict[str(base)] = ibases
ibases += 1
MV.basis_map.append(tmpdict)
igrade += 1
if MV.debug:
print 'basis strings =', MV.vbasis
print 'basis symbols =', MV.vsyms
print 'basis labels =', MV.basislabel
print 'basis =', MV.basis
print 'grades =', MV.nbasis
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()
mv.bladeflg = 0
for igrade in MV.n1rg:
coef = coefs[igrade]
base = bases[igrade]
if len(coef) > 0:
nbases = MV.nbasis[igrade]
mv.mv[igrade] = numpy.array(nbases*[ZERO], dtype=numpy.object)
nbaserg = range(len(base))
for ibase in nbaserg:
if igrade > 0:
k = MV.basis[igrade].index(base[ibase])
mv.mv[igrade][k] = coef[ibase]
else:
mv.mv[igrade] = numpy.array(
[coef[0]], dtype=numpy.object)
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
"""
if MV.bladeprint:
mv.convert_to_blades()
labels = MV.bladelabel
else:
if not mv.bladeflg:
labels = MV.basislabel
else:
labels = MV.bladelabel
mv.compact()
if isinstance(mv.mv[0], int):
value = ''
else:
value = (mv.mv[0][0]).__str__()
value = value.replace(' ', '')
dummy = sympy.Dummy('dummy')
for igrade in MV.n1rg[1:]:
if isinstance(mv.mv[igrade], numpy.ndarray):
j = 0
for x in mv.mv[igrade]:
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':
xstr = xstr.replace(
'_dummy*', '') + '*' + labels[igrade][j]
else:
xstr = xstr.replace('_dummy', labels[igrade][j])
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 = []
if MV.bladeprint:
mv.convert_to_blades()
labels = MV.bladelabel
else:
if not mv.bladeflg:
labels = MV.basislabel
else:
labels = MV.bladelabel
value = ''
for igrade in MV.n1rg:
tmp = []
if isinstance(mv.mv[igrade], numpy.ndarray):
j = 0
for x in mv.mv[igrade]:
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:] + ')'
value += xstr + labels[igrade][j]
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.bladeprint = 0
MV.tables_flg = 0
MV.str_mode = 0
MV.basisroot = ''
MV.index_offset = offset
if coords is 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.blade_table()
MV.inverse_blade_table()
MV.tables_flg = 1
isym = 0
MV.bvec = []
for name in MV.vbasis:
bvar = MV(value=isym, mvtype='basisvector', mvname=name)
bvar.bladeflg = 1
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]]
igrade = 1
while igrade <= MV.n:
base_index = MV.gabasis[igrade]
grade_bases = []
rgrade_bases = []
for index in base_index:
base = build_base(index, nbases)
base.simplify()
base.trigsimp()
rbase = build_base(index, rnbases, True)
rbase.simplify()
rbase.trigsimp()
grade_bases.append(base)
rgrade_bases.append(rbase)
igrade += 1
MV_bases.append(grade_bases)
MV_rbases.append(rgrade_bases)
#Calculate connection multivectors for geometric derivative
MV_connect = [[ZERO]]
igrade = 1
while igrade <= MV.n:
grade_connect = []
ibase = 0
for base in MV_bases[igrade]:
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()
grade_connect.append(sum)
ibase += 1
MV_connect.append(grade_connect)
igrade += 1
if debug:
print 'Curvilinear Bases: $' + bstr + ' = ' + bmhat + \
'_{i_{1}}\\W\\dots\\W' + bmhat + '_{i_{R}}$'
igrade = 1
for grade in MV_bases[1:]:
ibase = 0
for base in grade:
index = MV.gabasis[igrade][ibase]
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
igrade += 1
if debug_level == 5:
return
#Calculate representation of connection multivectors in curvilinear system
MV_Connect = [[ZERO]]
igrade = 1
while igrade <= MV.n:
grade_connect = []
ibase = 0
nbase = len(MV_bases[igrade])
if igrade < MV.n:
ibase = 0
p1base = len(MV_bases[igrade + 1])
m1base = len(MV_bases[igrade - 1])
while ibase < nbase:
Cm1 = numpy.array(m1base*[ZERO], dtype=numpy.object)
Cp1 = numpy.array(p1base*[ZERO], dtype=numpy.object)
C = MV_connect[igrade][ibase]
if igrade == 1:
X = C(0)
else:
X = MV()
jbase = 0
while jbase < m1base:
Cm1[jbase] = sympy.trigsimp((MV.inner_product(MV_rbases[igrade - 1][jbase], C))(), deep=True, recursive=True)
jbase += 1
jbase = 0
while jbase < p1base:
Cp1[jbase] = sympy.trigsimp((MV.inner_product(MV_rbases[igrade + 1][jbase], C))(), deep=True, recursive=True)
jbase += 1
X += MV(
(igrade - 1, Cm1), 'grade') + MV((igrade + 1, Cp1), 'grade')
X.simplify()
X.trigsimp()
grade_connect.append(X)
ibase += 1
else:
ibase = 0
m1base = len(MV_bases[igrade - 1])
while ibase < nbase:
Cm1 = numpy.array(m1base*[ZERO], dtype=numpy.object)
C = MV_connect[igrade][ibase]
jbase = 0
while jbase < m1base:
Cm1[jbase] = sympy.trigsimp((MV.inner_product(MV_rbases[igrade - 1][jbase], C))(), deep=True, recursive=True)
jbase += 1
X = MV()
X.mv[MV.n - 1] = Cm1
X.simplify()
X.trigsimp()
grade_connect.append(X)
ibase += 1
MV_Connect.append(grade_connect)
igrade += 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:
base.add_in_place(gr[ibasis][jbasis]*MV.bvec[jbasis])
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
igrade = 1
for grade in MV.Connect[1:]:
ibase = 0
for base in grade:
index = MV.gabasis[igrade][ibase]
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
igrade += 1
return
@staticmethod
def print_blades():
"""
Set multivector output to blade representation.
"""
MV.bladeprint = 1
return
@staticmethod
def print_bases():
"""
Set multivector output to base representation.
"""
MV.bladeprint = 0
return
@staticmethod
def multiplication_table():
"""
Calculate geometric product base multiplication table.
See reference 5 section 3 for details.
"""
MV.mtable = []
for igrade in MV.n1rg:
MV.mtable.append([])
for ibase in range(MV.nbasis[igrade]):
MV.mtable[igrade].append([])
if igrade == 0:
base1 = []
else:
base1 = MV.basis[igrade][ibase]
for jgrade in MV.n1rg:
MV.mtable[igrade][ibase].append([])
for jbase in range(MV.nbasis[jgrade]):
if jgrade == 0:
base2 = []
else:
base2 = MV.basis[jgrade][jbase]
base = base1 + base2
(coefs, bases) = MV.reduce_basis(base)
product = MV.convert(coefs, bases)
product.name = '(' + MV.basislabel[igrade][ibase] + \
')(' + MV.basislabel[jgrade][jbase] + ')'
MV.mtable[igrade][ibase][jgrade].append(product)
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):
bladeflg1 = mv1.bladeflg
bladeflg2 = mv2.bladeflg
if bladeflg1:
mv1.convert_from_blades()
if bladeflg2:
mv2.convert_from_blades()
for igrade in MV.n1rg:
gradei = mv1.mv[igrade]
if isinstance(gradei, numpy.ndarray):
for ibase in range(MV.nbasis[igrade]):
xi = gradei[ibase]
if xi != ZERO:
for jgrade in MV.n1rg:
gradej = mv2.mv[jgrade]
if isinstance(gradej, numpy.ndarray):
for jbase in range(MV.nbasis[jgrade]):
xj = gradej[jbase]
if xj != ZERO:
xixj = MV.mtable[igrade][ibase][jgrade][jbase].scalar_mul(xi*xj)
product.add_in_place(xixj)
product.bladeflg = 0
if bladeflg1:
mv1.convert_to_blades()
if bladeflg2:
mv1.convert_to_blades()
if bladeflg1 and bladeflg2:
product.convert_to_blades()
else:
if isinstance(mv1, MV):
product = mv1.scalar_mul(mv2)
else:
product = mv2.scalar_mul(mv1)
return(product)
@staticmethod
def wedge(igrade1, blade1, vector2, name=''):
"""
Calculate the outer product of a multivector blade
and a vector. See reference 5 section 5 for details.
"""
w12 = blade1*vector2
w21 = vector2*blade1
if igrade1 % 2 != 0:
w = w12 - w21
else:
w = w12 + w21
w.name = name
return(w*HALF)
@staticmethod
def blade_table():
"""
Calculate basis blades in terms of bases. See reference 5
section 5 for details. Used to convert from blade to base
representation.
"""
MV.btable = []
MV.bladelabel = []
basis_str = MV.basislabel[1]
for igrade in MV.n1rg:
MV.bladelabel.append([])
if igrade == 0:
MV.bladelabel[0].append('1')
tmp = [MV(value=ONE, mvtype='scalar', mvname='1')]
if igrade == 1:
tmp = []
for ibase in range(MV.nbasis[1]):
MV.bladelabel[1].append(basis_str[ibase])
tmp.append(MV(value=ibase,
mvtype='basisvector', mvname=basis_str[ibase]))
if igrade >= 2:
tmp = []
basis = MV.basis[igrade]
for blade in basis:
name = ''
for i in blade:
name += basis_str[i] + '^'
name = name[:-1]
MV.bladelabel[igrade].append(name)
lblade = MV.basis[igrade - 1].index(blade[:-1])
rblade = blade[-1]
igrade1 = igrade - 1
blade1 = MV.btable[igrade1][lblade]
vector2 = MV.btable[1][rblade]
b1Wv2 = MV.wedge(igrade1, blade1, vector2, name)
tmp.append(b1Wv2)
MV.btable.append(tmp)
if MV.debug:
print 'Blade Tabel:'
for grade in MV.btable:
for mv in grade:
print mv
print 'Blade Labels:'
print MV.bladelabel
return
@staticmethod
def inverse_blade_table():
"""
Calculate bases in terms of basis blades. See reference 5
section 5 for details. Used to convert from base to blade
representation.
"""
MV.ibtable = []
for igrade in MV.n1rg:
if igrade == 0:
tmp = [MV(value=ONE, mvtype='scalar', mvname='1')]
if igrade == 1:
tmp = []
for ibase in range(MV.nbasis[1]):
tmp.append(MV(value=ibase, mvtype='basisvector'))
if igrade >= 2:
tmp = []
iblade = 0
for blade in MV.btable[igrade]:
invblade = MV()
for i in range(igrade - 1):
invblade.mv[i] = -blade.mv[i]
invblade.mv[igrade] = +blade.mv[igrade]
invblade.bladeflg = 1
if igrade >= 4:
jgrade = igrade - 2
while jgrade > 1:
for ibase in range(MV.nbasis[jgrade]):
invblade.substitute_base(
jgrade, ibase, MV.ibtable[jgrade][ibase])
jgrade -= 2
invblade.name = MV.basislabel[igrade][iblade]
iblade += 1
tmp.append(invblade)
MV.ibtable.append(tmp)
if MV.debug:
print 'Inverse Blade Tabel:'
for grade in MV.ibtable:
for mv in grade:
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()
product.bladeflg = 1
mv1.convert_to_blades()
mv2.convert_to_blades()
for igrade1 in MV.n1rg:
if isinstance(mv1.mv[igrade1], numpy.ndarray):
pg1 = mv1.project(igrade1)
for igrade2 in MV.n1rg:
igrade = igrade1 + igrade2
if igrade <= MV.n:
if isinstance(mv2.mv[igrade2], numpy.ndarray):
pg2 = mv2.project(igrade2)
pg1pg2 = pg1*pg2
product.add_in_place(pg1pg2.project(igrade))
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()
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)
@staticmethod
def addition(mv1, mv2):
"""
MV.addition(mv1,mv2) calculates the sum
of the multivectors mv1 and mv2 (mv1+mv2).
"""
sum = MV()
if isinstance(mv1, MV) and isinstance(mv2, MV):
if mv1.bladeflg or mv2.bladeflg:
mv1.convert_to_blades()
mv2.convert_to_blades()
sum.bladeflg = 1
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):
if mv1.bladeflg or mv2.bladeflg:
mv1.convert_to_blades()
mv2.convert_to_blades()
diff.bladeflg = 1
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
'grade' [int i, symbol array A] X.grade(i) = A
[int i, string s]
'vector' symbol array A X.grade(1) = A
string s
'grade2' symbol array A X.grade(2) = A
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]
self.bladeflg = 0 # 1 for blade expansion
self.puregrade = 1
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 is 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:
value = MV.pad_zeros(value, MV.nbasis[1])
self.mv[1] = numpy.array(value, dtype=numpy.object)
if mvtype == 'grade2':
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:
value = MV.pad_zeros(value, MV.nbasis[2])
self.mv[2] = numpy.array(value, dtype=numpy.object)
if mvtype == 'grade':
igrade = value[0]
coefs = value[1]
if isinstance(coefs, str): # Most general pure grade multivector
base_symbol = coefs
coefs = []
bases = MV.basis[igrade]
if igrade == 0:
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)
self.mv[igrade] = numpy.array(coefs, dtype=numpy.object)
else:
self.mv[igrade] = coefs
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
for grade in MV.n1rg:
if grade % 2 == 0:
symbol_str = ''
if grade != 0:
symbol_lst = []
for base in MV.basis[grade]:
symbol = sympy.Symbol(
value + MV.construct_index(base))
symbol_lst.append(symbol)
self.mv[grade] = numpy.array(
symbol_lst, dtype=numpy.object)
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 != '':
for grade in MV.n1rg:
symbol_str = ''
if grade != 0:
symbol_lst = []
for base in MV.basis[grade]:
symbol = sympy.Symbol(
value + MV.construct_index(base))
symbol_lst.append(symbol)
self.mv[grade] = numpy.array(
symbol_lst, dtype=numpy.object)
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)
for grade in MV.n1rg:
if not isinstance(self.mv[grade], int):
if grade == 0:
coef = sympy.galgebra.latex_ex.LatexPrinter.str_basic(
self.mv[0][0])
if vars is 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:
for base in range(MV.nbasis[grade]):
coef = sympy.galgebra.latex_ex.LatexPrinter.str_basic(self.mv[grade][base])
if vars is None and MV.coords is not None:
self.mv[grade][
base] = sympy.Function(coef)(*MV.coords)
else:
self.mv[
grade][base] = sympy.Function(coef)(*vars)
@staticmethod
def construct_index(base):
index_str = ''
if len(base) == 0:
return('')
if MV.coords is 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
def max_grade(self):
"""
X.max_grade() is maximum grade of non-zero grades of X.
"""
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])
def set_coef(self, grade, base, value):
if isinstance(self.mv[grade], int):
self.mv[grade] = numpy.array(
MV.nbasis[grade]*[ZERO], dtype=numpy.object)
self.mv[grade][base] = value
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
def set_value(self, igrade, ibase, value):
if isinstance(self.mv[igrade], numpy.ndarray):
self.mv[igrade][ibase] = value
else:
self.mv[igrade] = numpy.array(
MV.nbasis[igrade]*[ZERO], dtype=numpy.object)
self.mv[igrade][ibase] = value
return
def add_in_place(self, mv):
"""
X.add_in_place(mv) increments multivector X by multivector
mv.
"""
if self.bladeflg or mv.bladeflg:
self.convert_to_blades()
mv.convert_to_blades()
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.
"""
if self.bladeflg or mv.bladeflg:
self.convert_to_blades()
mv.convert_to_blades()
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()
p.bladeflg = self.bladeflg
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()
n.bladeflg = self.bladeflg
for i in MV.n1rg:
if isinstance(self.mv[i], numpy.ndarray):
n.mv[i] = -self.mv[i]
return(n)
def __add_ab__(self, mv):
self.add_in_place(mv)
return
def __sub_ab__(self, mv):
self.sub_in_place(mv)
return
def __add__(self, mv):
"""See MV.addition(self,mv)"""
return(MV.addition(self, mv))
def __radd__(self, mv):
"""See MV.addition(mv,self)"""
return(MV.addition(mv, self))
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()
mv.bladeflg = self.bladeflg
mv.puregrade = self.puregrade
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()
div.bladeflg = self.bladeflg
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
def __call__(self, igrade=0, ibase=0):
"""
X(i,j) returns symbol in ith grade and jth base or blade of
multivector X.
"""
if not isinstance(self.mv[igrade], numpy.ndarray):
return(ZERO)
return(self.mv[igrade][ibase])
@staticmethod
def equal(mv1, mv2):
mv1.compact()
if isinstance(mv2, MV):
mv2.compact()
pure_grade = mv1.is_pure()
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
cpy.bladeflg = self.bladeflg
cpy.puregrade = self.puregrade
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)
def substitute_base(self, igrade, base, mv):
if not isinstance(self.mv[igrade], numpy.ndarray):
return
if isinstance(base, numpy.ndarray):
ibase = MV.basis[igrade].index(base)
else:
ibase = base
coef = self.mv[igrade][ibase]
if coef == ZERO:
return
self.mv[igrade][ibase] = ZERO
self.add_in_place(mv*coef)
return
def convert_to_blades(self):
"""
X.convert_to_blades(), inplace convert base representation
to blade representation. See reference 5 section 5.
"""
if self.bladeflg:
return
self.bladeflg = 1
for igrade in range(2, MV.n1):
if isinstance(self.mv[igrade], numpy.ndarray):
for ibase in range(MV.nbasis[igrade]):
coef = self.mv[igrade][ibase]
if (not coef == ZERO):
self.mv[igrade][ibase] = ZERO
self.add_in_place(MV.ibtable[igrade][ibase]*coef)
return
def convert_from_blades(self):
"""
X.convert_from_blades(), inplace convert blade representation
to base representation. See reference 5 section 5.
"""
if not self.bladeflg:
return
self.bladeflg = 0
for igrade in range(2, MV.n1):
if isinstance(self.mv[igrade], numpy.ndarray):
for ibase in range(MV.nbasis[igrade]):
coef = self.mv[igrade][ibase]
if (not coef == ZERO):
self.mv[igrade][ibase] = ZERO
self.add_in_place(MV.btable[igrade][ibase]*coef)
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
grades of X.
"""
if isinstance(r, int):
grade_r = MV()
if r > MV.n:
return(grade_r)
self.convert_to_blades()
if not isinstance(self.mv[r], numpy.ndarray):
return(grade_r)
grade_r.bladeflg = 1
grade_r.puregrade = 1
grade_r.mv[r] = +self.mv[r]
return(grade_r)
if isinstance(r, MV):
self.convert_to_blades()
r.convert_to_blades()
proj = MV()
for i in MV.n1rg:
if not isinstance(r.mv[i], int):
proj.mv[i] = self.mv[i]
proj.bladeflg = self.bladeflg
return(proj)
return(None)
def even(self):
"""
Even grade projection operator. For multivector X, X.even()
returns multivector of even grade components of X.
"""
egrades = MV()
self.convert_to_blades()
egrades.bladeflg = self.bladeflg
egrades.puregrade = self.puregrade
for igrade in range(0, MV.n1, 2):
egrades.mv[igrade] = +self.mv[igrade]
return(egrades)
def odd(self):
"""
Odd grade projection operator. For multivector X, X.odd()
returns multivector of odd grade components of X.
"""
ogrades = MV()
self.convert_to_blades()
ogrades.bladeflg = self.bladeflg
ogrades.puregrade = self.puregrade
for igrade in range(1, MV.n1, 2):
ogrades.mv[igrade] = +self.mv[igrade]
return(ogrades)
def rev(self):
"""
Revision operator. For multivector X, X.rev()
returns reversed multivector of X.
"""
revmv = MV()
self.convert_to_blades()
revmv.bladeflg = self.bladeflg
revmv.puregrade = self.puregrade
for igrade in MV.n1rg:
if isinstance(self.mv[igrade], numpy.ndarray):
if igrade < 2 or (not (((igrade*(igrade - 1))/2) % 2)):
revmv.mv[igrade] = +self.mv[igrade]
else:
revmv.mv[igrade] = -self.mv[igrade]
return(revmv)
def cse(self, grade):
cse_lst = []
if isinstance(self.mv[grade], numpy.ndarray):
for ibase in range(MV.nbasis[grade]):
if self.mv[grade][ibase] != ZERO:
cse_lst.append(sympy.cse(self.mv[grade][ibase]))
return(cse_lst)
def div(self, grade, divisor):
div_lst = []
if isinstance(self.mv[grade], numpy.ndarray):
for ibase in range(MV.nbasis[grade]):
if self.mv[grade][ibase] != ZERO:
div_lst.append(
self.mv[grade][ibase].as_coefficient(divisor))
return(div_lst)
# don't know which one is correctly named
def div(self):
return self.grad_int()
def collect(self, lst):
"""
Applies sympy collect function to each component
of multivector.
"""
for igrade in MV.n1rg:
if isinstance(self.mv[igrade], numpy.ndarray):
for ibase in range(MV.nbasis[igrade]):
if self.mv[igrade][ibase] != ZERO:
self.mv[igrade][ibase] = \
sympy.collect(self.mv[igrade][ibase], lst)
return
def sqrfree(self, lst):
"""
Applies sympy sqrfree function to each component
of multivector.
"""
for igrade in MV.n1rg:
if isinstance(self.mv[igrade], numpy.ndarray):
for ibase in range(MV.nbasis[igrade]):
if self.mv[igrade][ibase] != ZERO:
self.mv[igrade][ibase] = \
sympy.sqrfree(self.mv[igrade][ibase], lst)
return
def flatten(self):
flst = []
for igrade in MV.n1rg:
if isinstance(self.mv[igrade], int):
flst += MV.nbasis[igrade]*[ZERO]
else:
for coef in self.mv[igrade]:
flst.append(coef)
return(flst)
def subs(self, *args):
X = MV()
X.bladeflg = self.bladeflg
X.puregrade = self.puregrade
for igrade in MV.n1rg:
if isinstance(self.mv[igrade], numpy.ndarray):
X.mv[igrade] = numpy.array(
substitute_array(self.mv[igrade], *args))
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:
for igrade in MV.n1rg:
if not isinstance(self.mv[igrade], int):
for ibase in range(MV.nbasis[igrade]):
if expr1 != ZERO:
self.mv[igrade][ibase] = self.mv[
igrade][ibase].subs(expr1, expr2)
return
def simplify(self):
"""
Applies sympy simplify function
to each component of multivector.
"""
for igrade in MV.n1rg:
if isinstance(self.mv[igrade], numpy.ndarray):
for ibase in range(MV.nbasis[igrade]):
if self.mv[igrade][ibase] != ZERO:
self.mv[igrade][ibase] = \
sympy.simplify(self.mv[igrade][ibase])
return
def trigsimp(self):
"""
Applies sympy trigsimp function
to each component of multivector
using the recursive option.
"""
for igrade in MV.n1rg:
if isinstance(self.mv[igrade], numpy.ndarray):
for ibase in range(MV.nbasis[igrade]):
if self.mv[igrade][ibase] != ZERO:
self.mv[igrade][ibase] = \
sympy.trigsimp(self.mv[igrade]
[ibase], deep=True, recursive=True)
return
def cancel(self):
"""
Applies sympy cancel function
to each component of multivector.
"""
for igrade in MV.n1rg:
if isinstance(self.mv[igrade], numpy.ndarray):
for ibase in range(MV.nbasis[igrade]):
if self.mv[igrade][ibase] != ZERO:
self.mv[igrade][ibase] = \
sympy.cancel(self.mv[igrade][ibase])
return
def expand(self):
"""
Applies sympy expand function to each component
of multivector.
"""
for igrade in MV.n1rg:
if isinstance(self.mv[igrade], numpy.ndarray):
for ibase in range(MV.nbasis[igrade]):
if self.mv[igrade][ibase] != ZERO:
self.mv[igrade][ibase] = \
self.mv[igrade][ibase].expand()
return
def is_pure(self):
igrade = -1
ngrade = 0
self.compact()
for i in MV.n1rg:
if isinstance(self.mv[i], numpy.ndarray):
for base in self.mv[i]:
if base != 0:
igrade = i
ngrade += 1
break
if ngrade > 1:
return(-1)
if igrade == -1:
return(0)
return(igrade)
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()
igrade = 0
for grade in self.mv:
if not isinstance(grade, int):
D.mv[igrade] = MV.vdiff(grade, x)
igrade += 1
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)
def grad(self):
"""
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):
dD.add_in_place(rbase*iD)
if MV.curvilinear_flg: # Add Connection
igrade = 1
while igrade <= MV.n:
coefs = self.mv[igrade]
if type(coefs) != int:
for (coef, connect) in zip(coefs, MV.Connect[igrade]):
dD.add_in_place(coef*connect)
igrade += 1
return(dD)
def grad_ext(self):
"""
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):
dD.add_in_place(irbase ^ iD)
if MV.curvilinear_flg: # Add Connection
igrade = 1
while igrade <= MV.n:
coefs = self.mv[igrade]
if type(coefs) != int:
for (coef, connect) in zip(coefs, MV.Connect[igrade]):
if igrade < MV.n:
dD.add_in_place(coef*connect.project(igrade + 1))
igrade += 1
return(dD)
def curl(self):
return(self.grad_ext())
def grad_int(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):
dD.add_in_place(irbase | iD)
if MV.curvilinear_flg: # Add Connection
igrade = 1
while igrade <= MV.n:
coefs = self.mv[igrade]
if type(coefs) != int:
for (coef, connect) in zip(coefs, MV.Connect[igrade]):
dD.add_in_place(coef*connect.project(igrade - 1))
igrade += 1
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'))