reduce-historical/r38/packages/orthovec/orthovec.red

module orthovec;  % 3-D vector calculus package.

create!-package('(orthovec),'(contrib avector));

%          %========================================%
%          %              ORTHOVEC                  %
%          %========================================%
%          %     A 3-D VECTOR CALCULUS PACKAGE      %
%          %     USING ORTHOGONAL CURVILINEAR       %
%          %            COORDINATES                 %
%          %                                        %
%          %   copyright James W Eastwood,          %
%          %             Culham Laboratory,         %
%          %             Abingdon, Oxon.            %
%          %                                        %
%          %             February 1987              %
%          %                                        %
%          %   This new version differs from the    %
%          %   original version published in CPC,   %
%          %   47(1987)139-147 in the following     %
%          %   respects:                            %
%          %                                        %
%          %   *.+.,etc replaced by +,-,*,/         %
%          %   *unary vector +,-,/ introduced       %
%          %   *vector component selector _         %
%          %   *general tidy up                     %
%          %   *L'Hopitals rule in Taylor series    %
%          %   *extended division definition        %
%          %   *algebraic output of lisp vectors    %
%          %   *exponentiation of vectors           %
%          %   *vector extension of depend          %
%          %                                        %
%          %            Version 2                   %
%          %        All rights reserved             %
%          %     copyright James W Eastwood         %
%          %             June 1990                  %
%          %                                        %
%          %   This is a preliminary version of     %
%          %   the NEW VERSION of ORTHOVEC which    %
%          %   will be available from the Computer  %
%          %   Physics Communications Program       %
%          %   Library, Dept. of Applied Maths and  %
%          %   Theoretical Physics, The Queen's     %
%          %   University of Belfast, Belfast       %
%          %   BT7 1NN, Northern Ireland.           %
%          %   See any copy of CPC for further      %
%          %   details of the library services.     %
%          %                                        %
%          %========================================%
%          %       REDUCE 3.4 is assumed            %
%          %========================================%
%          
%
%
%-------------------------------------------------------------------
%                       INITIALISATION
%%

algebraic;

%select coordinate system
%========================
procedure vstart0;
begin scalar ctype;
write "Select Coordinate System by number";
write "1] cartesian";
write "2] cylindrical";
write "3] spherical";
write "4] general";
write "5] others";
%remove previous settings
clear u1,u2,u3,h1,h2,h3;
depend h1,u1,u2,u3;
depend h2,u1,u2,u3;
depend h3,u1,u2,u3;
nodepend h1,u1,u2,u3;
nodepend h2,u1,u2,u3;
nodepend h3,u1,u2,u3;
%select coordinate system
ctype := symbolic read();
if ctype=1 then
 << u1:=x;u2:=y;u3:=z;h1:=1;h2:=1;h3:=1 >>
 else if ctype=2 then
  << u1:=r;u2:=th;u3:=z;h1:=1;h2:=r;h3:=1 >>
  else if ctype=3 then
   << u1:=r;u2:=th;u3:=ph;h1:=1;h2:=r;h3:=r*sin(th) >>
   else if ctype=4 then
    << depend h1,u1,u2,u3;depend h2,u1,u2,u3;depend h3,u1,u2,u3 >>
    else << write "To define another coordinate system, give values ";
    write "to components u1,u2,u3 and give functional form or";
    write "DEPEND for scale factors h1,h2 and h3. For example,";
write "to set up paraboloidal coords u,v,w type in:-";
write "u1:=u;u2:=v;u3:=w;h1:=sqrt(u**2+v**2);h2:=h1;h3:=u*v;">>;
write "coordinate type = ",ctype;
write "coordinates = ",u1,",",u2,",",u3;
write "scale factors = ",h1,",",h2,",",h3;
return
end$
let vstart=vstart0()$

%give access to lisp vector procedures
%=======================================
symbolic operator putv,getv,mkvect;
flag('(vectorp), 'direct);
flag('(vectorp), 'boolean);

