Vector math
Maarten van Reeuwijk
maarten at remove_this_ws.tn.tudelft.nl
Fri Mar 26 05:07:00 EST 2004
Josiah Carlson wrote:
> Check numarray:
> http://www.stsci.edu/resources/software_hardware/numarray
> and Numeric Python:
> http://www.pfdubois.com/numpy/
I bumped into the same problem, but couldn't find those features, so I wrote
them myself using the Numeric package. Below follows an example where I
rotate the z-axis 45 degrees.
>>> from manips import *
>>> from Numeric import *
>>> a = array([1,1,1.])
>>> base = [array([1,0,0]), array([0,1,0]), array([0,0,1])]
>>> newbase = rotatebase(base, [0,0,45 / 180. * 3.14])
>>> changebase(a, newbase)
array([ 1.41421345e+00, 5.63088062e-04, 1.00000000e+00])
>>> a = transpose(array([ [1,1,1.],[-1,-1,-1.]]))
>>> transpose(changebase(a, newbase)
... )
array([[ 1.41421345e+00, 5.63088062e-04, 1.00000000e+00],
[ -1.41421345e+00, -5.63088062e-04, -1.00000000e+00]])
the 'transpose' commands in the second example are needed because the
changebase algorhytm expects an array containing only x, y and z values.
Below follow the two routines.
HTH, Maarten
def changebase(vecfield, newbase):
""" Changes the base of the vectorfield. newbase is a list containing
the
new basevectors (each is a 1D-array). vecfield is a list containing
the fields for the individual vector components fields. The lenght of
the vecfield and newbase lists must be identical. The new field is
calculated as:
*
u = u b
i j i,j
Clearly, the basevectors are the new principal directions relative
to the CURRENT coordinate system.
"""
vf = vecfield
nb = newbase
nf = []
if size([vf[0]]) == 1:
isfield = False
tc = vf.typecode()
else:
isfield = True
tc = vf[0].typecode()
tmp = array(vf[0], tc, savespace = 1)
for bcomp in range(len(newbase)):
if isfield:
tmp[:] = 0
else:
tmp = 0
for vcomp in range(len(newbase)):
tmp = tmp + nb[bcomp][vcomp] * vf[vcomp]
nf.append(array(tmp))
if not isfield:
nf = array(nf, tc, savespace = 1).flat
return nf
def rotatebase(base, rotation):
""" Rotates the base (given as a list of lists or array objects) over a
angle given by rotation. The transformation happens according to:
_ _ _ _
| 1 0 0 | | cos b 0 sin b |
Rx = | 0 cos a -sin a | ; Ry = | 0 1 0 |
|_ 0 sin a cos a _| |_-sin b 0 cos b _|
_ _
| cos c -sin c 0 |
Rz = | sin c cos c 0 |
|_ 0 0 1_|
z
^
|
|
| Rx : rotation along x-axis (a pos from y
-> z)
|--------> y Ry : rotation along y-axis (b pos from z
-> x)
/ Rz : rotation along z-axis (c pos from x
-> y)
/
x
newbase = Rx * Ry * Rz * base (all matrix multiplications)
"""
alpha, beta, gamma = rotation
sa = sin(alpha); ca = cos(alpha)
sb = sin(beta); cb = cos(beta)
sg = sin(gamma); cg = cos(gamma)
rot = array([[cb * cg , -cb * sg ,
sb ], \
[sa * sb * cg + ca * sg, -sa * sb * sg + ca * cg, -sa *
cb], \
[-ca * sb * cg + sa * sg, ca * sb * sg + sa * cg, ca *
cb ]])
newbase = []
for bv in base:
newbase.append(matrixmultiply(rot, array(bv)))
return newbase
--
===================================================================
Maarten van Reeuwijk Thermal and Fluids Sciences
Phd student dept. of Multiscale Physics
www.ws.tn.tudelft.nl Delft University of Technology
More information about the Python-list
mailing list