module avector; % Vector algebra and calculus package.
create!-package('(avector),'(contrib avector));
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% %
% Copyright notice %
% ---------------- %
% %
% Author : David Harper %
% Computer Laboratory, %
% University of Liverpool %
% P.O. Box 147 %
% Liverpool L69 3BX %
% ENGLAND %
% %
% (adh@maths.qmw.ac.uk) %
% %
% Date : 29 February 1988 %
% %
% Title : Vector algebra and calculus package %
% %
% Copyright (c) David Harper 1988 %
% %
% %
% Permission is granted to any individual or institution to %
% use, copy or re-distribute this software so long as it is %
% not sold for profit, provided that this copyright notice %
% is retained and the file is not altered. %
% %
% %
% End of copyright notice %
% %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% *************** This code is designed to operate *************** %
% *************** with version 3.4 of REDUCE and *************** %
% *************** the Standard Lisp interpreter. *************** %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% %
% VECTOR DECLARATIONS MODULE %
% %
% This section contains the routines required to interface the %
% vector package to REDUCE. The most important routine is the %
% vector predicate function VECP which tests an expression to %
% determine whether it must be evaluated using the routine %
% VECSM*. %
% %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
fluid '(!*coords !*csystems !*vectorfunctionlist
!*vtrace !*vectortracelevel!*);
switch vtrace;
symbolic procedure vecp u;
begin scalar x;
return
if null u or numberp u then nil
else if atom u then (get(u,'rtype)='!3vector or threevectorp u)
else if threevectorp u then t
else if (atom(x:=car u) and get(x,'rtype)='!3vector)
then isvectorindex cadr u
else if flagp(x,'vectorfn) then t
else if (flagp(x,'varithop) or flagp(x,'vectormapping))
then hasonevector cdr u
else nil;
end;
% The following procedure checks for a vector with three components
symbolic procedure threevectorp u;
vectorp u and (upbv u = 2);
% The following procedure checks to ensure that the arg list contains
% at least one vector
symbolic procedure hasonevector u;
if null u then nil
else (vecp car u) or (hasonevector cdr u);
% The following procedure checks that the arg list consists entirely
% of vectors
symbolic procedure areallvectors u;
if null u then nil
else if null cdr u then vecp car u
else (vecp car u) and (areallvectors cdr u);
% The following function checks to see whether its argument is a valid
% vector index i.e. whether it evaluates to an integer 0,1 or 2 or
% is equal to one of the coordinate names
symbolic procedure isvectorindex u;
not null getvectorindex(u,t);
% The following function evaluates its argument to a vector index
% or NIL if it isn't an integer 0,1 or 2 or one of the coordinate
% names. Set FLG to true if hard-error is required for invalid
% argument.
symbolic procedure getvectorindex(u,flg);
begin scalar vindx;
vindx := u;
if not fixp vindx then vindx:=locate(vindx,!*coords);
if ((null vindx) or (fixp vindx and (vindx<0 or vindx>2)))
and flg then rerror(avector,1,list(u,"not a valid vector index"));
return vindx
end;
% This routine gives the position of an object in a list. The first
% object is numbered zero. Returns NIL if the object can't be found.
symbolic procedure locate(u,v);
if not (u memq v) then nil
else if u=car v then 0
else 1+locate(u,cdr v);
% We may as well define some utility operators here too.
symbolic smacro procedure first u; car u;
symbolic smacro procedure second u; cadr u;
symbolic smacro procedure third u; caddr u;
% Here we redefine getrtype1 and getrtype2 to handle vectors.
remflag('(getrtype1 getrtype2),'lose); % We must use these definitions.
symbolic procedure getrtype1 u;
if threevectorp u then '!3vector else nil;
symbolic procedure getrtype2 u;
begin scalar x;
% Next line is maybe only needed by EXCALC.
return if vecp u then '!3vector
else if (x:= get(car u,'rtype)) and (x:= get(x,'rtypefn))
then apply1(x,cdr u)
else if x := get(car u,'rtypefn) then apply1(x,cdr u)
else nil
end;
% The following function declares a list of objects as vectors.
symbolic procedure vec u;
begin scalar y;
for each x in u do
<< % Check that the identifier is not already declared as an array
% or matrix or function
if not atom x then write("Cannot declare ",x," as a vector")
else
<< y := gettype x;
if y memq '(array procedure matrix operator) then
write("Object ",x," has already been declared as ",y)
else makenewvector x
>>;
>>;
return nil
end;
deflist('((vec rlis)),'stat);
% This procedure actually creates a vector.
symbolic procedure makenewvector u;
begin
% write("Declaring ",U," a vector");terpri();
put(u,'rtype,'!3vector);
put(u,'avalue,list('vector,mkvect(2)));
return nil
end;
% Vector function declarations : these are the routines that link
% our new data type into the REDUCE system.
put('!3vector,'letfn,'veclet); % Assignment routine
put('!3vector,'name,'!3vector); % Our name for the data type
put('!3vector,'evfn,'!*vecsm!*); % Evaluation function
put('!3vector,'prifn,'vecpri!*); % Printing function
flag('(!3vector),'sprifn);
% Here we re-define the routine VARPRI to include 3-vectors. We
% add a single line to the routine.
symbolic procedure varpri(u,v,w);
begin scalar x;
%U is expression being printed
%V is the original form that was evaluated.
%W is an id that indicates if U is the first, only or last element
% in the current set (or NIL otherwise).
if null u then u := 0;
if !*nero and u=0 then return nil;
v := setvars v;
if vecp u then return vecpri(u,'mat);
if (x := getrtype u) and flagp(x,'sprifn)
then return if null v then apply1(get(get(x,'tag),'prifn),u)
else maprin list('setq,car v,u);
if w memq '(first only) then terpri!* t;
if !*fort then return fvarpri(u,v,w);
if v then u := 'setq . aconc(v,u);
maprin u;
if null w or w eq 'first then return nil
else if not !*nat then prin2!* "$";
terpri!*(not !*nat);
return nil
end;
symbolic procedure getavalue u;
(if x then cadr x else nil) where x=get(u,'avalue);
% The following routine prints out a vector in a neat way (cf.
% the way in which matrices are printed !)
symbolic procedure vecpri(u,x);
begin scalar y,v0,v1,v2,xx;
y:= if vectorp u then u else getavalue u;
xx := x;
if null y then return nil;
% if null x then xx := 'vec;
xx := 'vec;
v0 := getv(y,0);
v1 := getv(y,1);
v2 := getv(y,2);
v0 := aeval v0; v1 := aeval v1; v2 := aeval v2;
varpri(v0,list('setq, list(xx,aeval first !*coords),v0),'only);
varpri(v1,list('setq, list(xx,aeval second !*coords),v1),'only);
varpri(v2,list('setq, list(xx,aeval third !*coords),v2),'only);
terpri!* t;
end;
% This routine is the 'hook' for VARPRI
symbolic procedure vecpri!*(u,v,w);
vecpri(u,v);
symbolic procedure indexedvectorp u;
(vecp car u) and (isvectorindex cadr u);
put('!3vector,'setelemfn,'setvectorelement);
% The following function sets one element of a vector object
symbolic procedure setvectorelement(u,v);
begin scalar vindx;
vindx := getvectorindex(cadr u,t);
putv(getavalue car u,vindx,v);
return nil
end;
% If SETK is invoked with an vector atom as its first argument, then
% control will be passed to the routine VECLET
symbolic procedure veclet(u,v,utype,b,vtype);
begin
if zerop v then return setvectortozero(u,utype);
if not equal(vtype,'!3vector)
then rerror(avector,2,"RHS is not a vector");
if equal(utype,'!3vector) then put(u,'avalue,list('!3vector,v))
else if utype memq '(array matrix) then
rerror(avector,3,list(u,"already defined as ",utype))
else
<< % We force U to be a vector
vec u;
write("*** ",u," re-defined as vector"); terpri();
put(u,'avalue,list('vector,v))
>>;
return v
end;
% A quick and dirty way of declaring a null vector
symbolic procedure setvectortozero(u,utype);
begin scalar x;
x := mkvect 2;
for k:=0:2 do putv(x,k,aeval 0);
return veclet(u,x,utype,t,'!3vector)
end;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% %
% VECTOR EVALUATION MODULE %
% %
% This section contains the routines required to evaluate vector %
% expressions. The main routine, VECSM!*, is mainly table-driven. %
% If you wish to include your own routines then you should be %
% aware of the mechanism used to invoke vector evaluation. %
% %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
symbolic procedure !*vecsm!*(u,v); vecsm!* u;
% The following routine is the main vector evaluation procedure. It
% takes a single argument which is either a vector atom or a
% vector expression in SLISP function form. The value returned may
% be a vector or scalar.
% Note that if the argument is not a vector expression then an
% error results. This is true in particular if the argument is
% an expression in standard-quotient form i.e. a list whose CAR
% is *SQ.
!*vectortracelevel!* := 0;
symbolic procedure prtblanks n;
for k:=1:min(n,15) do write " ";
symbolic procedure vecsimp!* u;
if vecp u then vecsm!* u else u;
symbolic procedure vecsm!* u;
begin scalar y,vecop,vargs,v;
!*vectortracelevel!* := !*vectortracelevel!* + 1;
if !*vtrace then <<prtblanks(!*vectortracelevel!*);
write("VECSM called with args ",u);
terpri()>>;
if atom u then v := (if vectorp u then u else getavalue u)
else if threevectorp u then v := u
else if (atom(y:= car u) and get(y,'rtype)='!3vector) then
v := getv(getavalue y,getvectorindex(cadr u,t))
% Those were the simple cases. Now for the tricky operators
% Separate the operator from its operands
else <<
vecop := car u;
vargs := mapcar(cdr u,'vecsimp!*);
% Select a course of action dependent upon the operator
if y := get(vecop,'vectorfunction) then
% We must check (if the op is an arithmetic operator)
% to ensure that there are vectors in the argument list
<< if (flagp(vecop,'vectorfn) or (flagp(vecop,'varithop) and
hasonevector vargs)) then v := apply(y,list vargs)
else v := aeval append(list(vecop),vargs)
>>
else if flagp(vecop,'vectormapping) then
<< % Check that the argument is really a vector
y := car vargs;
v := if threevectorp y then vectorapply(vecop,y,cdr vargs)
else scalarapply(vecop,y,cdr vargs)
>>
else <<!*vectortracelevel!* := 0;
rerror(avector,4,
list(vecop,"is not a valid vector operator"))>>;
>>;
if !*vtrace then <<
y := threevectorp v;
prtblanks(!*vectortracelevel!*);
write(if y then "** Vector" else "** Scalar"," result is ",v);
terpri()>>;
!*vectortracelevel!* := !*vectortracelevel!* - 1;
return v
end;
% Now we define a function to declare a list of scalar functions as
% valid in vector mode too. This means that we can pass them vector
% arguments and get a vector result.
symbolic procedure vectormapping u;
flag(u,'vectormapping);
deflist('((vectormapping rlis)),'stat);% Deflist used for bootstrapping.
% We will allow the basic transcendental functions to be vector-valued.
% Then we can, for example, evaluate Sin of a vector.
vectormapping 'sin,'cos,'log,'exp,'tan,'asin,'atan,'sinh,'cosh,'tanh;
vectormapping 'quotient,'minus,'df,'int,'sqrt;
% We must put appropriate flags upon the arithmetic operators and
% vector operators too ...
flag('(sub minus difference quotient plus times expt),'varithop);
flag('(avec cross dot vmod grad div curl delsq),'vectorfn);
% We must now define the procedures to carry out vector algebra and
% calculus. They must be given a VECTORFUNCTION property
symbolic smacro procedure vectorfn(oper,vfn);
put(oper,'vectorfunction,vfn);
% Scalar-vector multiplication
vectorfn('times,'vectormultiply);
symbolic procedure vectormultiply vargs;
% This routine multiplies together a list made up of scalars,
% 3x3 matrices and a vector. Note that the combinations
% vector*vector and vector*matrix are illegal.
begin scalar lht,rht,lhtype,rhtype;
lht := aeval car vargs; % Begin with first multiplicand
for each v in cdr vargs do
<< % Deal with each multiplicand in turn
rht := if vecp v then vecsm!* v
else v;
lhtype := !*typeof lht;
rhtype := !*typeof rht;
lht :=
if not (lhtype='!3vector or rhtype='!3vector) then
aeval list('times,lht,rht)
else if lhtype='!3vector then
if null rhtype then vectorapply('times,lht,list rht)
else rerror(avector,5,"Illegal operation vec*vec or vec*mat")
else if null lhtype then vectorapply('times,rht,list lht)
else matrixtimesvector(lht,rht)
>>;
return lht
end;
% Multiplication of a vector by a 3x3 matrix from the left
symbolic procedure matrixtimesvector(mymat,myvec);
begin scalar rows,myrow,x;
if atom mymat and idp mymat and null getavalue mymat then
rerror(avector,6,"Unset matrix in vector multiplication");
rows := if idp mymat then cdr getavalue mymat else cdr mymat;
if not (length(rows)=3 and length(car rows)=3) then
rerror(avector,7,"Matrix must be 3x3 for vector multplication");
x := mkvect(2);
for k:=0:2 do
<< % Multiply out a row at a time
myrow := car rows;
putv(x,k,aeval list('plus,
list('times, first myrow, getv(myvec,0)),
list('times, second myrow, getv(myvec,1)),
list('times, third myrow, getv(myvec,2))));
rows := cdr rows
>>;
return x
end;
symbolic procedure !*typeof u;
getrtype u;
% if vecp u then '!3vector
% else if matp u then 'matrix
% else if arrayp u then 'array
% else nil;
% Vector addition
vectorfn('plus,'vectorplus);
symbolic procedure vectorplus vargs;
% Add an arbitrarily-long list of vectors
begin scalar x;
x := vecsm!* car vargs;
for each v in cdr vargs do x:=vectoradd(x,vecsm!* v);
return x
end;
symbolic procedure vectoradd(u,v);
% Add two vectors or two scalars
begin scalar x,uisvec,visvec;
uisvec := vecp u; visvec := vecp v;
if uisvec and visvec then
<< % Adding two vectors
x :=mkvect(2);
for k:=0:2 do putv(x,k,aeval list('plus, getv(u,k), getv(v,k)));
return x
>>
else if not (uisvec or visvec) then
<< % Adding two scalars
return aeval list('plus, u, v)
>>
else rerror(avector,8,"Type mismatch in VECTORADD");
end;
% Difference of two vectors
vectorfn('difference,'vectordiff);
symbolic procedure vectordiff vargs;
% Vector - Vector
begin scalar x,y;
x := vecsm!* car vargs;
y := vecsm!* list('minus,cadr vargs); % Negate the second operand
return vectoradd(x,y)
end;
% General case of a quotient involving vectors
vectorfn('quotient,'vectorquot);
symbolic procedure vectorquot vargs;
% This code deals with the cases
%
% (1) Vector / scalar
% (2) Vector / (scalar-valued vector expression)
% (3) Scalar / (scalar-valued vector expression)
%
begin scalar vdivisor,vdividend;
vdivisor := aeval cadr vargs;
if vecp vdivisor
then rerror(avector,9,"Attempt to divide by a vector");
vdividend := aeval car vargs;
if threevectorp vdividend then return vectorapply('quotient,
vdividend, list vdivisor)
else return aeval list('quotient, vdividend, vdivisor);
end;
% Vector cross product
vectorfn('cross,'vectorcrossprod);
symbolic procedure vectorcrossprod vargs;
begin scalar x,y,u0,u1,u2,v0,v1,v2,w0,w1,w2;
x := vecsm!* car vargs;
y := vecsm!* cadr vargs;
u0 := getv(x,0); u1 := getv(x,1); u2 := getv(x,2);
v0 := getv(y,0); v1 := getv(y,1); v2 := getv(y,2);
% Calculate each component of the cross product
w0 := aeval list('difference,
list('times,u1,v2),
list('times,u2,v1));
w1 := aeval list('difference,
list('times,u2,v0),
list('times,u0,v2));
w2 := aeval list('difference,
list('times,u0,v1),
list('times,u1,v0));
x := mkvect(2);
putv(x,0,w0); putv(x,1,w1); putv(x,2,w2);
return x
end;
% Vector modulus
vectorfn('vmod,'vectormod);
% There are two definitions due to the existence of a bug in the REDUCE
% code for SQRT : in the IBM version of REDUCE 3.3 an attempt to take
% SQRT of 0 results in an error, so I have coded round it.
% The first version which follows is the succinct version which will
% work if SQRT(0) doesn't give an error.
% symbolic procedure vectormod u;
% aeval list('sqrt,list('dot,car u,car u));
% This version is a little longer but it works even if SQRT(0) doesn't.
symbolic procedure vectormod u;
begin scalar v;
v := aeval list('dot, car u, car u);
if zerop v then return 0
else return aeval list('sqrt,v);
end;
% Vector dot product
vectorfn('dot,'vectordot);
symbolic procedure vectordot vargs;
begin scalar x,y,u0,u1,u2,v0,v1,v2;
x := car vargs;
y := cadr vargs;
u0 := getv(x,0); u1 := getv(x,1); u2 := getv(x,2);
v0 := getv(y,0); v1 := getv(y,1); v2 := getv(y,2);
% Calculate the scalar product
return aeval list('plus,
list('times,u0,v0),
list('times,u1,v1),
list('times,u2,v2))
end;
% Component-wise assignment of a vector (AVEC)
vectorfn('avec,'vectoravec);
symbolic procedure vectoravec vargs;
begin scalar x;
% Build a vector from the argument list
if not eqn(length(vargs),3) then
rerror(avector,10,"Incorrect number of args in AVEC");
x := mkvect(2);
putv(x,0,aeval first vargs);
putv(x,1,aeval second vargs);
putv(x,2,aeval third vargs);
return x
end;
% Gradient of a scalar
vectorfn('grad,'vectorgrad);
symbolic procedure vectorgrad vargs;
begin scalar x,y;
x := mkvect(2);
y := aeval car vargs;
putv(x,0,aeval list('quotient,
list('df,y,first !*coords),
!*hfac 0));
putv(x,1,aeval list('quotient,
list('df,y,second !*coords),
!*hfac 1));
putv(x,2,aeval list('quotient,
list('df,y,third !*coords),
!*hfac 2));
return x
end;
% Divergence of a vector
vectorfn('div,'vectordiv);
symbolic procedure vectordiv vargs;
begin scalar x,u0,u1,u2;
x := vecsm!* car vargs;
u0 := getv(x,0); u1 := getv(x,1); u2 := getv(x,2);
u0 := aeval list('times,u0,!*hfac 1,!*hfac 2);
u1 := aeval list('times,u1,!*hfac 0,!*hfac 2);
u2 := aeval list('times,u2,!*hfac 0,!*hfac 1);
x := aeval list('plus,
list('df,u0,first !*coords),
list('df,u1,second !*coords),
list('df,u2,third !*coords));
x := aeval list('quotient,x,list('times,
!*hfac 0,!*hfac 1,!*hfac 2));
return x
end;
% Curl of a vector
vectorfn('curl,'vectorcurl);
symbolic procedure vectorcurl vargs;
begin scalar x,u0,u1,u2,v0,v1,v2,w0,w1,w2;
x := vecsm!* car vargs;
u0 := aeval list('times,getv(x,0),!*hfac 0);
u1 := aeval list('times,getv(x,1),!*hfac 1);
u2 := aeval list('times,getv(x,2),!*hfac 2);
v0 := first !*coords; v1 := second !*coords; v2 := third !*coords;
x := mkvect(2);
w0 := aeval list('times,
list('difference,
list('df,u2,v1),
list('df,u1,v2)),
!*hfac 0);
w1 := aeval list ('times,
list('difference,
list('df,u0,v2),
list('df,u2,v0)),
!*hfac 1);
w2 := aeval list('times,
list('difference,
list('df,u1,v0),
list('df,u0,v1)),
!*hfac 2);
putv(x,0,w0); putv(x,1,w1); putv(x,2,w2);
x := aeval list('quotient,
x,
list('times, !*hfac 0, !*hfac 1, !*hfac 2));
return x
end;
% Del-squared (Laplacian) of a scalar or vector
vectorfn('delsq,'vectordelsq);
symbolic procedure vectordelsq vargs;
begin scalar x,y,v0,v1,v2,w0,w1,w2;
x := vecsm!* car vargs;
if vecp x then
% Cunning definition of Laplacian of a vector in terms of
% grad, div and curl
return aeval list('difference,
list('grad, list('div,x)),
list('curl, list('curl,x)))
else
<< % Laplacian of a scalar ... which simply requires lots of
% calculus
if null x then x := car vargs;
y := aeval list('times,!*hfac 0, !*hfac 1, !*hfac 2);
v0 := first !*coords;
v1 := second !*coords;
v2 := third !*coords;
w0 := aeval list('df,
list('quotient,
list('times,
!*hfac 1,
!*hfac 2,
list('df,x,v0)),
!*hfac 0),
v0);
w1 := aeval list('df,
list('quotient,
list('times,
!*hfac 2,
!*hfac 0,
list('df,x,v1)),
!*hfac 1),
v1);
w2 := aeval list('df,
list('quotient,
list('times,
!*hfac 0,
!*hfac 1,
list('df,x,v2)),
!*hfac 2),
v2);
return aeval list('quotient,
list('plus,w0,w1,w2),
y)
>>;
end;
% Vector substitution - definition of SUB as a VECTORFN
% function.
vectorfn('sub,'vectorsub);
% Now we have to define mapping for SUB. It's made a little complicated
% by the fact that the argument list for SUB has the interesting bit
% i.e. the vector, at the end.
symbolic procedure vectorsub vargs;
begin scalar subslist,vexpr,x;
vexpr := car reverse vargs; % That was the easy part !
% Now we need to get the rest of the list
subslist := reverse cdr reverse vargs; % Dirty, but effective !
x := mkvect(2);
for k:=0:2 do
putv(x,k,aeval append('(sub),
append(subslist,list getv(vexpr,k))));
return x
end;
% Component-wise application of a scalar operation to a vector
symbolic procedure vectorapply(vecop,v,args);
begin scalar vv,x,y;
x := mkvect(2);
vv := if not vectorp v then vecsm!* v
else v;
for k:=0:2 do
<< % Apply the operation to each component
y := getv(vv,k);
y := if null args then aeval list(vecop,y)
else aeval append(list(vecop,y),args);
putv(x,k,y);
>>;
return x
end;
% We need to define a dummy routine to apply a function to a scalar
symbolic procedure scalarapply(op,v,args);
if null args then aeval list(op,v)
else aeval append(list(op,v),args);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% %
% COORDINATE SYSTEM MODULE %
% %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% We begin with a function which declares the names of the coordinates
% to be used
symbolic procedure coordinates u;
begin scalar x;
if not eqn(length(u),3)
then rerror(avector,11,"Wrong number of args");
for each y in u do if (x := gettype y) and not(x eq 'operator) then
rerror(avector,12,"Name declared as coordinate is not a kernel");
remflag(!*coords,'reserved);
!*coords := u;
x := aeval list('avec,first u,second u,third u);
remflag('(coords),'reserved);
setk('coords,x);
flag('(coords),'reserved);
%flag(u,'reserved);
return u
end;
symbolic operator coordinates;
algebraic procedure scalefactors(h1,h2,h3);
begin
remflag('(dvolume hfactors),'reserved);
hfactors := avec(h1,h2,h3);
dvolume := h1*h2*h3;
flag('(dvolume hfactors),'reserved);
end;
% We define a procedure that extracts the n-th scale factor
symbolic procedure !*hfac n;
if not fixp n or n<0 or n>2 then rerror(avector,13,"Invalid index")
else getv(getavalue 'hfactors,n);
% Now we define two useful operators that allow us to define and
% refer to a coordinate system by name
symbolic procedure getcsystem u;
begin scalar x,y;
if not atom u
then rerror(avector,14,"Invalid name for coordinate system");
if not flagp(u,'coordinatesystem)
then rerror(avector,15,"Unknown system");
x := get(u,'coordinates);
y := get(u,'scalefactors);
if x and y then
<< % Put the coordinates and scalefactors in place
remflag(!*coords,'reserved);
!*coords := x;
remflag('(coords),'reserved);
setk('coords,aeval list('avec,first x, second x, third x));
flag('(coords),'reserved);
%flag(x,'reserved);
put('hfactors,'avalue,list('!3vector,y));
remflag('(dvolume),'reserved);
setk('dvolume,aeval list('times, !*hfac 0, !*hfac 1, !*hfac 2));
flag('(dvolume),'reserved);
return x
>>
else rerror(avector,16,"Incompletely specified coordinate system")
end;
symbolic procedure putcsystem u;
begin
if not atom u
then rerror(avector,17,"Invalid name for coordinate system");
flag(list u,'coordinatesystem);
put(u,'coordinates,!*coords);
put(u,'scalefactors,getavalue 'hfactors);
!*csystems := union(list u,!*csystems);
return u
end;
deflist('((coordinates rlis)),'stat);
!*coords := '(x y z);
!*csystems := nil;
% The following procedure calculates the derivative of a vector
% function of a scalar variable, including the scale factors in
% the coefficients.
% symbolic operator vecdf;
% Commented out by M.MacCallum or surfint fails
flag('(vecdf), 'vectorfn);
% To replace previous line - M. MacCallum
vectorfn('vecdf,'vectordf);
symbolic procedure vectordf u;
begin scalar v,idv,x;
v := vecsm!* car u;
idv := cadr u;
if not vecp v then rerror(avector,18,"First arg is not a vector");
if not atom idv then rerror(avector,19,"Second arg is not an atom");
x := mkvect(2);
for k:=0:2 do
<< % Calculate components one at a time
putv(x,k,aeval list('times,
!*hfac k,
list('df,
getv(v,k),
idv)))
>>;
return x
end;
% We define three popular curvilinear coordinate systems :
% Cartesian, spherical polar and cylindrical
algebraic;
vec coords,hfactors;
% flag('(coords hfactors),'reserved); % Interferes with EXCALC.
infix dot,cross;
precedence dot,*;
precedence cross,*;
coordinates x,y,z;
scalefactors(1,1,1);
putcsystem 'cartesian;
coordinates r,theta,phi;
scalefactors(1,r,r*sin(theta));
putcsystem 'spherical;
coordinates r,z,phi;
scalefactors(1,1,r);
putcsystem 'cylindrical;
% And we choose to use Cartesians initially ...
getcsystem 'cartesian;
% Extensions to REDUCE vector package
% Definite-integral routine ... trivially simple
algebraic procedure defint(fn,x,xlower,xupper);
begin scalar indefint;
indefint:=int(fn,x);
return sub(x=xupper,indefint)-sub(x=xlower,indefint)
end;
vectormapping 'defint; % DEFINT is now a vector function too
% Component-extraction utility - allows us to get components
% of vectors which are arguments to algebraic procedures
symbolic procedure component(v,n);
if not vecp v then rerror(avector,20,"Argument is not a vector")
else getv(vecsm!* v,n);
symbolic operator component,vecp;
algebraic procedure volintegral(fn,vlower,vupper);
begin scalar integrand,idpvar,xlower,xupper,kindex;
integrand := fn*hfactors(0)*hfactors(1)*hfactors(2);
for k:=0:2 do
<< % Perform each integration separately. The order of integration
% is determined by the control vector VOLINTORDER
kindex := volintorder(k);
idpvar := coords(kindex);
xlower := component(vlower,kindex);
xupper := component(vupper,kindex);
integrand := defint(integrand,idpvar,xlower,xupper);
>>;
return integrand
end;
% Define the initial setting of VOLINTORDER
volintorder := avec(0,1,2);
% Line integral
algebraic procedure lineint(v,curve,ivar);
begin scalar scalfn,vcomp,hcomp,dcurve;
scalfn := 0;
for k:=0:2 do
<< % Form the integrand
vcomp := component(v,k);
hcomp := hfactors(k);
dcurve := df(component(curve,k),ivar);
scalfn := scalfn + vcomp*hcomp*dcurve
>>;
scalfn := vecsub(coords,curve,scalfn) ; % Added by M. MacCallum
return int(scalfn,ivar)
end;
algebraic procedure deflineint(v,curve,ivar,ilb,iub);
begin scalar indfint;
indfint := lineint(v,curve,ivar);
return sub(ivar=iub,indfint)-sub(ivar=ilb,indfint)
end;
% Attempt to implement dot and cross as single-character infix
% operators upon vectors
symbolic;
% Cross-product is easy : we simply tell Reduce that up-arrow is a
% synonym for CROSS
newtok '((!^) cross);
% Dot is more difficult : the period (.) is already defined as the
% CONS function, and unfortunately REVAL1 spots this before it
% checks the type of the arguments, so declaring CONS to be
% VECTORMAPPING won't work. What is required is a hack to the
% routine that carries out CONS at SYMBOLIC level.
% We now redefine RCONS which is the function invoked when CONS is used
% in Reduce.
remflag('(rcons),'lose); % We must use this definition.
symbolic procedure rcons u;
begin scalar x,y;
argnochk ('cons . u);
if (y := getrtype(x := reval cadr u)) eq 'vector
then return mk!*sq simpdot u
% The following line was added to enable . to be used as vector product
% (Amended by M. MacCallum)
else if (y eq '!3vector)
then return apply('vectordot, list(mapcar(u, 'vecsimp!*)) )
else if not(y eq 'list) then typerr(x,"list")
else return 'list . reval car u . cdr x
end;
vectorfn('cons,'vectordot);
% Rest added by M. MacCallum
flag('(surfint vecsub),'vectorfn);
vectorfn('surfint,'vsurfint);
symbolic procedure vsurfint vargs;
begin scalar sivar1, sivar2, sivar3, sivar4, sivar5 ;
if not (length vargs = 8) then rerror(avector,21,
"Wrong number of args to SURFINT");
if not (vecp(sivar1 := car vargs) and
vecp(sivar2 := cadr vargs) and
idp car(sivar3 := cddr vargs) and
idp car(sivar4 := cdddr sivar3))
then rerror(avector,22,
"Wrong type(s) of arguments supplied to SURFINT");
sivar2 := vecsm!* sivar2 ;
sivar3 := reverse cdddr reverse sivar3 ;
sivar5 := aeval list('cross,
list('vecdf, sivar2, car sivar3),
list('vecdf, sivar2, car sivar4)) ;
sivar1 := vecsm!* sivar1 ;
sivar5 := aeval list('dot, sivar1, sivar5) ;
sivar5 := aeval list('vecsub,'coords,sivar2,sivar5);
return aeval append(list('defint,
append(list('defint, sivar5), sivar3)),
sivar4) ;
end ;
vectorfn('vecsub,'vvecsub);
symbolic procedure vvecsub vargs ;
begin scalar vsarg1, vsarg2, vsarg3;
if not (length vargs = 3) then rerror(avector,23,
"Wrong number of arguments to VECSUB");
if not (vecp car vargs and vecp cadr vargs) then
rerror(avector,24,
"First two arguments to VECSUBS must be vectors");
vsarg1 := vecsm!* car vargs;
vsarg2 := vecsm!* cadr vargs;
vsarg3 := caddr vargs;
if not vecp vsarg3 then vsarg3 := prepsq cadr vsarg3;
return aeval list('sub,
list('equal, !*a2k component(vsarg1, 0),
component(vsarg2, 0)),
list('equal, !*a2k component(vsarg1, 1),
component(vsarg2, 1)),
list('equal, !*a2k component(vsarg1, 2),
component(vsarg2, 2)),
vsarg3);
end ;
endmodule;
end;