%-------------------------------------------------------------------
%                      INPUT-OUTPUT
%
%set a new vector
%===================
procedure svec(c1,c2,c3);
begin scalar a;a:=mkvect(2);
putv(a,0,c1);putv(a,1,c2);putv(a,2,c3); 
return a
end$

%output a vector
%===============
procedure vout(v);
begin;
if vectorp(v) then
  for j:=0:2 do write "[",j+1,"] ",getv(v,j)
  else write v;
return v
end$

%-------------------------------------------------------------------
%               REDEFINITION OF SOME STANDARD PROCEDURES
%
% Vector extension of standard definitions of depend and nodepend.

remflag('(depend nodepend),'lose);   % We must use these definitions.

symbolic procedure depend u;
begin scalar v,w; v:= !*a2k car u; 
for each x in cdr u do
if vectorp(v) then
 for ic:=0:upbv(v) do 
 <<if atom(w:=getv(v,ic)) and not numberp(w) then depend1(w,x,t)>>
 else depend1(car u,x,t)
end$

symbolic procedure nodepend u;
begin scalar v,w;
rmsubs();
v:= !*a2k car u; 
for each x in cdr u do
if vectorp(v) then
 for ic:=0:upbv(v) do 
 <<if atom(w:=getv(v,ic)) and not numberp(w) then depend1(w,x,nil)>>
 else depend1(car u,x,nil)
end $

%
%-------------------------------------------------------------------
%                      ALGEBRAIC OPERATIONS
%
%define symbols for vector algebra
%=====================================
newtok '(( !+ ) vectoradd);
newtok '(( !- ) vectordifference);
newtok '((!> !< ) vectorcross);
newtok '(( !* ) vectortimes);
newtok '(( !/ ) vectorquotient);
newtok '(( !_ ) vectorcomponent);
newtok '(( !^ ) vectorexpt);
%
%define operators
%================
operator vectorminus,vectorplus,vectorrecip;
infix vectoradd,vectordifference,vectorcross,vectorexpt,
      vectorcomponent,vectortimes,vectorquotient,dotgrad;
precedence vectoradd,<;
precedence vectordifference,vectoradd;
precedence dotgrad,vectordifference;
precedence vectortimes,dotgrad;
precedence vectorcross,vectortimes;
precedence vectorquotient,vectorcross;
precedence vectorexpt,vectorquotient;
precedence vectorcomponent,vectorexpt;

deflist( '(
   (vectordifference vectorminus)
   (vectoradd vectorplus)
   (vectorquotient vectorrecip)
   (vectorrecip vectorrecip)
 ), 'unary)$

