from sympy.vector.coordsysrect import CoordSysCartesian
from sympy.vector.dyadic import Dyadic
from sympy.vector.vector import Vector, BaseVector
from sympy.vector.scalar import BaseScalar
from sympy import sympify, diff, integrate, S, simplify
[docs]def express(expr, system, system2=None, variables=False):
"""
Global function for 'express' functionality.
Re-expresses a Vector, Dyadic or scalar(sympyfiable) in the given
coordinate system.
If 'variables' is True, then the coordinate variables (base scalars)
of other coordinate systems present in the vector/scalar field or
dyadic are also substituted in terms of the base scalars of the
given system.
Parameters
==========
expr : Vector/Dyadic/scalar(sympyfiable)
The expression to re-express in CoordSysCartesian 'system'
system: CoordSysCartesian
The coordinate system the expr is to be expressed in
system2: CoordSysCartesian
The other coordinate system required for re-expression
(only for a Dyadic Expr)
variables : boolean
Specifies whether to substitute the coordinate variables present
in expr, in terms of those of parameter system
Examples
========
>>> from sympy.vector import CoordSysCartesian
>>> from sympy import Symbol, cos, sin
>>> N = CoordSysCartesian('N')
>>> q = Symbol('q')
>>> B = N.orient_new_axis('B', q, N.k)
>>> from sympy.vector import express
>>> express(B.i, N)
(cos(q))*N.i + (sin(q))*N.j
>>> express(N.x, B, variables=True)
-sin(q)*B.y + cos(q)*B.x
>>> d = N.i.outer(N.i)
>>> express(d, B, N) == (cos(q))*(B.i|N.i) + (-sin(q))*(B.j|N.i)
True
"""
if expr == 0 or expr == Vector.zero:
return expr
if not isinstance(system, CoordSysCartesian):
raise TypeError("system should be a CoordSysCartesian \
instance")
if isinstance(expr, Vector):
if system2 is not None:
raise ValueError("system2 should not be provided for \
Vectors")
#Given expr is a Vector
if variables:
#If variables attribute is True, substitute
#the coordinate variables in the Vector
system_list = []
for x in expr.atoms(BaseScalar, BaseVector):
if x.system != system:
system_list.append(x.system)
system_list = set(system_list)
subs_dict = {}
for f in system_list:
subs_dict.update(f.scalar_map(system))
expr = expr.subs(subs_dict)
#Re-express in this coordinate system
outvec = Vector.zero
parts = expr.separate()
for x in parts:
if x != system:
temp = system.rotation_matrix(x) * parts[x].to_matrix(x)
outvec += matrix_to_vector(temp, system)
else:
outvec += parts[x]
return outvec
elif isinstance(expr, Dyadic):
if system2 is None:
system2 = system
if not isinstance(system2, CoordSysCartesian):
raise TypeError("system2 should be a CoordSysCartesian \
instance")
outdyad = Dyadic.zero
var = variables
for k, v in expr.components.items():
outdyad += (express(v, system, variables=var) *
(express(k.args[0], system, variables=var) |
express(k.args[1], system2, variables=var)))
return outdyad
else:
if system2 is not None:
raise ValueError("system2 should not be provided for \
Vectors")
if variables:
#Given expr is a scalar field
system_set = set([])
expr = sympify(expr)
#Subsitute all the coordinate variables
for x in expr.atoms(BaseScalar):
if x.system != system:
system_set.add(x.system)
subs_dict = {}
for f in system_set:
subs_dict.update(f.scalar_map(system))
return expr.subs(subs_dict)
return expr
[docs]def curl(vect, coord_sys):
"""
Returns the curl of a vector field computed wrt the base scalars
of the given coordinate system.
Parameters
==========
vect : Vector
The vector operand
coord_sys : CoordSysCartesian
The coordinate system to calculate the curl in
Examples
========
>>> from sympy.vector import CoordSysCartesian, curl
>>> R = CoordSysCartesian('R')
>>> v1 = R.y*R.z*R.i + R.x*R.z*R.j + R.x*R.y*R.k
>>> curl(v1, R)
0
>>> v2 = R.x*R.y*R.z*R.i
>>> curl(v2, R)
R.x*R.y*R.j + (-R.x*R.z)*R.k
"""
return coord_sys.delop.cross(vect).doit()
[docs]def divergence(vect, coord_sys):
"""
Returns the divergence of a vector field computed wrt the base
scalars of the given coordinate system.
Parameters
==========
vect : Vector
The vector operand
coord_sys : CoordSysCartesian
The cooordinate system to calculate the divergence in
Examples
========
>>> from sympy.vector import CoordSysCartesian, divergence
>>> R = CoordSysCartesian('R')
>>> v1 = R.x*R.y*R.z * (R.i+R.j+R.k)
>>> divergence(v1, R)
R.x*R.y + R.x*R.z + R.y*R.z
>>> v2 = 2*R.y*R.z*R.j
>>> divergence(v2, R)
2*R.z
"""
return coord_sys.delop.dot(vect).doit()
[docs]def gradient(scalar, coord_sys):
"""
Returns the vector gradient of a scalar field computed wrt the
base scalars of the given coordinate system.
Parameters
==========
scalar : SymPy Expr
The scalar field to compute the gradient of
coord_sys : CoordSysCartesian
The coordinate system to calculate the gradient in
Examples
========
>>> from sympy.vector import CoordSysCartesian, gradient
>>> R = CoordSysCartesian('R')
>>> s1 = R.x*R.y*R.z
>>> gradient(s1, R)
R.y*R.z*R.i + R.x*R.z*R.j + R.x*R.y*R.k
>>> s2 = 5*R.x**2*R.z
>>> gradient(s2, R)
10*R.x*R.z*R.i + 5*R.x**2*R.k
"""
return coord_sys.delop(scalar).doit()
[docs]def is_conservative(field):
"""
Checks if a field is conservative.
Paramaters
==========
field : Vector
The field to check for conservative property
Examples
========
>>> from sympy.vector import CoordSysCartesian
>>> from sympy.vector import is_conservative
>>> R = CoordSysCartesian('R')
>>> is_conservative(R.y*R.z*R.i + R.x*R.z*R.j + R.x*R.y*R.k)
True
>>> is_conservative(R.z*R.j)
False
"""
#Field is conservative irrespective of system
#Take the first coordinate system in the result of the
#separate method of Vector
if not isinstance(field, Vector):
raise TypeError("field should be a Vector")
if field == Vector.zero:
return True
coord_sys = list(field.separate())[0]
return curl(field, coord_sys).simplify() == Vector.zero
[docs]def is_solenoidal(field):
"""
Checks if a field is solenoidal.
Paramaters
==========
field : Vector
The field to check for solenoidal property
Examples
========
>>> from sympy.vector import CoordSysCartesian
>>> from sympy.vector import is_solenoidal
>>> R = CoordSysCartesian('R')
>>> is_solenoidal(R.y*R.z*R.i + R.x*R.z*R.j + R.x*R.y*R.k)
True
>>> is_solenoidal(R.y * R.j)
False
"""
#Field is solenoidal irrespective of system
#Take the first coordinate system in the result of the
#separate method in Vector
if not isinstance(field, Vector):
raise TypeError("field should be a Vector")
if field == Vector.zero:
return True
coord_sys = list(field.separate())[0]
return divergence(field, coord_sys).simplify() == S(0)
[docs]def scalar_potential(field, coord_sys):
"""
Returns the scalar potential function of a field in a given
coordinate system (without the added integration constant).
Parameters
==========
field : Vector
The vector field whose scalar potential function is to be
calculated
coord_sys : CoordSysCartesian
The coordinate system to do the calculation in
Examples
========
>>> from sympy.vector import CoordSysCartesian
>>> from sympy.vector import scalar_potential, gradient
>>> R = CoordSysCartesian('R')
>>> scalar_potential(R.k, R) == R.z
True
>>> scalar_field = 2*R.x**2*R.y*R.z
>>> grad_field = gradient(scalar_field, R)
>>> scalar_potential(grad_field, R)
2*R.x**2*R.y*R.z
"""
#Check whether field is conservative
if not is_conservative(field):
raise ValueError("Field is not conservative")
if field == Vector.zero:
return S(0)
#Express the field exntirely in coord_sys
#Subsitute coordinate variables also
if not isinstance(coord_sys, CoordSysCartesian):
raise TypeError("coord_sys must be a CoordSysCartesian")
field = express(field, coord_sys, variables=True)
dimensions = coord_sys.base_vectors()
scalars = coord_sys.base_scalars()
#Calculate scalar potential function
temp_function = integrate(field.dot(dimensions[0]), scalars[0])
for i, dim in enumerate(dimensions[1:]):
partial_diff = diff(temp_function, scalars[i + 1])
partial_diff = field.dot(dim) - partial_diff
temp_function += integrate(partial_diff, scalars[i + 1])
return temp_function
[docs]def scalar_potential_difference(field, coord_sys, point1, point2):
"""
Returns the scalar potential difference between two points in a
certain coordinate system, wrt a given field.
If a scalar field is provided, its values at the two points are
considered. If a conservative vector field is provided, the values
of its scalar potential function at the two points are used.
Returns (potential at point2) - (potential at point1)
The position vectors of the two Points are calculated wrt the
origin of the coordinate system provided.
Parameters
==========
field : Vector/Expr
The field to calculate wrt
coord_sys : CoordSysCartesian
The coordinate system to do the calculations in
point1 : Point
The initial Point in given coordinate system
position2 : Point
The second Point in the given coordinate system
Examples
========
>>> from sympy.vector import CoordSysCartesian, Point
>>> from sympy.vector import scalar_potential_difference
>>> R = CoordSysCartesian('R')
>>> P = R.origin.locate_new('P', R.x*R.i + R.y*R.j + R.z*R.k)
>>> vectfield = 4*R.x*R.y*R.i + 2*R.x**2*R.j
>>> scalar_potential_difference(vectfield, R, R.origin, P)
2*R.x**2*R.y
>>> Q = R.origin.locate_new('O', 3*R.i + R.j + 2*R.k)
>>> scalar_potential_difference(vectfield, R, P, Q)
-2*R.x**2*R.y + 18
"""
if not isinstance(coord_sys, CoordSysCartesian):
raise TypeError("coord_sys must be a CoordSysCartesian")
if isinstance(field, Vector):
#Get the scalar potential function
scalar_fn = scalar_potential(field, coord_sys)
else:
#Field is a scalar
scalar_fn = field
#Express positions in required coordinate system
origin = coord_sys.origin
position1 = express(point1.position_wrt(origin), coord_sys,
variables=True)
position2 = express(point2.position_wrt(origin), coord_sys,
variables=True)
#Get the two positions as substitution dicts for coordinate variables
subs_dict1 = {}
subs_dict2 = {}
scalars = coord_sys.base_scalars()
for i, x in enumerate(coord_sys.base_vectors()):
subs_dict1[scalars[i]] = x.dot(position1)
subs_dict2[scalars[i]] = x.dot(position2)
return scalar_fn.subs(subs_dict2) - scalar_fn.subs(subs_dict1)
[docs]def matrix_to_vector(matrix, system):
"""
Converts a vector in matrix form to a Vector instance.
It is assumed that the elements of the Matrix represent the
measure numbers of the components of the vector along basis
vectors of 'system'.
Parameters
==========
matrix : SymPy Matrix, Dimensions: (3, 1)
The matrix to be converted to a vector
system : CoordSysCartesian
The coordinate system the vector is to be defined in
Examples
========
>>> from sympy import ImmutableMatrix as Matrix
>>> m = Matrix([1, 2, 3])
>>> from sympy.vector import CoordSysCartesian, matrix_to_vector
>>> C = CoordSysCartesian('C')
>>> v = matrix_to_vector(m, C)
>>> v
C.i + 2*C.j + 3*C.k
>>> v.to_matrix(C) == m
True
"""
outvec = Vector.zero
vects = system.base_vectors()
for i, x in enumerate(matrix):
outvec += x * vects[i]
return outvec
def _path(from_object, to_object):
"""
Calculates the 'path' of objects starting from 'from_object'
to 'to_object', along with the index of the first common
ancestor in the tree.
Returns (index, list) tuple.
"""
if from_object._root != to_object._root:
raise ValueError("No connecting path found between " +
str(from_object) + " and " + str(to_object))
other_path = []
obj = to_object
while obj._parent is not None:
other_path.append(obj)
obj = obj._parent
other_path.append(obj)
object_set = set(other_path)
from_path = []
obj = from_object
while obj not in object_set:
from_path.append(obj)
obj = obj._parent
index = len(from_path)
i = other_path.index(obj)
while i >= 0:
from_path.append(other_path[i])
i -= 1
return index, from_path
def orthogonalize(*vlist, **kwargs):
"""
Takes a sequence of independent vectors and orthogonalizes them
using the Gram - Schmidt process. Returns a list of
orthogonal or orthonormal vectors.
Parameters
==========
vlist : sequence of independent vectors to be made orthogonal.
orthonormal : Optional parameter
Set to True if the the vectors returned should be orthonormal.
Default: False
Examples
========
>>> from sympy.vector.coordsysrect import CoordSysCartesian
>>> from sympy.vector.vector import Vector, BaseVector
>>> from sympy.vector.functions import orthogonalize
>>> C = CoordSysCartesian('C')
>>> i, j, k = C.base_vectors()
>>> v1 = i + 2*j
>>> v2 = 2*i + 3*j
>>> orthogonalize(v1, v2)
[C.i + 2*C.j, 2/5*C.i + (-1/5)*C.j]
References
==========
.. [1] https://en.wikipedia.org/wiki/Gram-Schmidt_process
"""
orthonormal = kwargs.get('orthonormal', False)
if not all(isinstance(vec, Vector) for vec in vlist):
raise TypeError('Each element must be of Type Vector')
ortho_vlist = []
for i, term in enumerate(vlist):
for j in range(i):
term -= ortho_vlist[j].projection(vlist[i])
# TODO : The following line introduces a performance issue
# and needs to be changed once a good solution for issue #10279 is found.
if simplify(term).equals(Vector.zero):
raise ValueError("Vector set not linearly independent")
ortho_vlist.append(term)
if orthonormal:
ortho_vlist = [vec.normalize() for vec in ortho_vlist]
return ortho_vlist
