share/vector/vector.mac

   1
   2 /* *****  This package may be broken.  Please use the VECT
   3 package on the SHARE directory. - JPG 6/5/78  ***** */
   4
   5 vector:list$
   6 for f in ["grad","div","curl","laplacian","curlgrad","graddiv",
   7     "divcurl","curlcurl"] do prefix(f,112,expr,expr)$
   8 infix("cross",112,112,expr,expr,expr)$
   9 infix("dotdel",108,108,expr,expr,expr)$
  10 nofix("christoffel",expr)$
  11
  12 dimension():=if dimension=3 then [1,2,3] else [1,2]$
  13
  14 type(arg):=if list=  /* check required form of answer */
  15     (if listp(arg) then  /* operation performed on a vector */
  16         (for element in arg do  /* check each argument */
  17             if listp(element)  /* an argument is a list */
  18                 then return(list))  /* return a list */
  19         else vector)  /* operation performed on a scalar */
  20     then result  /* return a list */
  21     else apply('matrix,[result])  /* return a matrix */$
  22
  23 _coordsystem&&  coordsystem(sys):=
  24     (if sys=rectangular then
  25         (coodvar:[x,y],
  26         scalefactor:[1,1])
  27     else if sys=polar then
  28         (coordvar:[r,th],
  29         scalefactor:[1,r])
  30     else if sys=cartesian then
  31         (coordvar:[x,y,z],
  32         scalefactor:[1,1,1])
  33     else if sys=cylindrical then
  34         (coordvar:[r,ph,z],
  35         scalefactor:[1,r,1])
  36     else if sys=spherical then
  37         (coordvar:[r,th,ph],
  38         scalefactor:[1,r,r*sin(th)])
  39     else (coordvar:read("coordinate variables"),
  40         scalefactor:read("scale factors")),
  41     dimension:length(coordvar),
  42     coordsystem:sys)$
  43
  44 coordsystem(cartesian)$
  45
  46 _cross&&  (a cross b) := if dimension=3 then block([result],
  47         result:[a[2]*b[3]-a[3]*b[2],
  48             a[3]*b[1]-a[1]*b[3],
  49             a[1]*b[2]-a[2]*b[1]],
  50         type([a,b]))
  51     else  /* 2 dimensional case */
  52         if nonscalarp(a) then
  53         (if nonscalarp(b) then  /* vector x vector */
  54             a[1]*b[2]-a[2]*b[1]
  55         else block([result],  /* vector x scalar */
  56             result:[a[2]*b,
  57                 -a[1]*b],
  58             type([a])))
  59     else block([result],  /* scalar x vector */
  60         result:[-a*b[2],
  61             a*b[1]],
  62         type([b]))$
  63
  64 _grad&&  (grad s) := block([result],
  65     result:map(lambda([i],
  66         diff(s,coordvar[i])/scalefactor[i]),
  67         dimension()),
  68     type(vector))$
  69
  70 _div&&  (div v) := if dimension=3 then
  71         (diff(scalefactor[2]*scalefactor[3]*v[1],coordvar[1])+
  72         diff(scalefactor[3]*scalefactor[1]*v[2],coordvar[2])+
  73         diff(scalefactor[1]*scalefactor[2]*v[3],coordvar[3]))
  74         /scalefactor[1]/scalefactor[2]/scalefactor[3]
  75     else  /* 2 dimensional case */
  76         (diff(scalefactor[2]*v[1],coordvar[1])
  77         +diff(scalefactor[1]*v[2],coordvar[2]))
  78         /scalefactor[1]/scalefactor[2]$
  79
  80 _curl&&  (curl a) := if dimension=3 then block([result],
  81         result:[(diff(scalefactor[3]*a[3],coordvar[2])
  82                 -diff(scalefactor[2]*a[2],coordvar[3]))
  83                 /scalefactor[2]/scalefactor[3],
  84             (diff(scalefactor[1]*a[1],coordvar[3])
  85                 -diff(scalefactor[3]*a[3],coordvar[1]))
  86                 /scalefactor[3]/scalefactor[1],
  87             (diff(scalefactor[2]*a[2],coordvar[1])
  88                 -diff(scalefactor[1]*a[1],coordvar[2]))
  89                 /scalefactor[1]/scalefactor[2]],
  90         type([a]))
  91     else  /* 2 dimensional case */
  92         if nonscalarp(a) then block([result],
  93             result:(diff(scalefactor[2]*a[2],coordvar[1])
  94                 -diff(scalefactor[1]*a[1],coordvar[2]))
  95                 /scalefactor[1]/scalefactor[2],
  96             type([a]))
  97         else block([result],  /* scalar argument */
  98             result:[diff(a,coordvar[2])/scalefactor[2],
  99                 -diff(a,coordvar[1])/scalefactor[1]],
 100             type(vector))$
 101
 102 _laplacian&&  (laplacian a) := if nonscalarp(a) then grad div a -curl curl a
 103     else if dimension=3 then
 104         (diff(diff(a,coordvar[1])*scalefactor[2]
 105         *scalefactor[3]/scalefactor[1],coordvar[1])
 106         +diff(diff(a,coordvar[2])*scalefactor[3]
 107         *scalefactor[1]/scalefactor[2],coordvar[2])
 108         +diff(diff(a,coordvar[3])*scalefactor[1]
 109         *scalefactor[2]/scalefactor[3],coordvar[3]))
 110         /scalefactor[1]/scalefactor[2]/scalefactor[3]
 111     else  /* 2 dimensional case */
 112         (diff(diff(a,coordvar[1])*scalefactor[2]
 113         /scalefactor[1],coordvar[1])
 114         +diff(diff(a,coordvar[2])*scalefactor[1]
 115         /scalefactor[2],coordvar[2]))/scalefactor[1]/scalefactor[2]$
 116
 117 _dotdel&&  (v dotdel b) := if  nonscalarp(b) then block([result, b:flatten(args(b)), v:flatten(args(v))],
 118         result:if member('christsym, arrays)
 119             then  /* use christoffel symbols */
 120             map(lambda([j],
 121                 sum((diff(b[i]*scalefactor[j],
 122                     coordvar[i])-sum(b[k]*scalefactor[k]
 123                     *christsym[k,j,i],k,1,dimension))
 124                     *v[i]/scalefactor[i],i,1,dimension)
 125                     /scalefactor[j]),dimension())
 126         else  /* vector b, no christoffel symbols */
 127             map(lambda([j],
 128                 sum(diff(b[i]*scalefactor[j],coordvar[i])
 129                     *v[i]/scalefactor[i],i,1,dimension)
 130                     /scalefactor[j]),dimension()),
 131         type([v,b]))
 132     else block([result],  /* scalar b case */
 133         result:if member('christsym, arrays)
 134             then /* use christoffel symbols */
 135             sum((diff(b,coordvar[i])-b*christsym[1,1,i])
 136                 *v[i]/scalefactor[i],i,1,dimension)
 137         else /* no christoffel symbols */
 138             sum(diff(b,coordvar[i])*v[i]
 139                 /scalefactor[i],i,1,dimension),
 140         type([v]))$
 141
 142 _christoffel&&  christoffel := (array(christsym,3,3,3),
 143     christsym[i,j,k]:=0,
 144     for i thru 3 do
 145         (christsym[i,i,i]:diff(scalefactor[i],coordvar[i])
 146             /scalefactor[i],
 147         for j thru 3 do if j#i then
 148             (christsym[i,j,i]:christsym[i,i,j]:diff(scalefactor[i],
 149                 coordvar[j])/scalefactor[i],
 150             christsym[j,i,i]:-diff(scalefactor[i],
 151                 coordvar[j])*scalefactor[i]/scalefactor[j]^2)))$
 152
 153 (curlgrad s) := 0$
 154
 155 (graddiv v) := block([result],
 156     result:div v,
 157     result:map(lambda([i],
 158         diff(result,coordvar[i])/scalefactor[i]),
 159         dimension()),
 160     type(vector))$
 161
 162 (divcurl v) := 0$