deflist('((vectorminus vectorplus) (vectorrecip vectortimes)),
		 'alt)$

%extract component of a vector
%=============================
procedure vectorcomponent(v,ic);
if vectorp(v) then
  if ic=1 or ic=2 or ic=3 then getv(v,ic-1)
    else rerror(orthovec,1,"Incorrect component number")
  else rerror(orthovec,2,"Not a vector")$

%
%add vector or scalar pair v1 and v2
%===================================
procedure vectoradd(v1,v2);
begin scalar v3;
if vectorp(v1) and vectorp(v2) then
  <<v3:=mkvect(2);
  for ic:=0:2 do putv(v3,ic,plus(getv(v1,ic),getv(v2,ic)))>>
  else
  if not(vectorp(v1)) and not(vectorp(v2)) then v3:=plus(v1, v2)
     else rerror(orthovec,3,"Incorrect args to vector add");
return v3
end$

%unary plus
%==========
procedure vectorplus(v);v$
%
%negate vector or scalar v
%=========================
procedure vectorminus(v); 
begin scalar v3;
if vectorp(v) then
  <<v3:=mkvect(2);
  for ic:=0:2 do putv(v3,ic,minus(getv(v,ic)))>>
  else v3:=minus(v);
return v3
end$

%scalar or vector subtraction
%============================
procedure vectordifference(v1,v2);(v1 + vectorminus(v2))$

%dot product or scalar times
%===========================
procedure vectortimes(v1,v2);
begin scalar v3;
if vectorp(v1) and vectorp(v2) then
  v3:= for ic:=0:2 sum times(getv(v1,ic),getv(v2,ic))
  else
  if not(vectorp(v1)) and not(vectorp(v2)) then
    v3:=times(v1 , v2 )
    else if vectorp(v1) and not(vectorp(v2)) then
      <<v3:=mkvect(2);
      for ic:=0:2 do putv(v3,ic,times(getv(v1,ic),v2)) >>
      else
      <<v3:=mkvect(2);
      for ic:=0:2 do putv(v3,ic,times(getv(v2,ic),v1)) >>;
return v3
end$

%vector cross product
%====================
procedure vectorcross(v1,v2);
begin scalar v3;
if vectorp(v1) and vectorp(v2) then
  <<v3:=mkvect(2);
  putv(v3,0,getv(v1,1)*getv(v2,2)-getv(v1,2)*getv(v2,1));
  putv(v3,1,getv(v1,2)*getv(v2,0)-getv(v1,0)*getv(v2,2));
  putv(v3,2,getv(v1,0)*getv(v2,1)-getv(v1,1)*getv(v2,0))>>
  else rerror(orthovec,4,"Incorrect args to vector cross product");
return v3
end$

%vector division
%===============
procedure vectorquotient(v1,v2);
if vectorp(v1) then
 if vectorp(v2) then
  quotient (v1*v2,v2*v2)
  else v1*recip(v2)
 else if vectorp(v2) then
  v1*v2*recip(v2*v2)
  else quotient(v1,v2)$

procedure vectorrecip(v);
if vectorp(v) then
 v*recip(v*v)
 else recip(v)$

%length of vector
%================
procedure vmod(v);sqrt(v * v)$

%vector exponentiation
%=====================
procedure vectorexpt(v,n);
if vectorp(v) then expt(vmod(v),n) else expt(v,n)$

%-------------------------------------------------------------------
%                      DIFFERENTIAL OPERATIONS
%

%div
%===
procedure div(v);
if vectorp(v) then
  (df(h2*h3*getv(v,0),u1)+df(h3*h1*getv(v,1),u2)
    +df(h1*h2*getv(v,2),u3))/h1/h2/h3
  else rerror(orthovec,5,"Incorrect arguments to div")$

%grad
%====
procedure grad(s);
begin scalar v;
v:=mkvect(2);
if vectorp(s) then 
  rerror(orthovec,6,"Incorrect argument to grad")
  else << putv(v,0,df(s,u1)/h1);
          putv(v,1,df(s,u2)/h2);
          putv(v,2,df(s,u3)/h3) >>;
return v
end$

%curl
%====
procedure curl(v);
begin scalar v1;
v1:=mkvect(2);
if vectorp(v) then
  << putv(v1,0,(df(h3*getv(v,2),u2)-df(h2*getv(v,1),u3))/h2/h3);
     putv(v1,1,(df(h1*getv(v,0),u3)-df(h3*getv(v,2),u1))/h3/h1);
     putv(v1,2,(df(h2*getv(v,1),u1)-df(h1*getv(v,0),u2))/h1/h2) >>
  else rerror(orthovec,7,"Incorrect argument to curl");
return v1
end$

%laplacian
%=========
procedure delsq(v);
if vectorp(v) then (grad(div(v)) - curl(curl(v))) else div(grad(v))$

%differentiation
%===============
procedure vdf(v,x);
begin scalar v1;
if vectorp(x) then
  rerror(orthovec,8,"Second argument to VDF must be scalar")
  else if vectorp(v) then
    <<v1:=mkvect(2);for ic:=0:2 do putv(v1,ic,vdf(getv(v,ic),x)) >>
    else v1:=df(v,x);
return v1
end$

%v1.grad(v2)
%===========
procedure dotgrad(v1,v2);
if vectorp(v1) then
  if vectorp(v2) then
    (1/2)*(grad(v1 * v2) + v1 * div(v2)  - div(v1) * v2
          - (curl(v1 >< v2) + v1 >< curl(v2) - curl(v1) >< v2 ))
    else v1 * grad(v2)
  else rerror(orthovec,9,"Incorrect arguments to dotgrad")$

%3-D Vector Taylor Expansion about vector point
%==============================================
procedure vtaylor(vex,vx,vpt,vorder);
%note: expression vex, variable vx, point vpt and order vorder
%      are any legal mixture of vectors and scalars
begin scalar vseries;
if vectorp(vex) then
  <<vseries:=mkvect(2);
  for ic:=0:2 do putv(vseries,ic,vptaylor(getv(vex,ic),vx,vpt,vorder))>>
  else vseries:=vptaylor(vex,vx,vpt,vorder);
return vseries
end$

%Scalar Taylor expansion about vector point
%==========================================
procedure vptaylor(sex,vx,vpt,vorder);
%vector variable
if vectorp(vx) then
  if vectorp(vpt) then
%vector order
    if vectorp(vorder) then
      taylor( taylor( taylor( sex,
      getv(vx,0), getv(vpt,0), getv(vorder,0) ),
      getv(vx,1), getv(vpt,1), getv(vorder,1) ),
      getv(vx,2), getv(vpt,2), getv(vorder,2) )
      else
      taylor( taylor( taylor( sex,
      getv(vx,0), getv(vpt,0), vorder),
      getv(vx,1), getv(vpt,1), vorder),
      getv(vx,2), getv(vpt,2), vorder)
    else rerror(orthovec,10,"VTAYLOR: vector VX mismatches scalar VPT")
%scalar variable
  else if vectorp(vpt) then
    rerror(orthovec,11,"VTAYLOR: scalar VX mismatches vector VPT")
    else if vectorp(vorder) then
      rerror(orthovec,12,"VTAYLOR: scalar VX mismatches vector VORDER")
      else taylor(sex,vx,vpt,vorder)$
 
%Scalar Taylor expansion of ex wrt x about point pt to order n
%=============================================================
procedure taylor(ex,x,pt,n);
begin scalar term,series,dx,mfac;
if numberp n then <<
  mfac:=1;dx:=x-pt;term:=ex;
  series:= limit(ex,x,pt) +
  for k:=1:n sum limit((term:=df(term,x)),x,pt)*(mfac:=mfac*dx/k) >>
  else rerror(orthovec,13,
	      "Truncation orders of Taylor series must be integers");
return series
end$
%
%limiting value of exression ex as x tends to pt
%===============================================
procedure limit(ex,x,pt);
begin scalar lim,denex,numex;
%polynomial
lim:=if (denex:=den(ex))=1 then sub(x=pt,ex)
 else
%zero denom rational
 if sub(x=pt,denex)=0 then
%l'hopital's rule
  << if sub(x=pt,(numex:=num(ex)))=0 then
   limit(df(numex,x)/df(denex,x),x,pt)
%singular
   else rerror(orthovec,14,"Singular coefficient found by LIMIT")>>
%nonzero denom rational
  else sub(x=pt,ex);
return lim
end$
%
%-------------------------------------------------------------------
%                      INTEGRAL OPERATIONS
%
% Vector Integral
%================
procedure vint(v,x);
begin scalar v1;
if vectorp(x) then
  rerror(orthovec,15,"Second argument to VINT must be scalar")
  else if vectorp(v) then
    <<v1:=mkvect(2);for ic:=0:2 do putv(v1,ic,int(getv(v,ic),x)) >>
    else v1:=int(v,x);
return v1
end$

%Definite Vector Integral
%========================
procedure dvint(v,x,xlb,xub);
begin scalar integr,intval;
if vectorp(xlb) or vectorp(xub) 
  then rerror(orthovec,16,"Limits to DVINT must be scalar")
  else if vectorp(v) then
    <<intval:=mkvect(2);
    for ic:=0:2 do <<integr:=int(getv(v,ic),x);
    putv(intval,ic,sub(x=xub,integr)-sub(x=xlb,integr)) >> >>
    else
    <<integr:=int(v,x);intval:=sub(x=xub,integr)-sub(x=xlb,integr)>>;
return intval
end$

%Volume Integral 
%===============
procedure volint(v);
begin scalar v1;
if vectorp(v) then
  <<v1:=mkvect(2);for ic:=0:2 do putv(v1,ic,volint(getv(v,ic))) >>
  else v1:= int( int( int(v*h1*h2*h3,u1),u2),u3);
return v1
end$

%Definite Volume Integral
%========================
procedure dvolint(v,vlb,vub,n);
begin scalar v1,intgrnd;
if vectorp(vlb) and vectorp(vub) then
  <<intgrnd:= (h1*h2*h3) * v;
  v1:= if n=1 then
    dvint(dvint(dvint(intgrnd,
    u1,getv(vlb,0),getv(vub,0)),
    u2,getv(vlb,1),getv(vub,1)),
    u3,getv(vlb,2),getv(vub,2) )
      else if n=2 then
      dvint(dvint(dvint(intgrnd,
      u3,getv(vlb,2),getv(vub,2)),
      u1,getv(vlb,0),getv(vub,0)),
      u2,getv(vlb,1),getv(vub,1) )
        else if n=3 then
        dvint(dvint(dvint(intgrnd,
        u2,getv(vlb,1),getv(vub,1)),
        u3,getv(vlb,2),getv(vub,2)),
        u1,getv(vlb,0),getv(vub,0) )
          else if n=4 then
          dvint(dvint(dvint(intgrnd,
          u1,getv(vlb,0),getv(vub,0)),
          u3,getv(vlb,2),getv(vub,2)),
          u2,getv(vlb,1),getv(vub,1) )
            else if n=5 then
            dvint(dvint(dvint(intgrnd,
            u2,getv(vlb,1),getv(vub,1)),
            u1,getv(vlb,0),getv(vub,0)),
            u3,getv(vlb,2),getv(vub,2) )
              else 
              dvint(dvint(dvint(intgrnd,
              u3,getv(vlb,2),getv(vub,2)),
              u2,getv(vlb,1),getv(vub,1)),
              u1,getv(vlb,0),getv(vub,0)) >>
  else rerror(orthovec,17,"Bounds to DVOLINT must be vectors");
return v1
end$

%Scalar Line Integral
%====================
procedure lineint(v,vline,tt);
if vectorp(v) and vectorp(vline) and not vectorp(tt) then
 int(sub( u1=getv(vline,0), u2=getv(vline,1), u3=getv(vline,2), 
  getv(v,0) * df(getv(vline,0),tt) *  h1 +
  getv(v,1) * df(getv(vline,1),tt) *  h2 +
  getv(v,2) * df(getv(vline,2),tt) *  h3 ) , tt)
 else rerror(orthovec,18,"Incorrect arguments to LINEINT")$

%Definite Scalar Line Integral
%=============================
procedure dlineint(v,vline,tt,tlb,tub);
begin scalar integr,intval;
if vectorp(tlb) or vectorp(tub) 
  then rerror(orthovec,19,"Limits to DLINEINT must be scalar")
  else <<integr:=lineint(v,vline,tt);
  intval:=sub(tt=tub,integr)-sub(tt=tlb,integr)>>;
return intval
end$
 
%
%-------------------------------------------------------------------
%        SET DEFAULT COORDINATES TO CARTESIAN
%
% write "Cartesian coordinates selected by default";
% write "If you wish to change this then type VSTART";
% write "and follow the instructions given.";
% write "u1,u2,u3 are reserved for coordinate names";
% write "h1,h2,h3 are reserved for scale factor names";
ctype:=1$u1:=x$u2:=y$u3:=z$h1:=1$h2:=1$h3:=1$
% write "coordinate type = ",ctype;
% write "coordinates = ",u1,",",u2,",",u3;
% write "scale factors = ",h1,",",h2,",",h3;

%-------------------------------------------------------------------

endmodule;

end;