# Vector math

Maarten van Reeuwijk maarten at remove_this_ws.tn.tudelft.nl
Fri Mar 26 11:07:00 CET 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.

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

```