Vector math

[Maarten van Reeuwijk]

> 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