Check-ins Using Download Hex

Artifact 5e1e57e494dee6627126f87322c7dec12a0f7ccd2482fcbf6d11ba2027c26d37:

module expres$

% Author: R. Liska

% Version REDUCE 3.6    05/1991

global '(!*outp)$                  % declarations for 3.4
fluid '(!*wrchri orig!* posn!*)$

switch wrchri$

global '(olddimension!* dimension!* coordindx!* cyclic1!* cyclic2!*
         sfprod!* nscal!*)$

flag('(share),'eval)$    % So that SHARE recognized by FASL.

share olddimension!*,dimension!*,coordindx!*,cyclic1!*,cyclic2!*,
         sfprod!*,nscal!*$

!*precise := nil;   % Needed in this module.

nscal!*:=0$

put('tensor,'tag,'tens)$
put('tensor,'fn,'tensfn)$
put('tensor,'evfn,'expres)$
put('tens,'prifn,'tenspri)$
flag('(tensor),'sprifn)$
put('tens,'setprifn,'settenspri)$
put('tensor,'typeletfn,'tenslet)$

symbolic procedure ptensor x$
   'tensor$

symbolic procedure poptensor u$
   if null u then 'tensor else nil$

deflist('((tensor ptensor)
          (tensop poptensor)
          (df getrtypecar)
          (diff getrtypecar)
           (!& ptensor)
          (!# ptensor)
          (!? ptensor)
          (grad ptensor)
          (div ptensor)
          (lapl ptensor)
          (curl ptensor)
          (vect ptensor)
          (dyad ptensor)
          (dirdf ptensor)),'rtypefn)$
put('cons,'rtypefn,'getrtypecons)$
put('rcons,'evfn,'evrcons)$
remprop('cons,'psopfn)$

symbolic procedure getrtypecons u$
if getrtypecar u eq 'tensor then 'tensor
  else 'rcons$

symbolic procedure evrcons(u,v)$
rcons cdr u$

symbolic procedure tensor u$
for each a in u do
  <<put(a,'rtype,'tensop)$
    put(a,'avalue,list('tensor,mktensor(0 , (0 . 1)))) >>$

deflist('((tensor rlis)),'stat)$

symbolic procedure tenslet(u,v,typu,b,typv)$
if not atom u then lprie list(" Non atom tensor variable ",u)
  else if b then
     <<if not eqcar(v,'tensor) then v:= mktensor(0,
             if eqcar(v,'!*sq) then cadr v else simp!* v)$
       rmsubs()$
       put(u,'rtype,'tensop)$
       put(u,'avalue,list('tensor,v)) >>
  else
    <<remprop(u,'rtype)$
      remprop(u,'avalue) >>$

%======================================================================
%  Data structure for tensor quantities
%  ====================================
%  (tensor nr rnk val car !*sqvar!*)
%  nr   -  integer, should be equal to actual nscal!*, otherwise
%          the quantity has been defined in previous coor. system
%          number of coordinate system
%  rnk  -  integer, 0,1,2
%          rank of the tensor
%          0 - scalar
%          1 - vertor
%          2 - dyad (matrix)
%  val  -  value
%          s.q.                   for rnk = 0
%          list of s.q.s          for rnk = 1
%          list of lists of s.q.s for rnk = 2
%  !*sqvar!* used in resimplification routine
%======================================================================
%  Smacro definitions for access of data structure subparts
%======================================================================

smacro procedure tensrnk u$
% determines rank from cddr of datastructure
car u$

smacro procedure tensval u$
% determines value from cddr of datastructure
cadr u$

symbolic procedure mktensor(rnk,u)$
'tensor . nscal!* . rnk . u . if !*resubs then !*sqvar!* else list nil$

symbolic procedure settenspri(u,v)$
if not atom u then lprie list(" Non-atom tensor variable ",u)
  else <<prin2!* u$
         prin2!* get('setq,'prtch)$
         tenspri v >>$

symbolic procedure tenspri u$
begin
  scalar rnk$
  u:=cddr u$
  rnk:=car u$
  u:=cadr u$
  if rnk=0 then
      <<pmaprin u$
        terpri!* t >>
    else if rnk=1 then
      <<tenspri1 u$
        orig!*:=0$
        terpri!* t >>
    else if rnk=2 then
      <<prin2!* " ( "$
        tenspri1 car u$
        for each i in cdr u do
          <<prin2!* " , "$
            terpri!* t$
            tenspri1 i >>$
        prin2!* " ) "$
        orig!*:=0$
        terpri!* t >>
    else lprie list(" Can't print tensor ",u," with rank ",rnk)
end$

symbolic procedure tenspri1 u$
<<prin2!* " ( "$
  orig!*:=posn!*$
  pmaprin car u$
  for each i in cdr u do
    <<prin2!* " , "$
      terpri!* t$
      pmaprin i >>$
  orig!*:=orig!* - 3$
  prin2!* " ) " >>$

symbolic procedure pmaprin u$
maprin(!*outp:=prepsq!*  u)$

symbolic procedure updatedimen()$
if olddimension!* = dimension!* then t
  else
    <<if dimension!* =2 then <<coordindx!* :='(2 1)$
                               cyclic1!* :='(1 2)$
                               cyclic2!* :='(2 1) >>
        else if dimension!* =3 then
          <<coordindx!* :='(3 2 1)$
            cyclic1!* :='(2 3 1)$
            cyclic2!* :='(3 1 2) >>
        else lprie list(" Can't handle dimension = ",dimension!*)$
      olddimension!* := dimension!* >>$

symbolic procedure expres(expn,v)$
express expn$

symbolic procedure resimptens u$
mktensor(caddr u,if caddr u=0 then resimp cadddr u
                   else if caddr u=1 then for each a in cadddr u collect
                                                      resimp a
                   else if caddr u=2 then
                     for each a in cadddr u collect
                       for each b in a collect resimp b
                   else lprie list("Can't handle tensor ",u,
                                   " with rank ",caddr u))$

symbolic procedure express expn$
begin
  scalar lst,matrx,rnk,opexpress$
  if not atom expn then go to op$
  if get(expn,'rtype) eq 'tensop and (rnk:=get(expn,'avalue)) and
        car rnk memq '(tensor tensop) and (rnk:=cadr rnk) then return
      if cadr rnk=nscal!* then
          if car cddddr rnk then rnk
            else resimptens rnk
        else lprie list(" You must rebind tensor ",expn,
            " in the new coordinate system")$
  return mktensor(0,simp!* expn)$
op:if car expn = 'vect then return mktensor
       (1,testdim1 for each a in cdr expn collect simp!* a)
     else if car expn = 'dyad then return mktensor
       (2,testdim2 for each a in cdr expn collect
                     for each b in a collect simp!* b)
     else if car expn eq 'tensor then return
       if cadr expn=nscal!* then
           if car cddddr expn then expn
             else resimptens expn
         else lprie list(" You must rebind tensor ",expn,
             " in the new coordinate system")$
  lst:=for each a in cdr expn collect cddr express a$
  if (opexpress:=get(car expn,'express)) then
      <<rnk:=eval(opexpress . list mkquote lst)$
        return mktensor(car rnk,cadr rnk)>>$
  if get(car expn,'simpfn) then return mktensor(0,simp(
      car expn . for each a in lst collect
                       if car a=0 then list('!*sq,cdr a,nil)
                         else typerr(expn," well formed scalar ") ))$
  lst:=for each a in lst collect
               if car a=0 then prepsq cdr a
                 else typerr(expn," well formed tensor")$
  return mktensor(0,!*k2q(car expn.lst))
end$

procedure testdim1 u$
if length u=dimension!* then u
  else <<lprie "Bad number of vector components"$u>>$

procedure testdim2 u$
begin
  scalar x$
  if length u = dimension!* then t
    else go to er$
  x:=u$
a:if length car u = dimension!* then t
    else go to er$
  x:=cdr x$
  if x then go to a$
  return u$
er:lprie "Bad number of dyad components"$
  return u
end$

%======================================================================
%  Procedures in EXPRESS properties of operators are returning
%  (rnk val), their argument is list of (rnk val)

symbolic procedure vectors arglist$
for each i in arglist do
  <<put(i,'rtype,'tensop)$
    put(i,'simpfn,'simpiden)$
    put(i,'avalue,list('tensop,
                  mktensor(1,reverse
                             for each a in coordindx!* collect
                               simp list(i,a) )) )>>$

deflist('((vectors rlis)),'stat)$

symbolic procedure dyads arglist$
for each i in arglist do
  <<put(i,'rtype,'tensop)$
    put(i,'simpfn,'simpiden)$
    put(i,'avalue,list('tensop,
                        mktensor(2,reverse
                                   for each a in coordindx!* collect
                                     reverse
                                     for each b in coordindx!* collect
                                       simp list(i,a,b)))) >>$

deflist('((dyads rlis)),'stat)$

symbolic procedure plusexpress u$
begin
  scalar z$
  z:=car u$
a:u:=cdr u$
  if null u then return z$
  z:=plus2ex(z,car u)$
  go to a
end$

put('plus,'express,'plusexpress)$

symbolic procedure plus2ex(x,y)$
begin
  scalar mtx,mty,slx,sly,rnk,ans,ans1$
  rnk:=tensrnk x$
  if not(rnk=tensrnk y) then lprie "Tensor mishmash"$
  if rnk=0 then return list(rnk,addsq(cadr x,cadr y))
    else if rnk=1 then
             <<slx:=tensval x$
               sly:=tensval y$
               while slx do
                 <<ans:=addsq(car slx,car sly) . ans$
                   slx:=cdr slx$
                   sly:=cdr sly>>$
               ans:= list(1,reverse ans) >>
    else if rnk=2 then
             <<mtx:=tensval x$
               mty:=tensval y$
               while mtx do
                 <<slx:=car mtx$
                   sly:=car mty$
                   ans1:=nil$
                   while slx do
                     <<ans1:=addsq(car slx,car sly) . ans1$
                       slx:=cdr slx$
                       sly:=cdr sly>>$
                   ans:=reverse ans1 . ans$
                   mtx:=cdr mtx$
                   mty:=cdr mty>>$
               ans:=list(2,reverse ans) >>$
  return ans
end$

symbolic procedure timesexpress u$
begin
  scalar z$
  z:=car u$
a:u:=cdr u$
  if null u then return z$
  z:=times2ex(z,car u)$
  go to a
end$

put('times,'express,'timesexpress)$

symbolic procedure times2ex(x,y)$
begin
  scalar rnkx,rnky$
  rnkx:=tensrnk x$
  rnky:=tensrnk y$
  return
    if rnkx=0 then list(rnky,times0ex(tensval x,tensval y,rnky))
      else if rnky=0 then list(rnkx,times0ex(tensval y,tensval x,rnkx))
      else lprie " Tensor mishmash "
end$

symbolic procedure times0ex(x,y,rnk)$
if rnk=0 then multsq(x,y)
  else if rnk=1 then
    for each a in y collect multsq(x,a)
  else if rnk=2 then
    for each a in y collect
      for each b in a collect multsq(x,b)
  else lprie " Tensor mishmash "$

symbolic procedure minusexpress expn$
timesexpress list(list(0,cons(-1,1)),car expn)$

put('minus,'express,'minusexpress)$

symbolic procedure differenceexpress expn$
plusexpress list(car expn,minusexpress list cadr expn)$

put('difference,'express,'differenceexpress)$

symbolic procedure quotientexpress expn$
if tensrnk cadr expn = 0 then
    times2ex(list(0,simp!* list('recip,prepsq tensval cadr expn)),
             car expn)
  else lprie " Tensor mishmash "$

put('quotient,'express,'quotientexpress)$

symbolic procedure exptexpress expn$
if tensrnk car expn=0 and tensrnk cadr expn = 0 then
    list(0,simp!* list('expt,
                       prepsq tensval car expn,
                       prepsq tensval cadr expn))
  else lprie " Tensor mishmash "$

put('expt,'express,'exptexpress)$

symbolic procedure recipexpress expn$
if tensrnk car expn = 0 then
    list(0,simp!* list('recip,
                       prepsq tensval car expn))
  else lprie " Tensor mishmash "$

put('recip,'express,'recipexpress)$

symbolic procedure inprodexpress expn$
begin
  scalar arg1,arg2,rnk1,rnk2$
  arg1:=tensval car expn$
  arg2:=tensval cadr expn$
  rnk1:=tensrnk car expn$
  rnk2:=tensrnk cadr expn$
  return
    if rnk1=1 then inprod1ex(arg1,arg2,rnk2)
      else if rnk1=2 then inprod2ex(arg1,arg2,rnk2)
      else lprie " Tensor mishmash "
end$

put('cons,'express,'inprodexpress)$

symbolic procedure inprod1ex(x,y,rnk)$
begin
  scalar lstx,lsty,mty,z,zz$
  lstx:=x$
  lsty:=y$
  if rnk=1 then
      <<z:=nil . 1$
        while lstx do
          <<z:=addsq(z,multsq(car lstx,car lsty))$
            lstx:=cdr lstx$
            lsty:=cdr lsty>>$
        z:=list(0,z)>>
    else if rnk=2 then
      <<z:=nil$
        lsty:=mty:=copy2(y,t)$
        while car mty do
          <<lstx:=x$
            lsty:=mty$
            zz:=nil . 1$
            while lstx do
              <<zz:=addsq(zz,multsq(car lstx,caar lsty))$
                rplaca(lsty,cdar lsty)$
                lsty:=cdr lsty$
                lstx:=cdr lstx>>$
            z:=nconc(z,list zz) >>$
        z:=list(1,z)>>
    else lprie " Tensor mishmash "$
  return z
end$

symbolic procedure inprod2ex(x,y,rnk)$
begin
  scalar mtx,z$
  mtx:=x$
  if rnk=1 then
      while mtx do
        <<z:=nconc(z,cdr inprod1ex(car mtx,y,1))$
          mtx:=cdr mtx>>
    else if rnk=2 then
      while mtx do
        <<z:=nconc(z,cdr inprod1ex(car mtx,copy2(y,t),2))$
          mtx:=cdr mtx>>
    else lprie " Tensor mishmash "$
  return list(rnk,z)
end$

symbolic procedure outexpress expn$
begin
  scalar x,y,z$
  x:=tensval car expn$
  y:=tensval cadr expn$
  if tensrnk car expn=1 and tensrnk cadr expn=1 and null cddr expn then
      for each i in x do z:=(for each a in y collect multsq(a,i) ) . z
    else lprie list(" Outer product of ",expn)$
  return list(2,reverse z)
end$

put('!&,'express,'outexpress)$
flag('(!&),'tensfn)$

symbolic procedure copy2(x,p)$
if null x then nil else
if p then copy2(car x,nil) . copy2(cdr x,t)
  else car x . copy2(cdr x,nil)$

symbolic procedure listar(arg,j)$
if j=1 then car arg
  else if j=2 then cadr arg
  else if j=3 then caddr arg
  else lprie list(" LISTAR ",arg,j)$

symbolic procedure listarsq(arg,j)$
prepsq listar(arg,j)$

symbolic procedure dinprod expn$
begin
  scalar x,y,z,xx,yy$
  x:=tensval car expn$
  y:=copy2(tensval cadr expn,t)$
  z:=nil . 1$
  if not(tensrnk car expn=2 and tensrnk cadr expn=2 and null cddr expn)
      then lprie list(" D-scalar product of ",expn)$
a:if null x and null y then go to d
    else if null x or null y then go to er$
  xx:=car x$
  yy:=car y$
b:if null xx and null yy then go to c
    else if null xx or null yy then go to er$
  z:=addsq(z,multsq(car xx,car yy))$
  xx:=cdr xx$
  yy:=cdr yy$
  go to b$
c:x:=cdr x$
  y:=cdr y$
  go to a$
d:return list(0,z)$
er:lprie list(" EXPRESS error ",expn," D-S dyads of dif. size")
end$

put('!#,'express,'dinprod)$
put('hash,'express,'dinprod)$
put('hash,'rtypefn,'ptensor)$

symbolic procedure antisymsum(u,v)$
if dimension!* = 2 then difmul(car u,cadr u,cadr v,car v)
  else if dimension!* = 3 then list
      (difmul(cadr u,caddr u,caddr v,cadr v),
       difmul(caddr u,car u,car v,caddr v),
       difmul(car u,cadr u,cadr v,car v))
  else lprie list(" ANTISYMSUM ",u,v)$

symbolic procedure difmul(a,b,c,d)$
% A*C-B*D$
addsq(multsq(a,c),negsq multsq(b,d))$

symbolic procedure vectprod expn$
begin
  scalar x,y,rnx,rny$
  x:=tensval car expn$
  y:=tensval cadr expn$
  rnx:=tensrnk car expn$
  rny:=tensrnk cadr expn$
  if rnx=1 and rny=1 then return
      list(dimension!* - 2,antisymsum(x,y))
    else if rnx=2 and rny=1 then return
      list(dimension!* - 1,for each a in x collect antisymsum(a,y) )
    else if rnx=1 and rny=2 then return
      list(dimension!* - 1,
           if dimension!*=3 then
               tp1 copy2(for each a in tp1(copy2(y,t)) collect
                         antisymsum(x,a),t)
             else for each a in tp1(copy2(y,t)) collect
                         antisymsum(x,a) )
    else lprie list(" VECTPROD of ",expn)
end$

put('!?,'express,'vectprod)$

algebraic operator diff$

symbolic procedure gradexpress expn$
begin
  scalar arg,vt,ans,row,z$
  arg:=tensval car expn$
  vt:=tensrnk car expn$
  if vt=0 then
      for each i in coordindx!* do
        ans:=simp!* list('quotient,
                       list('diff,
                            list('!*sq,arg,nil),
                            getel list('coordinats,i)),
                       getel list('sf,i)) . ans
    else if vt=1 then
      for each i in coordindx!* do
        <<row:=nil$
          for each j in coordindx!* do
            <<z:=list('diff,
                      listarsq(arg,j),
                      getel list('coordinats,i))$
              z:=list list('quotient,
                           z,
                           getel list('sf,i))$
              for each k in coordindx!* do
                z:=list('times,
                        getel list('christoffel,i,j,k),
                        listarsq(arg,k)) . z$
              row:=simp!*('plus.z) . row>>$
          ans:=row . ans>>
    else lprie list(" GRAD of ",expn)$
  return list(vt+1,ans)
end$

put('grad,'express,'gradexpress)$

symbolic procedure divexpress expn$
begin
  scalar arg,vt,ans,z$
  arg:=tensval car expn$
  vt:=tensrnk car expn$
  if vt=1 then
      <<for each i in coordindx!* do
          z:=list('diff,
                  list('times,
                       sfprod!*,
                       listarsq(arg,i),
                       list('recip,
                            getel list('sf,i))),
                  getel list('coordinats,i)) . z$
          z:='plus . z$
          z:=list('quotient,
                  z,
                  sfprod!*)$
          ans:=simp!* z>>
    else if vt=2 then
      for each i in coordindx!* do
        <<z:=nil$
          for j:=1:dimension!* do
            <<z:=list('quotient,
                      list('diff,
                           list('times,
                                listarsq(listar(arg,j),i),
                                sfprod!*,
                                list('recip,
                                     getel list('sf,j))),
                           getel list('coordinats,j)),
                      sfprod!*) . z$
              for l:=1:dimension!* do
                z:=list('times,
                        listarsq(listar(arg,j),l),
                        getel list('christoffel,j,i,l)) . z>>$
          ans:=simp!*('plus.z) . ans>>
    else lprie list(" DIV of ",expn)$
  return list(vt-1,ans)
end$

put('div,'express,'divexpress)$

symbolic procedure laplexpress expn$
begin
  scalar arg,vt,ans$
  arg:=tensval car expn$
  vt:=tensrnk car expn$
  if vt=0 then
      <<for i:=1:dimension!* do
          ans:=list('diff,
                    list('times,
                         sfprod!*,
                         list('expt,
                              getel list('sf,i),
                              -2),
                         list('diff,
                              list('!*sq,arg,nil),
                              getel list('coordinats,i))),
                    getel list('coordinats,i)).ans$
        ans:=list(0,simp!* list('quotient,
                              'plus . ans,
                              sfprod!*))>>
    else if vt=1 then
        ans:=divexpress list gradexpress expn
    else lprie list(" LAPLACIAN of ",expn)$
  return ans
end$

put('lapl,'express,'laplexpress)$

symbolic procedure curlexpress expn$
begin
  scalar arg,vt,ans,ic1,ic2$
  arg:=tensval car expn$
  vt:=tensrnk car expn$
  if vt=1 then
      for each i in (if dimension!* = 3 then coordindx!*
                       else '(1) ) do
        <<ic1:=listar(cyclic1!*,i)$
          ic2:=listar(cyclic2!*,i)$
          ans:=simp!* list('times,
                         getel list('sf,i),
                         list('recip,sfprod!*),
                         list('difference,
                              list('diff,
                                   list('times,
                                        getel list('sf,ic2),
                                        listarsq(arg,ic2)),
                                   getel list('coordinats,ic1)),
                              list('diff,
                                   list('times,
                                        getel list('sf,ic1),
                                        listarsq(arg,ic1)),
                                   getel list('coordinats,ic2))))
                                                            . ans>>
    else lprie list(" CURL of ",expn)$
  return (if dimension!* = 3 then list(1,ans)
            else list(0,car ans))
end$

put('curl,'express,'curlexpress)$
flag('(cons grad div lapl curl tens vect dyad dirdf !& !# !?)
     ,'tensfn)$

symbolic procedure exscalval u$
begin
  scalar fce,args$
  fce:=car u$
  args:=for each a in cdr u collect cddr express a$
  fce:=eval(get(fce,'express) . list mkquote args)$
  if car fce=0 then return cadr fce
    else typerr(u," is not scalar ")$
  return (nil . 1)
end$

algebraic$

infix #,?,&$
precedence .,**$
precedence #,.$
precedence ?,#$
precedence &,?$

symbolic flag('(cons !# !? div lapl curl dirdf),'full)$
symbolic for each a in '(cons !# !? div lapl curl dirdf)
   do put(a,'simpfn,'exscalval)$

symbolic procedure scalefactors transf$
begin
  scalar var$
  dimension!*:=car transf$
  transf:=cdr transf$
  if dimension!*=2 then
      <<var:=cddr transf$
        transf:=list( car transf,cadr transf)>>
    else if dimension!*:=3 then
      <<var:=cdddr transf$
        transf:=list(car transf,cadr transf,caddr transf)>>
    else lprie list(" Can't handle dimension = ",dimension!*)$
  if dimension!*=length var then t
    else lprie list(" Transformation ",transf,var)$
  for i:=1:dimension!* do
    setel(list('coordinats,i),listar(var,i))$
  for row:=1:dimension!* do
    for col:=1:dimension!* do
      setel(list('jacobian,row,col),
            aeval list('df,
                       listar(transf,col),
                       getel list('coordinats ,row)))$
  updatedimen()$
  rscale()
end$

deflist('((scalefactors rlis)),'stat)$

flag('(remd),'eval);
remd 'jacobian; remprop('jacobian,'opfn);   % For bootstrapping.

array jacobian(3,3),coordinats (3),sf(3),christoffel(3,3,3)$

procedure rscale()$
begin
  sfprod!*:=1$
  nscal!*:=nscal!* + 1$
  for row:=1:dimension!* do
    <<for col:=1:(row-1) do
        <<sf(row):=gcov(row,col)$
          if sf(row)=0 then nil
            else write " WARNING: Coordinate system is nonorthogonal"
                       ," unless following simplifies to zero: ",
                       sf(row) >>$
      sf(row):=sqrt gcov(row,row)$
      sfprod!*:=sfprod!* *sf(row)>>$
  on nero$
  for i:=1:dimension!* do for j:=1:dimension!* do
    for k:=1:dimension!* do begin christoffel(i,j,k):=
      ((if i=j then df(sf(j),coordinats (k)) else 0)
      -(if i=k then df(sf(k),coordinats (j)) else 0))
      /(sf(j)*sf(k))$
      if wrchri(a)=0 then write christoffel(i,j,k):=
          christoffel(i,j,k)   end$
  off nero
end$

procedure gcov(j,k)$
for l:=1:dimension!* sum jacobian(j,l)*jacobian(k,l)$

symbolic$

symbolic procedure simpwrchri u$
if !*wrchri then nil . 1
  else 1 . 1$

put('wrchri,'simpfn,'simpwrchri)$

symbolic procedure rmat$
'dyad . cdr matstat()$

symbolic procedure formdyad(u,v,m)$
'list . mkquote 'dyad . cddr formmat(u,v,m)$

put('dyad,'stat,'rmat)$
put('dyad,'formfn,'formdyad)$

symbolic procedure dirdfexpress expn$
begin
  scalar arg,vt,direc,ans,z,dj,di,argj,sfj,sfi,cooj$
  arg:=cadr expn$
  vt:=tensrnk arg$
  direc:=car expn$
  if not (tensrnk direc=1) then return
      lprie list (" Direction in DIRDF is not a vector ",expn)$
  if vt=0 then return inprodexpress list (direc,
                         gradexpress list arg)$
  arg:=tensval arg$
  direc:=tensval direc$
  if not(vt=1) then return
      lprie list (" Argument of DIRDF is dyadic ",expn)$
  for each i in coordindx!* do
    <<z:=nil$
      di:=listarsq(direc,i)$
      sfi:=getel list('sf,i)$
      for j:=1:dimension!* do
        <<dj:=listarsq(direc,j)$
          argj:=listarsq(arg,j)$
          sfj:=getel list('sf,j)$
          cooj:=getel list('coordinats,j)$
          z:=list('times,
                  dj,
                  list('recip,sfj),
                  list('diff,
                       listarsq(arg,i),
                       cooj)) . z$
          z:=list('times,
                  di,
                  argj,
                  list('recip,sfi),
                  list('recip,sfj),
                  list('df,sfi,cooj)) . z$
          z:=list('minus,
                  list('times,
                       dj,
                       argj,
                       list('recip,sfi),
                       list('recip,sfj),
                       list('df,
                            sfj,
                            getel list('coordinats,i)))) . z>>$
      z:='plus . z$
      z:=simp!* z$
      ans:=z . ans >>$
  return list(1,ans)
end$

put('dirdf,'express,'dirdfexpress)$

symbolic procedure dfexpress expn$
begin
  scalar arg,vt,rest$
  arg:=tensval car expn$
  vt:=tensrnk car expn$
  rest:=cdr expn$
  rest:=for each a in rest collect
          if tensrnk a=0 then
              if atom tensval a then tensval a
                else if cdr tensval a=1 and numberp car tensval a
                       then car tensval a
                else !*q2k tensval a
            else lprie list(" Bad arg of DF ",expn)$
  if vt=0 then return list(0,simpdf(list('!*sq,arg,t) . rest))
    else if vt=1 then return list(1,for each a in arg collect
                                       simpdf(list('!*sq,a,t) . rest) )
    else if vt=2 then return list(2,for each a in arg collect
                                      for each b in a collect
                                        simpdf(list('!*sq,b,t) . rest) )
    else lprie list(" Bad tensor in DF ",expn)
end$

put('df,'express,'dfexpress)$

symbolic procedure diffexpress expn$
begin
  scalar arg,vt,rest$
  arg:=tensval car expn$
  vt:=tensrnk car expn$
  rest:=cdr expn$
  rest:=for each a in rest collect
             if tensrnk a=0 then
                 if atom tensval a then tensval a
                   else if cdr tensval a=1 and numberp car tensval a
                    then car tensval a
                   else !*q2k tensval a
               else lprie list(" Bad arg of DIFF ",expn)$
  if vt=0 then return list(0,simp('diff . (prepsq arg . rest)))
    else if vt=1 then return list(1,for each a in arg collect
                                    simp('diff . (prepsq a . rest)) )
    else if vt=2 then return list(2,for each a in arg collect
                                      for each b in a collect
                                        simp('diff . (prepsq b . rest)))
    else lprie list(" Bad tensor in DIFF ",expn)
end$

put('diff,'express,'diffexpress)$

remprop('diff,'number!-of!-args);  % Until we understand what's up.

algebraic$

scalefactors 3,x,y,z,x,y,z$

endmodule;

end;
REDUCE Historical
[ REDUCE Sourceforge Project | Historical SVN Repository | GitHub Mirror | SourceHut Mirror | NotABug Mirror | Chisel Mirror | Chisel RSS ]