share/contrib/vector3d/vector3d.mac

   1 /* File name: vector3d.mac
   2  * Package providing vector algebra and differential calculus operations
   3  * on three dimensional vectors in orthogonal coordinate systems.
   4  * Implements cross product of two vectors, div, curl, and laplacian
   5  * of vector functions of the coordinates. and grad and laplacian
   6  * of scalar functions of the coordinates.  The vector arguments must
   7  * be three element lists.
   8  * Walter Eastes: June, 2004
   9  * Copyright updated: July, 2007
  10  */
  11
  12 /*
  13 Copyright (C) 2004-2007 by Walter Eastes
  14
  15 These routines are free software; you can redistribute them and/or modify
  16 them under the terms of the GNU General Public License as published by
  17 the Free Software Foundation; either version 2 of the License, or
  18 (at your option) any later version.
  19 */
  20
  21 coordsys(sys):=
  22 (  if sys=cartesian then (scalef : [1,1,1], coordvar : [x,y,z])
  23    else if sys=cylindrical then (scalef : [1,r,1], coordvar : [r,ph,z])
  24    else if sys=spherical then (scalef : [1,r*sin(ph),r], coordvar : [r,th,ph])
  25    else (scalef : read("Scale factors"), coordvar : read("Coordinates"))
  26 )$
  27
  28 coordsys(cartesian)$
  29
  30 cross(a,b) :=
  31    [ a[2]*b[3]-a[3]*b[2], a[3]*b[1]-a[1]*b[3], a[1]*b[2]-a[2]*b[1] ] $
  32
  33 div(v) :=
  34 (  diff(scalef[2]*scalef[3]*v[1], coordvar[1]) +
  35    diff(scalef[3]*scalef[1]*v[2], coordvar[2]) +
  36    diff(scalef[1]*scalef[2]*v[3], coordvar[3])
  37 ) / (scalef[1]*scalef[2]*scalef[3]) $
  38
  39 curl(a) :=
  40 [  (diff(scalef[3]*a[3], coordvar[2]) -
  41     diff(scalef[2]*a[2], coordvar[3])) / (scalef[2]*scalef[3]),
  42
  43    (diff(scalef[1]*a[1], coordvar[3]) -
  44     diff(scalef[3]*a[3], coordvar[1])) / (scalef[3]*scalef[1]),
  45
  46    (diff(scalef[2]*a[2], coordvar[1]) -
  47     diff(scalef[1]*a[1], coordvar[2])) / (scalef[1]*scalef[2])
  48 ]$
  49
  50 grad(f) :=
  51 [  diff(f,coordvar[1]) / scalef[1],
  52    diff(f,coordvar[2]) / scalef[2],
  53    diff(f,coordvar[3]) / scalef[3]
  54 ]$
  55
  56 laplacian(v) :=
  57    if listp(v) then grad(div(v)) - curl(curl(v))
  58    else
  59    (  diff(diff(v,coordvar[1]) * scalef[2]*scalef[3]/scalef[1], coordvar[1]) +
  60       diff(diff(v,coordvar[2]) * scalef[3]*scalef[1]/scalef[2], coordvar[2]) +
  61       diff(diff(v,coordvar[3]) * scalef[1]*scalef[2]/scalef[3], coordvar[3])
  62    ) / (scalef[1]*scalef[2]*scalef[3])
  63 $
  64