reduce-historical/r34/lib/fide1.red

COMMENT                     ***** NOTE *****;

% For fastloading of this package please follow these steps:
%    faslout fide1;
%    in "fide1.red";
%    faslend;
%    load "fide1";   % used to need gentran, but no longer
%    faslout fide;
%    in "fide.red"$
%    faslend;

% For running this package the matrix, gentran, fide1 and fide packages
% have to be loaded.  However, loading fide is enough to make this
% happen, since matrix, gentran and fide1 are automatically loaded.
% Version 1.1 of the FIDE package is the result of porting the FIDE
% package version 1.0.0 to REDUCE 3.4. Some addition presented in the
% LINBAND module.

off echo$
load!-package 'matrix;

%***********************************************************************
%*****                                                             *****
%*****     P a c k a g e      F I D E    Ver. 1.1   17/05/1991     *****
%*****                                                             *****
%***********************************************************************
%**                                                                   **
%**FInite difference method for partial Differential Equation systems **
%**                                                                   **
%***********************************************************************
%**  (C) 1991, Richard Liska                                          **
%**      Faculty of Nuclear Science and Physical Engineering          **
%**      Technical University of Prague                               **
%**      Brehova 7                                                    **
%**      115 19 Prague 1                                              **
%**      Czechoslovakia                                               **
%**      Email: Richard Liska <tjerl@cspuni12.bitnet>                 **
%**                                                                   **
%**  This package can be distributed through REDUCE Network Library.  **
%***********************************************************************

% The FIDE package consists of the following modules:
%
%  EXPRES  for transforming PDES into any orthogonal coordinate system.
%  IIMET   for discretization of PDES by integro-interpolation method.
%  APPROX  for determining the order of approximation of difference
%          scheme
%  CHARPOL for calculation of amplification matrix and characteristic
%          polynomial of difference scheme, which are needed in Fourier
%          stability analysis.
%  HURWP   for polynomial roots locating necessary in verifying the von
%          Neumann stability condition.
%  LINBAND for generating the block of FORTRAN code, which solves a
%          system of linear algebraic equations with band matrix
%          appearing quite often in difference schemes.
%
% These modules can also be used separately.


%***********************************************************************
%*****                                                             *****
%*****      M O D U L E      E X P R E S                           *****
%*****                                                             *****
%***********************************************************************

module expres$
% Author: R. Liska
% Program EXPRES, Version REDUCE 3.4    05/1991

symbolic$

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

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

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!*$

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 mapcar(cdr expn,function simp!*))
     else if car expn = 'dyad then return mktensor
       (2,testdim2 mapcar(cdr expn,function (lambda a$
                                   mapcar(a,function simp!*))))
     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:=mapcar(cdr expn,function (lambda a$ 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 . mapcar(lst,function(lambda a$
                       if car a=0 then list('!*sq,cdr a,nil)
                         else typerr(expn," well formed scalar "))) ))$
  lst:=mapcar(lst,function(lambda a$
               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 lstx,lsty,mtx,z,zz$
  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:=mapcar(y,function(lambda a$
                                              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,mapcar(x,function(lambda a$
                                    antisymsum(a,y))))
    else if rnx=1 and rny=2 then return
      list(dimension!* - 1,
           if dimension!*=3 then
               tp1 copy2(mapcar(tp1 copy2(y,t),
                         function(lambda a$ antisymsum(x,a))),t)
             else mapcar(tp1 copy2(y,t),
                         function(lambda a$ 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:=mapcar(cdr u,function(lambda a$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)$

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) >>$
      write 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:=mapcar(rest,function(lambda a$
                               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,mapcar(arg,function(
         lambda a$simpdf(list('!*sq,a,t) . rest))) )
    else if vt=2 then return list(2,mapcar(arg,function(
         lambda a$mapcar(a,function(
           lambda b$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:=mapcar(rest,function(lambda a$
                               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,mapcar(arg,function(
         lambda a$simp('diff . (prepsq a . rest))) ))
    else if vt=2 then return list(2,mapcar(arg,function(
         lambda a$mapcar(a,function(
           lambda b$simp('diff . (prepsq b . rest))) ))))
    else lprie list(" Bad tensor in DIFF ",expn)
end$

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

algebraic$

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

endmodule$

%***********************************************************************
%*****                                                             *****
%*****      M O D U L E      I I M E T                             *****
%*****                                                             *****
%***********************************************************************

module iimet$
% Author: R. Liska
% Program IIMET, Version REDUCE 3.4  05/1991$

symbolic$

global '(cursym!* !*val dimension!*)$
fluid '(!*exp alglist!*)$

symbolic procedure array u$
begin
  scalar msg,erfg$
  msg:=!*msg$
  !*msg:=nil$
  erfg:=erfg!*$
  erfg!*:=nil$
  arrayfn('symbolic,
          mapcar (u,function(lambda a$car a . sub1lis cdr a)))$
  erfg!*:=erfg$
  !*msg:=msg
end$

symbolic procedure sub1lis u$
   if null u then nil else ((car u - 1) . sub1lis cdr u)$

sfprod!*:=1$

global'(date!*!*)$

date!*!*:= "IIMET Ver 1.1  17/05/91"$

put('version,'stat,'rlis)$

put('diff,'simpfn,'simpiden)$

global '(coords!* icoords!* dvars!* grids!* given!* same!*
         difml!* iobjs!* !*twogrid !*eqfu !*fulleq !*centergrid)$

switch twogrid,eqfu,fulleq,centergrid$

!*twogrid:=t$      %  Given functions can be on both grids.
!*eqfu:=nil$       %  During pattern matching the given and
                   %  looked for functions are different.
!*fulleq:=t$       %  Optimalization is performed on both sides of PDE.
!*centergrid:=t$   %  Centers of grid cells are in points I
                   %  (otherwise in I+1/2).

icoords!*:='(i j k l m n i1 j1 k1 l1 m1 n1)$
%  Indices which are given implicit.

procedure coordfn$
% Stat procedure of the COORDINATES statement, which defines indexes
% of coordinates.
begin
  scalar cor,icor$
  flag('(into),'delim)$
  cor:=remcomma xread nil$
  remflag('(into),'delim)$
  if cursym!* eq 'into then
      icor:=remcomma xread nil
    else if cursym!* eq '!*semicol!* then
      icor:=icoords!*
    else return symerr('coordfn,t)$
  return list('putcor,
              mkquote cor,
              mkquote icor)
end$

put('coordinates,'stat,'coordfn)$
flag('(putcor),'nochange)$

procedure putcor(u,v)$
begin
  scalar j$
  j:=1$
  coords!*:=u$
  while u do
    <<if not idp car u then        msgpri
          (" Coordinate ",car u," must be identifier",nil,'hold)$
      if not(idp car v or fixp car v) then        msgpri
          (" Index ",car v," must be identifier or integer",nil,'hold)$
      put(car u,'index,list car v)$
      put(car v,'coord,list car u)$
      put(car u,'ngrid,j)$
      j:=j+1$
      put(car u,'simpfn,'simpiden)$
      u:=cdr u$
      v:=cdr v >>
end$

procedure tcar u$
if pairp u then car u
  else u$

procedure grid u$
%  Procedure definning the statement GRID.
eval list(get(car u,'grid),
          mkquote cdr u)$

put('grid,'stat,'rlis)$
put('uniform,'grid,'gridunif)$

procedure gridunif u$
flag(u,'uniform)$

procedure dependence u$
%  Procedure definning the statemnt DEPENDENCE.
begin
  scalar x,y,z,gg,l,te,yy,y1,yl$
  if null coords!* then rederr
       " Coordinates have not been defined yet"$
  gg:=explode '!*grid$
  l:=list(length coords!* + 1)$
a:x:=car u$
  y:=car x$
  if idp y then if not y memq dvars!* then dvars!*:=y . dvars!*
                  else nil
    else return msgpri(" Variable ",y," must be identifier",nil,
                       'hold)$
  z:=cdr x$
  x:=car z$
  if not numberp x then go to b$
  if x=1 then apply('vectors,list y)
    else if x=2 then apply('dyads,list y)
    else if x=0 then t
    else return errpri2(car u,'hold)$
  z:=cdr z$
b:yl:=nil$
  yy:=explode y$
  te:=aeval y$
  if eqcar(te,'tensor) then te:=caddr te
    else te:=nil$
  if te=1 then
      for i:=1:dimension!* do
        <<y1:=intern compress append(yy,explode i)$
          yl:=y1 . yl$
          setk1(list(y,i),y1,t)$
          x:=intern compress append(explode y1,gg)$
              %  The name of an array filled with grids of Y(I)
          put(y1,'grid,x)$
          eval list('array,mkquote list(x . l)) >>
    else if te=2 then
      for i:=1:dimension!* do
        for j:=1:dimension!* do
          <<y1:=intern compress append(yy,append(explode i,
                                                 explode j))$
            yl:=y1 .  yl$
            setk1(list(y,i,j),y1,t)$
            x:=intern compress append(explode y1,gg)$
              %  The name of an array filled with grids of Y(I)
            put(y1,'grid,x)$
            eval list('array,mkquote list(x . l)) >>
    else
      <<yl:=y . yl$
        x:=intern compress append(yy,gg)$
        put(y,'grid,x)$
        eval list('array,mkquote list(x . l)) >>$
  for each a in yl do put(a,'simpfn, 'simpiden)$
  put(y,'names,reverse yl)$
  if te member '(1 2) then
      <<put(y,'kkvalue,get(y,'kvalue))$
        remprop(y,'kvalue) >>$
  for each v in z do
    if v memq coords!* then
        for each w in yl do depend1(w,v,t)
      else msgpri(" Identifier ",v," is not coordinate",nil,'hold)$
  u:=cdr u$
  if u then go to a$
  return nil
end$

put('dependence,'stat,'rlis)$

procedure given u$
begin
  scalar x,xnam$
a:x:=car u$
  xnam:=get(x,'names)$
  if not idp x then msgpri
       (" Variable ",x," must be identifier",nil,'hold)
    else if xnam then given!* := union(xnam,given!*)
    else msgpri
       (" Identifier ",x," is not variable",nil,'hold)$
  u:=cdr u$
  if u then go to a$
  return nil
end$

put('given,'stat,'rlis)$

procedure cleargiven$
<<remflag(given!*,'noeq)$
  given!*:=nil >>$

put('cleargiven,'stat,'endstat)$
flag('(cleargiven),'eval)$

newtok'(( !. !. ) isgr)$
algebraic infix ..$
grids!* := '(one half)$

procedure trlis$
% Stat procedure of the statement ISGRID.
begin
  scalar x$
  put('!*,'newnam,'tims)$
  x:=rlis()$
  remprop('!*,'newnam)$
  return x
end$

procedure formtr(u,vars,mode)$
list('isgrid,mkquote cdr u)$

procedure isgrid u$
% Proceudre definning the statement ISGRID.
begin
  scalar x,y,z,z1,te,gd,lz,lz1$
a:x:=car u$
  y:=car x$
  x:=cdr x$
  if not(y memq dvars!*) then return msgpri
      (" Identifier ",y," is not variable",nil,'hold)$
  if null x then go to er$
  te:=aeval y$
  te:=if eqcar(te,'tensor) then caddr te
        else nil$
  if (te=1 and null atom x and indexp car x and gridp cdr x) or
     (te=2 and null atom x and indexp car x and null atom cdr x
      and indexp cadr x and gridp cddr x) or
     ((te=0 or null te) and null atom x and gridp x) then t
    else go to er$
  if te=1 then
      <<z:=car x$
        x:=cdr x$
        lz:=nil$
        if numberp z then lz:=z . lz
          else for i:=1:dimension!* do lz:=i . lz$
        for each a in lz do mapc(x,function(lambda b$
            setel(list(get(cadr assoc(list(y,a),
                                      get(y,'kkvalue)),
                           'grid),
                       car b),
                  cadr b . nil)))  >>
    else if te=2 then
      <<z:=car x$
        z1:=cadr x$
        x:=cddr x$
        lz:=lz1:=nil$
        if numberp z then lz:=z . lz
          else for i:=1:dimension!* do lz:=i . lz$
        if numberp z1 then lz1:=z1 . lz1
          else for i:=1:dimension!* do lz1:=i . lz1$
        for each a in lz do for each b in lz1 do
          mapc(x,function(lambda c$
            setel(list(get(cadr assoc(list(y,a,b),
                                      get(y,'kkvalue)),
                           'grid),
                       car c),
                  cadr c . nil))) >>
    else mapc(x,function(lambda c$
        setel(list(get(y,'grid),car c),cadr c . nil)))$
  u:=cdr u$
  if u then go to a$
  return nil$
er:errpri2(car u,'hold)
end$

put('isgrid,'stat,'trlis)$
put('isgrid,'formfn,'formtr)$

procedure indexp u$
u eq 'tims or (numberp u and 0<u and u<dimension!*)$

procedure gridp u$
null atom u and
  begin
    scalar x$
  a:x:=car u$
    if null atom x and car x eq 'isgr and null atom cdr x
       and cadr x memq coords!* and null atom cddr x and
       caddr x memq grids!* then rplaca(u,cdr x)
      else return nil$
    x:=car u$
    rplaca(x,get(car x,'ngrid))$
    u:=cdr u$
    if u then go to a$
    return t
  end$

procedure grideq u$
begin
  scalar x,y,z$
  x:=u$
a:y:=car x$
  if not car y memq dvars!* then return msgpri(
      "Identifier ",car y," is not variable",nil,'hold)$
  z:=cdr y$
b:if not atom car z and caar z eq 'isgr and cadar z memq coords!* and
      caddar z memq grids!* then put(car y,cadar z,caddar z)
    else return errpri2(u,'hold)$
  z:=cdr z$
  if z then go to b$
  x:=cdr x$
  if x then go to a$
  return nil
end$

put('grideq,'stat,'rlis)$

fluid '(coord unvars sunvars interpolp novars ncor nvar intp icor gvar
        hi hip1 hip2 him1 him2 lhs rhs lsun lun xxgrid resar)$

%  GVAR   -  grid on which is actual equation integrated
%  NVAR   -  number of actual equation in IIM21
%  NCOR   -  number of actual coordinate
%  UNVARS -  list of looked for functions
%  NOVARS -  list of functions controlled by the SAME statement
%  SUNVARS - list of looked for functions, which are not controlled by
%            SAME, and of given functions, which are not controlled by
%            SAME and which can be only on one grid (if OFF TWOGRID)
%            in case ON TWOGRID given functions are not in SUNVARS
%            distinguishing of given functions defined in ISGRID is done
%            in procedures TWOGRIDP and GETGRID
%  LSUN    - length of  SUNVARS-1
%  INTERPOLP-flag for MATCHLINEAR (for OFF EXP), has value T for calling
%            from INTERPOL and NIL for calling from SIMPINTT - for
%            NGETVAR in SIMPINTT
%  RESAR  -  the name of an array which will be filled by the result

procedure zero u$
nil . 1$

procedure ngetvar(u,v)$
if atom u then get(u,v)
  else if atom car u then get(car u,v)
  else nil$

procedure ngrid u$
if u eq 'one then 'none
  else if u eq 'half then 'nhalf
  else nil$

procedure gnot u$
% Procedure inverts denotation of half-integer and integer grid
if u='one then 'half
  else if u='half then 'one
  else nil$

procedure same u$
begin
  scalar x,y,z,bo$
  if null same!* then <<same!*:=list u$
                        return nil >>$
  x:=same!*$
a:y:=car x$
  z:=u$
  bo:=nil$
  while z and not bo do
    <<if car z memq y then bo:=t$
      z:=cdr z >>$
  if bo then go to b$
  x:=cdr x$
  if x then go to a$
  same!*:= u . same!*$
  return nil$
b:rplaca(x,union(y,u))$
  return nil
end$

put('same,'stat,'rlis)$

procedure clearsame$
same!*:=nil$

put('clearsame,'stat,'endstat)$
flag('(clearsame),'eval)$

procedure mksame$
begin
  scalar x,y,z,yy,bo$
  x:=expndsame()$
a:y:=car x$
  yy:=y$
  while yy and not(car yy memq unvars) do yy:=cdr yy$
  if null yy then
       <<msgpri
      (" Same declaration ",nil,list(y,
       " doesn't contain unknown variable"),nil,'hold)$
         go to b >>$
  if y neq yy then
      <<z:=car y$         % The looked for function on the first place
        rplaca(y,car yy)$
        rplaca(yy,z) >>$
  z:=car y$
  yy:=cdr y$
  put(z,'sames,yy)$
  novars:=union(novars,yy)$
  for each a in cdr y do
  % Testing if A has appeared in the statement DEPENDENCE
    if not get(a,'grid) then msgpri
        (" Identifier ",a," is not variable",nil,'hold)
      else put(a,'same,z)$
  for i:=1:length coords!* do
    <<yy:=y$
      bo:=nil$
      while yy and not bo do     % Test on ISGRID
        <<bo:=getel list(get(car yy,'grid),i)$
          yy:=cdr yy >>$
      if bo then filgrid(y,bo,i) >>$
b:x:=cdr x$
  if x then go to a$
  sunvars:=setdiff(unvars,novars)$
  return unvars
end$

procedure filgrid(y,bo,i)$
% Filling up after finding ISGRID according to SAME
begin
  scalar yy,bg$
  yy:=y$
  while yy do
    <<bg:=getel list(get(car yy,'grid),i)$
      if null bg then setel(list(get(car yy,'grid),i),bo)
        else if bg eq bo then t
        else msgpri
          (" Same declaration ",nil,list(y,
           " doesn't correspond to isgrid declarations"),nil,'hold)$
       yy:=cdr yy>>
end$

procedure expndsame$
% Extending SAME!* by new identifiers for vectors and tensors
begin
  scalar x,y,sam$
  x:=same!*$
a:y:=mapcan(car x,function(lambda a$copy1 get(a,'names)))$
  sam:=y . sam$
  x:=cdr x$
  if x then go to a$
  return sam
end$

procedure copy1 u$
if null u then nil
  else if atom u then u
  else car u . copy1 cdr u$

procedure nrsame$
% Changing the numbering of variables according to SAME
for each a in sunvars do
  begin
    scalar x,nx$
    x:=get(a,'sames)$
    if x then
        <<nx:=get(a,'nrvar)$
          for each b in x do put(b,'nrvar,nx) >>$
    return nil
  end$

procedure iim u$
% Procedure defines the statement IIM
begin
  scalar xx,xxx,be,beb1,val,twogr$
  iim1 u$
  iobjs!*:=append(unvars,append(coords!*,given!*))$
  val:=!*val$
  !*val:=nil$
  novars:=sunvars:=nil$
  if same!* then mksame() else sunvars:=unvars$
  twogr:=!*twogrid$
  xxx:=setdiff(given!*,novars)$
  if !*twogrid then
      if null xxx then !*twogrid:=nil
        else flag(xxx,'twogrid)
    else sunvars:=union(sunvars,xxx)$
  flag(given!*,'noeq)$
  xxx:=0$                %  Numbering of variables and equation
  for each a in sunvars do
    <<put(a,'nrvar,xxx)$
      xxx:=xxx+1 >>$
  if same!* then nrsame()$
  xxx:=0$
  for each a in unvars do
    <<put(a,'nreq,xxx)$
      xxx:=xxx+1 >>$
  lun:=length unvars-1$
  lsun:=length sunvars-1$
  eval list('array,mkquote list('!*f2 . add1lis list(lun,lsun,1)))$
  xxx:=coords!*$
d:coord:=car xxx$
  icor:=tcar get(coord,'index)$
  difml!*:=nil$
  for i:=0:10 do difml!*:=append(difml!*,
       for each a in getel list('difm!*,i) collect
           if (xx:=atsoc(coord,cdr a)) then car a . cdr xx
             else if (xx:=atsoc('all,cdr a)) then car a . cdr xx
             else nil )$
  difml!*:=mapcon(difml!*,function(lambda a$if null car a then nil
                                              else list car a ))$
  if !*twogrid then difml!*:=
      for each a in difml!* collect
        if (xx:=caadr a) and (!*eqfu or memq(caar xx,'(f g))) then
            (car a . extdif(cdr a,nil))
          else a$
  be:=iim2 ()$
  iim21 be$
  if car be then beb1:=iim22 be
    else beb1:=list(car be,cadr be,car be)$
  if not fixp intp then msgpri(" INTP after heuristic search ",
            nil,list("is not a number, INTP=",intp),nil,nil)$
  if not(intp=0) then iim3 beb1$
  iim4 ()$
  xxx:=cdr xxx$
  if xxx then go to d$
  iim5 ()$
  for each a in '(rtype avalue dimension) do remprop('!*f2,a)$
  !*val:=val$ !*twogrid:=twogr$
  return nil
end$

procedure extdif(x,lg)$
% Performs corrections of diff. schemes for given functions on
% two grids - everytime chooses the scheme with minimal penalty.
% LG - list of all terms from (U V W F G), which has been in X
%      already changed and choosen.
begin
  scalar olds,news,y,gy,xx,lgrid,gg,g1$
  lgrid:=get('difm!*,'grids)$
  gy:=caar x$
  gg:=gy$
  for each a in lg do gg:=delete(atsoc(a,gg),gg)$
  if gg then gg:=caar gg
    else return x$
  x:=mapcar(x,function(lambda a$a))$
a:xx:=x$
  y:=car x$
  gy:=car y$
  g1:=atsoc(gg,gy)$
  gy:=(gg . gnot cdr g1) . delete(g1,gy)$
  gy:=acmemb(gy,lgrid)$
  while cdr xx and not eqcar(cadr xx,gy) do xx:=cdr xx$
  if cdr xx then
      <<olds:=y . (cadr xx . olds)$
        y:=if cadr y > cadadr xx
              or (cadr y=cadadr xx
                     and sublength caddr y > sublength car cddadr xx)
             then cadr xx
             else y$
        rplacd(xx,cddr xx) >>
    else olds:=y . olds$
  gy:=car y$
  g1:=atsoc(gg,gy)$
  gy:=delete(g1,gy)$
  if null gy then t
    else if (xx:=acmemb(gy,lgrid)) then gy:=xx
    else nconc(lgrid, list gy)$
  y:=gy . cdr y$
  news:=y . news$
  x:=cdr x$
  if x then go to a$
  if(xx:=caar news) and (!*eqfu or memq(caar xx,'(f g))) then
      <<olds:=extdif(olds,gg . lg)$
        news:=extdif(news,lg) >>$
  return nconc(olds,news)
end$

procedure sublength u$
if atom u then 0
  else length u + sublengthca u$

procedure sublengthca u$
if null u then 0
  else sublength car u + sublengthca cdr u$

procedure iim1 u$
% Checks the syntax of the IIM statement, calculates scalar PDEs,
% vector and tensor equations are expanded to scalar components.
begin
  scalar x,xx,e,te,exp$
  terpri()$
  prin2t"*****************************"$
  prin2 "*****      Program      *****          "$
  prin2t date!*!*$
  prin2t"*****************************"$
  exp:=!*exp$
  !*exp:=t$
  rhs:=lhs:=unvars:=nil$
  if null coords!* then return lprie " Coordinates defined not yet"$
  if null dvars!* then return lprie " Variables defined not yet"$
  for each v in dvars!* do
    if eqcar((te:=aeval v),'tensor) and caddr te member '(1 2) then
        put(v,'kvalue,get(v,'kkvalue))$
  if atom u or not idp car u then return errpri2(u,'hold)$
  resar:=car u$
  u:=cdr u$
a:if atom u or atom cdr u then return errpri2(u,'hold)$
  x:=car u$
  if not idp x then return msgpri
      (" Parameter ",x," must be identifier",nil,'hold)
    else if not(x memq dvars!*) then return msgpri
      (" Identifier ",x," is not variable",nil,'hold)
    else if x memq unvars then return msgpri
      (" Variable ",x," has second appearance",nil,'hold)
    else if x memq given!* then return msgpri
      (" Variable ",x," cannot be declared given",nil,'hold)
    else unvars:=x . unvars$
  e:=cadr u$
  if not eqexpr e then return msgpri
      (" Parameter ",e," must be equation",nil,'hold)
    else e:=aeval list('times,
                       list('difference,cadr e,caddr e),
                       sfprod!*)$
  if atom e then return msgpri
      (" Equation ",e," must be P.D.E.",nil,'hold)$
  te:=aeval x$
  te:=if eqcar(te,'tensor) then caddr te
        else nil$
  if(te=1 and car e eq 'tensor and caddr e=1) or
    (te=2 and car e eq 'tensor and caddr e=2) or
    (null te and car e eq 'tensor and caddr e=0) then
      e:=cadddr e   %  Necessary to carrect after change in EXPRESS
    else if null te and car e eq '!*sq then e:=cadr e
    else return msgpri
      (" Tensor order of",x," does not correspond to order of ",e,
       'hold)$
  lhs:=e . lhs$
  u:=cddr u$
  if u then go to a$
  for each v in dvars!* do
    if eqcar((te:=aeval v),'tensor) and caddr te member '(1 2) then
        remprop(v,'kvalue)$
b:x:=car unvars$
  e:=car lhs$  %  Transformation of vectors and tensor into components
  te:=aeval x$
  te:=if eqcar(te,'tensor) then caddr te
        else nil$
  if te=1 then rhs:=append(e,rhs)
    else if te=2 then for each a in reverse e do
                        rhs:=append(a,rhs)
    else rhs:=e . rhs$
  xx:=append(get(x,'names),xx)$  %  Add the checking if given equation
  unvars:=cdr unvars$            %  solves given variable (tensor var.)
  lhs:=cdr lhs$
  if unvars then go to b$
  unvars:=xx$
  lhs:=rhs$
  put('diff,'simpfn,'zero)$      %  Splitting left and right hand side
  alglist!*:=nil . nil$          %  All derivatives go to left h.s.
  rhs:=mapcar(rhs,function resimp)$
  put('diff,'simpfn,'simpiden)$
  alglist!*:=nil . nil$
  x:=lhs$
  xx:=rhs$
  terpri()$
  prin2t "      Partial Differential Equations"$
  prin2t "      =============================="$
  terpri()$
c:rplaca(xx,negsq car xx)$
  rplaca(x,prepsq!* addsq(car x,car xx))$
  rplaca(xx,prepsq!* car xx)$
  maprin car x$
  prin2!* "    =    "$
  maprin car xx$
  terpri!* t$
  rplaca(x,prepsq simp car x)$
  rplaca(xx,prepsq simp car xx)$
  x:=cdr x$
  xx:=cdr xx$
  if x then go to c$
  terpri()$
  x:=length lhs-1$
  if x=0 then eval list('array, mkquote list(resar . add1lis list(1)))
    else eval list('array,mkquote list(resar . add1lis list(x,1)))$
  !*exp:=exp$
  return nil
end$

procedure iim2$
% Defines the steps of the grid, splits variables to free and predefined
% grid in actual coordinate.
begin
  scalar b,e,xx,dihalf,dione,dihalfc$
  e:=append(explode 'h,explode coord)$
  e:=intern compress e$
  if flagp(coord,'uniform) then hi:=hip1:=him1:=him2:=e
    else <<put(e,'simpfn,'simpiden)$
           xx:=tcar get(coord,'index)$
           him1:=list(e,list('difference,xx,1))$
           him2:=list(e,list('difference,xx,2))$
           hi:=list(e,xx)$
           hip2:=list(e,list('plus,xx,2))$
           hip1:=list(e,list('plus,xx,1)) >>$
  dihalf:=list(
        'di . list('quotient,
                   list('plus,him1,hi),
                   2),
        'dim1 . him1,
        'dip1 . hi,
        'dim2 . list('quotient,
                     list('plus,him2,him1),
                     2),
        'dip2 . list('quotient,
                     list('plus,hi,hip1),
                     2))$
  dihalfc:=list(
        'di . list('quotient,
                   list('plus,hip1,hi),
                   2),
        'dim1 . hi,
        'dip1 . hip1,
        'dim2 . list('quotient,
                     list('plus,hi,him1),
                     2),
        'dip2 . list('quotient,
                     list('plus,hip2,hip1),
                     2))$
  dione:=list(
        'di . hi,
        'dim1 . list('quotient,
                     list('plus,him1,hi),
                     2),
        'dip1 . list('quotient,
                     list('plus,hi,hip1),
                     2),
        'dim2 . him1,
        'dip2 . hip1)$
  put('steps,'one,
      ('i . icor) . (if !*centergrid then dione else dihalf))$
  put('steps,'half,
      ('i . list('plus,
                 icor,
                 '(quotient 1 2))) . (if !*centergrid then dihalfc
                                        else dione))$
  ncor:=get(coord,'ngrid)$     % Number of the COODR coordinate
  e:=nil$
  for each a in sunvars do     % Splitting of variables with predefined
                               % grid.
    if (xx:=getel list(get(a,'grid),ncor))
            and car xx then e:=a . e$
  b:=setdiff(sunvars,e)$
  return list(b,e)
end$

procedure filfree(var,vgrid,freelst,pgr,peq)$
begin
  scalar x,nx,grn,nv,ng,ngrn,g1,g2,saml,bsam,asam,egrid$
  x:=ngetvar (var,'nrvar)$
c:put(var,pgr,vgrid)$
  egrid:=vgrid$
  if flagp(var,'noeq) then go to d$
  nx:=ngetvar (var,'nreq)$
% calulating in which point will be the euation for VAR integrated
  if egrid:=get(var,coord) then go to a
    else egrid:=vgrid$
  put('f2val,'free,'f2vzero)$
  if (g1:=f2eval(nx,x,0)) > (g2:=f2eval(nx,x,1)) then
      egrid:=gnot vgrid$
  if not g1=g2 then go to a$
  put('f2val,'free,'f2vmin)$
  if(g1:=f2eval(nx,x,0)) > (g2:=f2eval(nx,x,1)) then
      egrid:=gnot vgrid$
  if not g1=g2 then go to a$
  put('f2val,'free,'f2vmax)$
  if (g1:=f2eval(nx,x,0)) > (g2:=f2eval(nx,x,1)) then
      egrid:=gnot vgrid$
a:put(var,peq,egrid)$
% Penalties for free variables in the equation for VAR
  grn:=gnot egrid$
  ng:=ngrid egrid$
  ngrn:=ngrid grn$
  for each a in freelst do
    <<nv:=ngetvar (a,'nrvar)$
      asam:=a . get(a,'sames)$
      for each aa in asam do put(aa,pgr,egrid)$
      put(a,ng,cfplus2(get(a,ng),getel list('!*f2,nx,nv,0)))$
      for each aa in asam do put(aa,pgr,grn)$
      put(a,ngrn,cfplus2(get(a,ngrn),getel list('!*f2,nx,nv,1)))$
      for each aa in asam do remprop(aa,pgr) >>$
  if bsam then go to d$
  saml:=get(var,'sames)$
  bsam:=t$
d:if null saml then go to b$
  var:=car saml$
  saml:=cdr saml$
  go to c$
b:return egrid
end$

procedure f2eval(i,j,k)$
eval getel list('!*f2,i,j,k)$

procedure f2plus(i,j,k,l)$
% Procedure fills F2(I,J,K) with the number F2(I,J,K)+L$
begin
  scalar ma,x,y$
  if pairp l then
      if length car l=2 and cadr l=caddr l then l:=cadr l
        else if length l=3 and cadr l=caddr l and cadr l=cadddr l and
                cadr l=car cddddr l then l:=cadr l$
  ma:=list('!*f2,i,j,k)$
  x:=getel ma$
  if numberp l then
      if numberp x then setel(ma,x+l)
        else rplaca(x,car x+l)
    else
      if numberp x then setel(ma,list(x,l))
        else if (y:=assoc(car l,cdr x)) then
                 <<rplaca(cdr y,cadr y+cadr l)$
                   rplaca(cddr y,caddr y+caddr l)$
                   if cddar l then
                       <<rplaca(cdddr y,cadddr y + cadddr l)$
                         rplaca(cddddr y,car cddddr y+car cddddr l)>>>>
               else rplacd(x,(l . cdr x))
end$

procedure f2var u$
% Forms the elements of array !*F2 into the form
%               (FPLUS <CISLO> {(F2VAL U V N1 N2)})
if numberp u then u
  else ('fplus .
            car u . mapcar(cdr u,function(lambda a$
                      list('f2val,car a,cdr a))))$

fexpr procedure f2val x$
% Evaluates the expression (F2VAL ...)
begin
  scalar us,ns,u,v,w,n1,n2,n3,n4,gu,gv,gw$
  us:=car x$
  ns:=cadr x$
  u:=car us$
  v:=cadr us$
  n1:=car ns$
  n2:= cadr ns$
  gu:=get(u,xxgrid)$
  gv:=get(v,xxgrid)$
  if cddr us then
      <<w:=caddr us$
        gw:=get(w,xxgrid)$
        n3:=caddr ns$
        n4:=cadddr ns>>$
  return
    if w then
        if gu and gv and gw then
            if gu eq gv and gu eq gw then n1
              else if gu eq gv then n2
              else if gu eq gw then n3
              else if gv eq gw then n4
              else apply(get('f2val,'free),list(us,ns))
          else if gu and gv then
            if gu eq gv then aplf2val(u,w,n1,n2)
              else aplf2val(u,w,n3,n4)
          else if gu and gw then
            if gu eq gw then aplf2val(u,v,n1,n3)
              else aplf2val(u,v,n2,n4)
          else if gv and gw then
            if gv eq gw then aplf2val(u,v,n1,n4)
              else aplf2val(u,v,n2,n3)
          else apply(get('f2val,'free),list(us,ns))
      else if gu and gv then
        if gu eq gv then n1
          else n2
      else apply(get('f2val,'free),list(us,ns))
end$

procedure aplf2val(u,v,n1,n2)$
apply(get('f2val,'free),list(list(u,v),list(n1,n2)))$

fexpr procedure fplus u$
% Evaluates the expression (FPLUS ...)
begin
  scalar x,y,z$
% U:=CAR U$
  y:=car u$
a:u:=cdr u$
  z:=eval car u$
  if numberp z then y:=y+z
    else x:=z . x$
  if cdr u then go to a$
  return
    if x then ('fplus . y . x)
      else y
end$

procedure cfplus2(u,v)$
% Adds the expressions of type (FPLUS ...).
% Destroys U, does not destroy V.
begin
  scalar f2v$
  f2v:=get('f2val,'free)$
  put('f2val,'free,'f2vunchange)$
  if not fixp u then u:=eval u$
  if not fixp v then v:=eval v$
  put('f2val,'free,f2v)$
  return
    if fixp u then
        if fixp v then (u + v)
          else ('fplus . (cadr v+u) . cddr v)
      else if numberp v then <<rplaca(cdr u,cadr u+v)$ u>>
      else <<rplaca(cdr u,cadr u+cadr v)$
             nconc(u,cddr v) >>
end$

procedure f2vunchange(us,ns)$
list('f2val,us,ns)$

procedure f2vzero(us,ns)$
0$

procedure f2vplus(us,ns)$
eval('plus . ns)$

procedure f2vmax(us,ns)$
eval('max . ns)$

procedure f2vmin(us,ns)$
eval('min . ns)$

put('f2val,'fselect,'f2vplus)$
put('f2val,'fgrid,'f2vmin)$

procedure iim21 u$
% Fills the array !*F2 according to the system of PDE and penalties
% given.
% Fills the properties NONE,NHALF (FONE,FHALF) of free variables.
% According to predefined variables filles the properties XGRID and EQ
% of predefined variables.
begin
  scalar b,e,lh,lhe,xx,rh,bdef$
  b:=car u$            % Free vars
  e:=cadr u$           % Predefined vars
  for i:=0:lun do
    for j:=0:lsun do   % Filling the array F2
      <<setel(list('!*f2,i,j,0),0)$
        setel(list('!*f2,i,j,1),0) >>$
  nvar:=0$             % Number of actual variable
  lh:=lhs$
  rh:=rhs$
  interpolp:=nil$
  put('intt,'simpfn,'simpiden)$
a:lhe:=car lh$         % Actual equation
  lhe:=formlnr list('intt,lhe,coord)$
  rplaca(lh,lhe)$
  bdef:=t$
  for each var in sunvars do
    if get(var,coord) then t
      else bdef:=nil$
  if null b and bdef then go to c$
      % If there are no free variables it is not necessary to fill
      % the array F2 - no optimalization is necessary -> To use this
      % statement, we have to test if we know in which point (over
      % which interval) will all equation be integrated (discretized).
  put('intt,'simpfn,'simpintt)$
  alglist!*:=nil . nil$
  simp lhe$
  put('intt,'simpfn,'simpiden)$
  if !*fulleq then     % Optimalizatioon is performed for both sides of
      <<lhe:=car rh$   % equations, otherwise only for left hand side.
        lhe:=formlnr list('intt,lhe,coord)$
        put('intt,'simpfn,'simpintt)$
        simp lhe$
        alglist!*:=nil . nil$
        put('intt,'simpfn,'simpiden)$
        rh:=cdr rh >>$
c:nvar:=nvar+1$
  lh:=cdr lh$
  if lh then go to a$
  for i:=0:lun do
    for j:=0:lsun do
      for k:=0:1 do
        setel(list('!*f2,i,j,k),f2var getel list('!*f2,i,j,k))$
  xxgrid:='xgrid$
  for each a in b do
    <<put(a,'none,0)$
      put(a,'nhalf,0) >>$
  for each a in e do    % Predefined variables
    filfree(a,car getel list(get(a,'grid),ncor),b,'xgrid,'eq)$
  % Predefined penalties
  intp:=0$
  for each a in e do
    if a memq unvars then
        <<xx:=ngetvar(a,'nreq)$
          for each c in e do
            if get(c,'xgrid) eq get(a,'eq) then intpfplus(xx,c,0)
              else intpfplus(xx,c,1) >>$
  for each a in b do
    <<put(a,'fone,get(a,'none))$
      put(a,'fhalf,get(a,'nhalf)) >>$
  return nil
end$

procedure iim22 u$
begin
  scalar b,e,bb,b1,b2,x,xx,x1,nv,g1,g2$
  b:=car u$
  e:=cadr u$
  bb:=b$             % Heuristic determination of grids for
                     % variables from B
f:x:=car bb$         % Chose the next variable X
  put('f2val,'free,get('f2val,'fselect))$
  xx:=abs(eval get(x,'none)-eval get(x,'nhalf))$
  b2:=cdr bb$
  while b2 do
    if xx<(x1:=abs(eval get(car b2,'none)-eval get(car b2,'nhalf)))
        then <<x:=car b2$
               xx:=x1$
               b2:=cdr b2 >>
      else b2:=cdr b2$
  b1:=x . b1$       %  List of variables subsequently choosen from B
  bb:=delete(x,bb)$ %  List of variables remaining in B
  put('f2val,'free,get('f2val,'fgrid))$
  put(x,'xgrid,'one)$
  g1:=eval get(x,'none)$
  put(x,'xgrid,'half)$
  g2:=eval get(x,'nhalf)$
  if g1>g2 then xx:='half
    else xx:='one$
  filfree(x,xx,bb,'xgrid,'eq)$
  intpgplus(x,xx)$
  for each ax in (x . get(x,'sames)) do
    if ax memq unvars then
    <<x1:=ngetvar(ax,'nreq)$
      for each a in append(b1,e) do
        if get(a,'xgrid)=get(ax,'eq) then intpfplus(x1,a,0)
          else intpfplus(x1,a,1) >>$
  if bb then go to f$
  return list(b,e,b1)
end$

procedure intpfplus(nx1,a,n)$
intp:=cfplus2(intp,getel list('!*f2,nx1,ngetvar(a,'nrvar),n))$

procedure intpgplus(a,ga)$
intp:=cfplus2(intp,get(a,ngrid ga))$

procedure iim3 u$
begin
  scalar b,e,b1,bb$
  prin2t" Backtracking needed in grid optimalization"$
  b:=car u$         % Free vars
  e:=cadr u$        % Predefined vars
  b1:=caddr u$
  for each a in b do        % Full search - bactracking
    <<put(a,'none,get(a,'fone))$
      put(a,'nhalf,get(a,'fhalf)) >>$
  xxgrid:='bxgrid$
  nbxgrid(e,'bxgrid,'beq,'xgrid,'eq)$
  put('f2val,'free,'f2vunchange)$
  varyback(b1,nil)$
  for each a in union(unvars,given!*) do
    <<remprop(a,'bxgrid)$
      remprop(a,'beq) >>$
  return nil
end$

procedure nbxgrid(u,ng,ne,og,oe)$
for each a in u do
  for each b in (a . get(a,'sames)) do
    <<put(b,ng,get(b,og))$
      put(b,ne,get(b,oe)) >>$

procedure varyback(bb,b1)$
% Performs full search of BB. B1 is B-BB. N is the number of
% interpolations performed up to now.
if null bb then
    begin
      scalar none,nhalf,n,eqg,i,j$
      n:=0$
      for each a in unvars do
        <<none:=nhalf:=0$
          i:=get(a,'nreq)$
          for each b in sunvars do
            <<j:=get(b,'nrvar)$
              none:=none + if get(b,'bxgrid) eq 'one then f2eval(i,j,0)
                             else f2eval(i,j,1)$
              nhalf:=nhalf + if get(b,'bxgrid) eq 'half
                               then f2eval(i,j,0)
                              else f2eval(i,j,1) >>$
          put(a,'beq,if (eqg:=get(a,coord)) then eqg
                       else if none<=nhalf then 'one
                       else 'half)$
          n:=n + if eqg then
                     if eqg eq 'one then none
                       else if eqg eq 'half then nhalf
                       else <<msgpri("VARYBACK ",eqg,nil,nil,'hold)$ 0>>
                   else if none<=nhalf then none
                   else nhalf >>$
      if n<intp then
          <<nbxgrid(b1,'xgrid,'eq,'bxgrid,'beq)$
            for each a in unvars do put(a,'eq,get(a,'beq))$
            intp:=n >>
    end
  else if intp=0 then t
  else <<varb(bb,b1,'one)$
         varb(bb,b1,'half) >>$

procedure varb(bb,b1,xx)$
% Subprocedure of VARYBACK procedure
% In BB are temporary free variables
% In B1 are temporary predefined variables (over BXGRID property)
begin
  scalar x$
  x:=car bb$
  for each a in (x . get(x,'sames)) do
    put(a,'bxgrid,xx)$
  return varyback(cdr bb,x . b1)
end$

procedure iim4$
begin
  scalar lh,rh,x,lhe,var$
  intp:=intp/6$
  prin2 intp$
  prin2 " interpolations are needed in "$
  prin2 coord$
  prin2t " coordinate"$
  for each a in unvars do
    <<prin2"  Equation for "$
      prin2 a$
      prin2" variable is integrated in "$
      prin2 get(a,'eq)$
      prin2t" grid point" >>$
  interpolp:=t$
  put('intt,'simpfn,'simpinterpol)$
  lh:=lhs$
  rh:=rhs$
  x:=unvars$
j:var:=car x$
  gvar:=get(var,'eq)$
  lhe:=car lh$
  alglist!*:=nil . nil$
  lhe:=prepsq simp lhe$
  rplaca(lh,lhe)$
  lhe:=car rh$
  lhe:=formlnr list('intt,lhe,coord)$
  lhe:=prepsq simp lhe$
  rplaca(rh,lhe)$
  x:=cdr x$
  lh:=cdr lh$
  rh:=cdr rh$
  if x then go to j$
  put('intt,'simpfn,'simpiden)$
  return lhs
end$

procedure iim5$
begin
  scalar lh,rh,val,nreq,ar$
  val:=!*val$
  !*val:=nil$
  for each a in union(union(unvars,sunvars),given!*) do
    <<remprop(a,'xgrid)$
      remprop(a,'eq)$
      remprop(a,'nreq)$
      remprop(a,'nrvar)$
      remprop(a,'sames)$
      remprop(a,'none)$
      remprop(a,'nhalf)$
      remprop(a,'fone)$
      remprop(a,'fhalf) >>$
  remflag(given!*,'twogrid)$
  remflag(given!*,'noeq)$
  terpri()$
  prin2t "         Equations after Discretization Using IIM :"$
  prin2t "         =========================================="$
  terpri()$
  lh:=lhs$
  rh:=rhs$
  nreq:=0$
k:rplaca(lh,prepsq!* simp!* car lh)$
  maprin car lh$
  prin2!* "   =   "$
  rplaca(rh,prepsq!* simp!* car rh)$
  maprin car rh$
  terpri!* t$
  terpri()$
  ar:=if null cdr lhs then list(resar,0) else list(resar,nreq,0)$
  setel(ar,car lh)$
  ar:=if null cdr lhs then list(resar,1) else list(resar,nreq,1)$
  setel(ar,car rh)$
  lh:=cdr lh$
  rh:=cdr rh$
  nreq:=nreq+1$
  if lh then go to k$
  !*val:=val$
  return nil
end$

put('iim,'stat,'rlis)$

array difm!*(10)$

procedure iscomposedof(x,objs,ops)$
if null x then nil
  else if atom x then if idp x then memq(x,objs)
                        else if fixp x then t
                        else nil
  else if idp car x and car x memq ops and cdr x then
            iscompos(cdr x,objs,ops)
  else nil$

procedure iscompos(x,objs,ops)$
if null x then t
  else if idp car x then car x memq objs and iscompos(cdr x,objs,
                                                       ops)
  else if numberp car x then iscompos(cdr x,objs,ops)
  else if atom car x then nil
  else if idp caar x then caar x memq ops and cdar x and
                iscompos(cdar x,objs,ops) and iscompos(cdr x,objs,ops)
  else nil$

global'(difconst!* diffuncs!*)$
difconst!* := '(i n di dim1 dip1 dim2 dip2)$
diffuncs!*:=nil$

procedure difconst u$
difconst!* := append(u,difconst!*)$

put('difconst,'stat,'rlis)$

procedure diffunc u$
<<diffuncs!*:=append(u,diffuncs!*)$
  for each a in u do put(a,'matcheq,'matchdfunc) >>$

put('diffunc,'stat,'rlis)$

procedure matchdfunc(u,v)$
begin
  scalar x,y$
  return
    if null u and null v then list t
      else if null u or null v then nil
      else if (x:=matcheq(car u,car v)) and (y:=matchdfunc(cdr u,cdr v))
        then union(x,y)
      else nil
end$

procedure difmatch u$
begin
  scalar l,gds,gdsf,pl,x,dx,y,z,coor$
  coor:=car u$
  if not atom coor then go to er$
  u:=cdr u$
  x:=car u$
  if not iscomposedof(x,'(u f x n v w g),
       append(diffuncs!*,'(diff times expt quotient recip)))then
      go to er$
  x:=prepsq simp x$
  l:=if atom x then 0 else length x$
  x:=x . nil$
  if null(y:=getel list('difm!*,l)) then setel(list('difm!*,l),list x)
    else if (z:=assoc(car x,y)) then x:=z
    else nconc(y,list x)$
  y:=cdr u$
a:gds:=nil$
  gdsf:=nil$
  if not eqexpr car y then go to b$
a1:if not(cadar y memq '(u v w f g) and caddar y memq grids!*) then
                                                        go to er$
  if cadar y memq '(f g) then gdsf:=(cadar y . caddar y) . gdsf
    else gds:=(cadar y . caddar y) . gds$
  y:=cdr y$
  if null y then go to er$
  if eqexpr car y then go to a1$
b:if not fixp car y then go to er$
  pl:=car y$
  y:=cdr y$
  if null y then go to er$
  if not iscomposedof(car y,difconst!*,append(diffuncs!*,'(u x f v w
     g plus minus difference times quotient recip expt)))then go to er$
  dx:=car y$
  y:=cdr y$
  gds:=nconc(gdsf,gds)$
  defdfmatch(x,gds,pl,list dx,coor)$
  if y then go to a$
  return nil$
er:errpri2(y,'hold)
end$

procedure defdfmatch(x,gds,pl,dx,coor)$
begin
  scalar y,z,yy$
  y:=get('difm!*,'grids)$
  if null y then put('difm!*,'grids,list gds)
    else if null gds then t
    else if (z:=acmemb(gds,y)) then gds:=z
    else nconc(y,list gds)$
  y:=cdr x$
  if y then
      if (yy:=atsoc(coor,y)) then
          if (z:=assoc(gds,cdr yy)) then
              <<rplacd(z,pl . dx)$
                msgpri(" Difference scheme for ",car x," redefined ",
                    nil,nil) >>
            else nconc(cdr yy,list(gds . (pl . dx)))
        else nconc(y,list(coor . list(gds . (pl . dx))))
    else rplacd(x,list(coor . list(gds . (pl . dx))))$
  return y
end$

deflist('((difmatch rlis) (cleardifmatch endstat)),'stat)$

procedure cleardifmatch$
for i:=0:10 do difm!*(i):=nil$

flag('(cleardifmatch),'eval)$

procedure acmemb(u,v)$
if null v then nil
  else if aceq(u,car v) then car v
  else acmemb(u,cdr v)$

procedure aceq(u,v)$
if null u then null v
  else if null v then nil
  else if car u member v then aceq(cdr u,delete(car u,v))
  else nil$

procedure matcheq(u,v)$
if null u or null v then nil
  else if numberp u then if u=v then list t else nil
  else if atom u then
    begin
      scalar x$
      x:=eval list(get(u,'matcheq),mkquote u,mkquote
                       (if atom v then list v else v))$
      return
        if x then x
          else if null !*exp and pairp v and car v memq
                  '(plus difference) then matchlinear(u,v)
          else nil
    end
  else if atom v then nil
  else if atom car u and car u eq car v then
      eval list(get(car u,'matcheq),mkquote cdr u,mkquote cdr v)
  else if null !*exp and car v memq'(plus difference)
       and car u eq 'diff then
      matchlinear(u,v)
  else nil$

algebraic operator matchplus$

fluid'(uu vv)$

procedure matchlinear(u,v)$
% Construction for OFF EXP and second and next coordinates
begin
  scalar x,uu,vv,alg$
  if not atom u then return
      if car u eq 'diff then matchlindf(u,v)
        else if car u eq 'times then matchlintimes(u,v)
        else nil$
  uu:=u$
  vv:='first$
  x:=formlnr list('matchplus,v,coord)$
  put('matchplus,'simpfn,'matchp)$
  alg:=alglist!*$
  alglist!*:=nil . nil$
  simp x$
  alglist!*:=alg$
  put('matchplus,'simpfn,'simpiden)$
  return
    if vv then list(u .  (if interpolp then v else vv))
      else nil
end$

procedure matchp y$
begin
  scalar x$
  if null vv then return(nil . 1)$
  x:=matcheq(uu,car y)$
  if null x then return
      begin
        vv:=nil$
        return(nil . 1)
      end$
  if vv eq 'first then return
      begin
        vv:=cdar x$
        return (nil . 1)
      end$
  if mainvareq(vv,cdar x) then return (nil . 1)$
  vv:=nil$
  return(nil . 1)
end$

unfluid '(uu vv)$

procedure mainvareq(x,y)$
if atom x then eq(x,y)
  else if car x memq iobjs!* then eq(car x,car y)
  else if car x memq '(diff expt) then
    (car y eq car x and mainvareq(cadr x,cadr y) and cddr x=cddr y)
  else nil$

procedure tlist x$
if atom x then list x
  else x$

procedure matchlindf(u,v)$
begin
  scalar x,y,b$
  x:=mapcar(cdr v,function fsamedf)$
  y:=cdar x$
  if null y then return nil$
  x:=for each a in x collect if y=cdr a then car a
                               else b:=t$
  if b then return nil$
  x:=(car v . x) . y$
  return matchdf(cdr u,x)
end$

procedure fsamedf u$
begin
  scalar x$
  return
    if atom u then nil . nil
      else if car u eq 'minus then
        <<x:=fsamedf cadr u$
          list('minus,car x) . cdr x >>
      else if car u eq 'diff then cadr u . cddr u
      else if car u eq 'times then
        begin
          scalar y,z$
          x:=cdr u$
        a:if null x or y=t then go to b$
          if numberp car x then z:=car x . z
            else if eqcar(car x,'diff) then
              <<if y then y:=t
                  else y:=cddar x$
                z:=cadar x . z >>
            else if depends(car x,coord) then y:=t
            else z:=car x . z$
          x:=cdr x$
          go to a$
        b:return if y=t then nil . nil
                   else ('times . z) . y
        end
    else nil . nil
end$

procedure matchlintimes(u,v)$
begin
  scalar x,y,z$
  y:=cadr v$
  if eqcar(y,'times) then y:=cdr y
    else if eqcar(y,'minus) and eqcar(cadr y,'times)
     then y:=  (-1) . cdadr y
    else return nil$
  x:=for each a in cdr v collect
       if eqcar(a,'times) then <<y:=intersect(y,cdr a)$
                                 a>>
         else if eqcar(a,'minus) and eqcar(cadr a,'times) then
           <<y:=intersect(y,cdadr a)$
             'times . ((-1) . cdadr a) >>
         else y:=nil$
  if null y then return nil$
  x:=for each a in x collect <<z:=setdiff(cdr a,y)$
                               if null cdr z then car z
                                 else 'times . z >>$
  x:=car v . x$
  return matchtimes(cdr u,x . y)
end$

procedure intersect(u,v)$
if null u or null v then nil
  else if member(car u,v) then car u . intersect(cdr u,v)
  else intersect(cdr u,v)$

procedure matchu(u,v)$
if car v memq unvars or (!*eqfu and car v memq given!*) then list(u . v)
  else if car v eq 'diff and not coord memq cddr v
                and matcheq(u,tlist cadr v)
           then list(u . (car v . (tlist cadr v . cddr v)))
  else if car v eq 'times then
% Product can be inside brackets or in DIFF
    begin
      scalar x,b1,vv$
      x:=for each a in cdr v collect a$   % To allow RPLACA
      vv:=car v . x$
      b1:=0$
      while x and b1<2 do
        <<if depends(car x,coord) then
              <<b1:=b1+1$
                if atom car x then rplaca(x,list car x) >>$
          x:=cdr x >>$
      return if b1=0 or b1>1 then nil
               else (u . vv)
    end
  else nil$

put('u,'matcheq,'matchu)$
put('v,'matcheq,'matchu)$
put('w,'matcheq,'matchu)$

procedure matchf(u,v)$
if car v memq given!* then list(u . v)
  else if car v eq 'diff and not coord memq cddr v
                and matchf(u,tlist cadr v)
           then list(u . (car v . (tlist cadr v . cddr v)))
  else nil$

put('f,'matcheq,'matchf)$
put('g,'matcheq,'matchf)$

procedure matchx(u,v)$
if car v eq coord then list t
  else nil$

put('x,'matcheq,'matchx)$

procedure matchtimes(u,v)$
begin
  scalar bool,bo,x,y,y1,asl$
  x:=u$
a:y:=t . v$
d:bool:=nil$
  while not bool and cdr y do
    <<y:=cdr y$
      bool:=matcheq(car x,car y)>>$
  if null bool then go to b$
  bo:=bool$
c: if not atom bo and not atom car bo then y1:=atsoc(caar bo,asl)
    else y1 := nil$
  if y1 and not y1=car bo then go to d$
  bo:=cdr bo$
  if bo then go to c$
  v:=delete(car y,v)$
  x:=cdr x$
  asl:=union(bool,asl)$
  if x then go to a$
  if v then return nil$
  return asl$
b:return if null cdr v and cdr x then
             if y:=matcheq('times . x,car v) then union(asl,y)
               else nil
           else nil
end$

put('times,'matcheq,'matchtimes)$

procedure matchexpt(u,v)$
if fixp cadr u then
    if cadr u=cadr v then matcheq(car u,car v)
      else nil
  else if cadr u='n then
      begin
        scalar x$
        x:=matcheq(car u,car v)$
        return if x then (('n . cadr v) . x)
                 else nil
      end
  else nil$

put('expt,'matcheq,'matchexpt)$

procedure matchquot(u,v)$
begin
  scalar man,mad$
  return if(man:=matcheq(car u,car v))
              and (mad:=matcheq(cadr u,cadr v)) then
             union(man,mad)
           else nil
end$

put('quotient,'matcheq,'matchquot)$

procedure matchrecip(u,v)$
matcheq(car u,car v)$

put('recip,'matcheq,'matchrecip)$

procedure matchdf(u,v)$
begin
  scalar x,asl,y$
  asl:=matcheq(car u,car v)$
  if null asl then return nil$
  y:=x:=append(cdr v,nil)$
  while x and car x neq coord do x:=cdr x$
  if null x then return nil
    else if null cddr u then
             if null cdr x or idp cadr x then go to df1
                 else return nil
      else if cdr x and caddr u=cadr x then t
      else return nil$
  rplacd(x,cddr x)$
df1:y:=delete(coord,y)$
  if null y or null interpolp then return asl
%  !!! Be aware !!! in mixed derivations of product
    else return list(car u . ('diff . (cdar asl . y)))
end$

put('diff,'matcheq,'matchdf)$

procedure finddifm u$
begin
  scalar x,v,asl,eqfu,b,bfntwo,bftwo1$
  eqfu:=!*eqfu$
  if eqfu then !*eqfu:=nil$
a:x:=t . difml!*$
  bftwo1:=bfntwo$
  bfntwo:=nil$
  if !*eqfu then b:=t$
  while cdr x and not asl do
    <<x:=cdr x$
      asl:=matcheq(caar x,u)$
% Test for given variables of type F, which have to be on one grid,
% if choosed substitution from DIFMATCH contains definition of F
% on grids.
      if asl then for each a in asl do
          if not atom a and car a memq'(f g) and (cadr a memq novars or
             (null !*twogrid and cadr a memq sunvars)) and
             null atsoc(car a,caadar x) then bfntwo:=t$
      if bfntwo then asl:=nil >>$
  !*eqfu:=eqfu$
  if null asl then
      if null b and eqfu then go to a
    else go to nm$
  return list(('x . coord) . delete(t,asl),cdar x)$
nm:if eqcar(u,'times) and null !*exp then
       <<!*exp:=t$
         alglist!*:=nil . nil$
         v:=prepsq simp u$
         v:=formlnr list('intt,v,coord)$
         !*exp:=nil$
         if null atom v and car v memq'(plus difference) then
             return ('special . simp v) >>$
  msgpri(" Matching of ",u," term not find ",nil,'hold)$
  if bfntwo or bftwo1 then
       lprie(" Variable of type F not defined on grids in DIFMATCH")$
  return nil
end$

procedure tdifpair x$
% From CDR ATSOC(.,ASL) makes an atom - free variable
<<if eqcar(x,'diff) then
      if eqcar(cadr x,'times) then
          <<x:=cdadr x$
            while x and not depends(car x,coord) do x:=cdr x$
            x:=car x>>
        else x:=cadr x$
  if pairp x then x:=car x$
  x >>$

procedure simpintt u$
begin
  scalar asl,agdsl,l,x,nv,y,x1,y1,nv1,n1,n2,nn1,nn2,
         x2,y2,nv2,n3,n4,n5,n6,lgrids,gds$
  u:=prepsq simp car u$
  if u=1
        then go to r$
  asl:=finddifm u$
  if null asl or eqcar(asl,'special) then go to r$
  agdsl:=cadr asl$              % List from DIFML!*
  asl:=car asl$                 % ASOC. list of assignments
  gds:=caar agdsl$
  l:=length gds$
  if l=0 then go to r$
a:y:=caar gds$
  x:=atsoc(y,asl)$
  if null x then go to er1$
  x:=tdifpair cdr x$
  if !*twogrid and flagp(x,'twogrid) then
      if l=1 then go to r
        else <<gds:=cdr gds$
               l:=l-1$
               go to a>>$
  nv:=ngetvar(x,'nrvar)$
  if l=1 then go to l1
    else go to l2$
l1:x:=assoc(list(y . 'one),agdsl)$
  if null x then go to er2$
  f2plus(nvar,nv,0,6*cadr x)$
  x:=assoc(list(y . 'half),agdsl)$
  if null x then go to er2$
  f2plus(nvar,nv,1,6*cadr x)$
  go to r$
l2:y1:=caadr gds$
  x1:=atsoc(y1,asl)$
  if null x1 then go to er1$
  x1:=tdifpair cdr x1$
  if !*twogrid and flagp(x1,'twogrid) then
      if l=2 then go to l1
        else <<gds:=cdr gds$
               l:=l-1$
               go to l2 >>$
  nv1:=ngetvar(x1,'nrvar)$
  lgrids:=get('difm!*,'grids)$
  if l=3 then go to l3
    else if l>3  then go to er$
l20:n1:=atsoc(acmemb(list(y . 'one,y1 . 'one),lgrids),agdsl)$
  n2:=atsoc(acmemb(list(y . 'one,y1 . 'half),lgrids),agdsl)$
  nn1:=atsoc(acmemb(list(y . 'half,y1 . 'half),lgrids),agdsl)$
  nn2:=atsoc(acmemb(list(y . 'half,y1 . 'one),lgrids),agdsl)$
  if n1 and n2 and nn1 and nn2 then t
    else go to er2$
  n1:=cadr n1$
  n2:=cadr n2$
  nn1:=cadr nn1$
  nn2:=cadr nn2$
l21:add2sint(nv,nv1,x,x1,n1,n2,nn1,nn2)$
  go to r$
l3:y2:=caaddr gds$
  x2:=atsoc(y2,asl)$
  if null x2 then go to er1$
  x2:=tdifpair cdr x2$
  if !*twogrid and flagp(x2,'twogrid) then go to l20$
  nv2:=ngetvar(x2,'nrvar)$
  n1:=atsoc(acmemb(list(y . 'one,y1 . 'one,y2 . 'one),lgrids),agdsl)$
  n2:=atsoc(acmemb(list(y . 'half,y1 . 'half,y2 . 'half),lgrids),agdsl)$
  nn1:=atsoc(acmemb(list(y . 'one,y1 . 'one,y2 . 'half),lgrids),agdsl)$
  nn2:=atsoc(acmemb(list(y . 'half,y1 . 'half,y2 . 'one),lgrids),agdsl)$
  n3:=atsoc(acmemb(list(y . 'one,y1 . 'half,y2 . 'one),lgrids),agdsl)$
  n4:=atsoc(acmemb(list(y . 'half,y1 . 'one,y2 . 'half),lgrids),agdsl)$
  n5:=atsoc(acmemb(list(y . 'one,y1 . 'half,y2 . 'half),lgrids),agdsl)$
  n6:=atsoc(acmemb(list(y . 'half,y1 . 'one,y2 . 'one),lgrids),agdsl)$
  if n1 and n2 and nn1 and nn2 and n3 and n4 and n5 and n6 then t
    else go to er2$
  n1:=cadr n1$
  n2:= cadr n2$
  nn1:=cadr nn1$
  nn2:=cadr nn2$
  n3:=cadr n3$
  n4:=cadr n4$
  n5:=cadr n5$
  n6:=cadr n6$
  if n1=nn1 and n2=nn2 and n3=n5 and n4=n6 then
      <<n2:=n3$
        nn1:=nn2$
        nn2:=n4$
        go to l21 >>
    else if n1=n3 and n2=n4 and nn1=n5 and nn2=n6 then
      <<n2:=nn1$
        nn1:=n4$
        x1:=x2$
        nv1:=nv2$
        go to l21 >>
    else if n1=n6 and n2=n5 and nn1=n4 and nn2=n3 then
      <<n2:=nn2$
        nn1:=n5:
        nn2:=n4$
        x:=x2$
        nv:=nv2$
        go to l21 >>$
  add3sint(nv,nv1,nv2,x,x1,x2,n1,n2,n3,n4,n5,n6,nn1,nn2)$
r:return (nil . 1)$
er:msgpri(nil,l," Free vars not yet implemented ",nil,'hold)$
  go to r$
er1:msgpri(" Failed matching of variables in  ",
        u,list(asl,agdsl),nil,'hold)$
  go to r$
er2:msgpri(" All grids not given for term  ",u,list(asl,agdsl),
                   nil,'hold)$
  go to r
end$

procedure add2sint(nv,nv1,x,x1,n1,n2,nn1,nn2)$
begin
% Enhansment for symmetries, when only one variable influence
  if n1=n2 and nn1=nn2 then
      <<f2plus(nvar,nv,0,6*n1)$
        f2plus(nvar,nv,1,6*nn1)$
        return >>
    else if n1=nn2 and n2=nn1 then
      <<f2plus(nvar,nv1,0,6*n1)$
        f2plus(nvar,nv1,1,6*n2)$
        return >>$
  n1:=3*n1$
  n2:=3*n2$
  nn1:=3*nn1$
  nn2:=3*nn2$
  x:=list(x,x1)$
  f2plus(nvar,nv,0,list(x,n1,n2))$
  f2plus(nvar,nv,1,list(x,nn1,nn2))$
  x:=reverse x$
  f2plus(nvar,nv1,0,list(x,n1,nn2))$
  f2plus(nvar,nv1,1,list(x,nn1,n2))$
  return
end$

procedure add3sint(nv,nv1,nv2,x,x1,x2,n1,n2,n3,n4,n5,n6,nn1,nn2)$
begin
  n1:=2*n1$
  n2:=2*n2$
  nn1:=2*nn1$
  nn2:=2*nn2$
  n3:=2*n3$
  n4:=2*n4$
  n5:=2*n5$
  n6:=2*n6$
  x:=list(x,x1,x2)$
  f2plus(nvar,nv,0,list(x,n1,nn1,n3,n5))$
  f2plus(nvar,nv,1,list(x,n2,nn2,n4,n6))$
  f2plus(nvar,nv1,0,list(x,n1,nn1,n4,n6))$
  f2plus(nvar,nv1,1,list(x,n2,nn2,n3,n5))$
  f2plus(nvar,nv2,0,list(x,n1,nn2,n3,n6))$
  f2plus(nvar,nv2,1,list(x,n2,nn1,n4,n5))$
  return
end$

procedure simpinterpol u$
begin
  scalar asl,agdsl,gds,x,y,xx,a$
  u:=prepsq simp car u$
  if eqcar(u,'diff) and not coord memq cddr u then
%  !!!! Be aware  !!!! could not work for mixed derivatives
      return <<asl:=!*exp$
               !*exp:=t$
               u:= simp formlnr('diff . (prepsq simp formlnr
                       list('intt,cadr u,coord) . cddr u))$
               !*exp:=asl$
               u>>$
  asl:=finddifm u$
  if null asl then return (nil . 1)
    else if eqcar(asl,'special) then return cdr asl$
  agdsl:=cadr asl$   %  Actual list from DIFML!*, contains definition
                     %  of grid, penalty and diff. scheme
  asl:=car asl$      %  Assoc. list of assignments of variables X,U,V,W
                     %  to actual variables
  if not gvar memq grids!* then go to erg$
  asl:=append(asl,get('steps,gvar))$   % Adding DIM1, DIP1 ... to assoc.
                                       % list
  if null caar agdsl then return simp sublap(asl,caddar agdsl)$  % For
  a:=caar agdsl$                       % DIFMATCH without def. grids
b:if null a then go to c$
  y:=caar a$
  x:=atsoc(y,asl)$
  if null x then go to er1$            % GDS is assoc. list of actual
  xx:=cdr x$                           % assignments of grids to
  x:=getgrid xx$                       % variables U, V
  if gvar eq 'half then x:=gnot x$
  if !*twogrid and twogridp xx then t
    else gds:=(y . x) . gds$
  a:=cdr a$
  go to b$
c:if null gds then go to a$      % For given functions which can be on
  y:=get('difm!*,'grids)$        % both grids
  x:=acmemb(gds,y)$              % Unique GDS
  if null x then go to er1$
  gds:=x$
a:x:=assoc(gds,agdsl)$
  if null x then go to erg$
  x:=caddr x$                    % Actual difference scheme
  return simp sublap(asl,x)$
er1:msgpri(" Failed matching of ",u,list(asl,agdsl,gds),nil,
                  'hold)$
  return (nil . 1)$
erg:msgpri(" Bad grids ",u,list(asl,agdsl,gds),nil,'hold)$
  return (nil . 1)
end$

procedure twogridp u$
% Checks if prefix form U can be on both grids
begin
  scalar x$
  return
    if atom u then
        if flagp(u,'twogrid) then
            if !*twogrid and u memq given!* and
               getel list(get(u,'grid),ncor) then nil
              else t
          else nil
      else if flagp(car u,'twogrid) then
        if !*twogrid and car u memq given!* and
           getel list(get(car u,'grid),ncor) then nil
          else t
      else if car u memq '(diff plus difference) then twogridp cadr u
      else if car u eq 'times then twogridpti cdr u
      else nil
end$

procedure twogridpti u$
begin
  scalar x$
a:x:=twogridp car u$
  if x then return x$
  u:=cdr u$
  if u then go to a$
  return nil
end$


procedure getgrid u$
begin
  scalar x$
  return
    if atom u then
        if x:=get(u,'xgrid) then x
          else if !*twogrid and u memq given!* and
                  (x:=getel list(get(u,'grid),ncor)) then car x
          else nil
      else if (x:=get(car u,'xgrid)) then x
      else if !*twogrid and car u memq given!* and
               (x:=getel list(get(car u,'grid),ncor)) then car x
      else if car u eq 'diff then
               if atom cadr u then getgrid cadr u
                %else if caadr u eq 'times then           % probably can
                %         if (x:=getgrid cadadr u) then x % be deleted
                %           else getgrid caddr cadr u     % !!!!!"!!!
                 else getgrid cadr u
      else if car u memq '(plus difference) then getgrid cadr u
      else if car u eq 'times then getgti cdr u
      else nil
end$

procedure getgti u$
begin
  scalar x$
a:x:=getgrid car u$
  if x then return x$
  u:=cdr u$
  if u then go to a$
  return nil
end$

procedure sublap(u,v)$
% U is assoc. list, V is pattern diff. scheme
% Performs substitution of assod. list into the diff. scheme
begin
  scalar x$
  return
    if null u or null v then v
      else if atom v then
        if numberp v then v
          else if x:=atsoc(v,u) then cdr x
          else v
      else if flagp(car v,'app) then sublap1(u,v)
      else (sublap(u,car v) . sublap(u,cdr v))
end$

flag('(u f v w x g),'app)$

procedure sublap1(u,v)$
begin
  scalar x,y$
  x:=atsoc(car v,u)$
  if null x then return msgpri(" Substitution for ",v," not find",
                               nil,'hold)$
  x:=cdr x$
  y:=mapcar(cdr v,function(lambda a$irev sublap(u,a)))$
  return
    if eqcar(x,'diff) then ('diff . (subappend(cadr x,y) . cddr x))
      else subappend(x,y)
end$

procedure subappend(x,y)$
if atom x then if x memq iobjs!* and depends(x,coord) then (x . y)
                 else x
  else if car x memq iobjs!* and depends(car x,coord) then append(x,y)
         else (car x . mapcar(cdr x,function(lambda a$
                                     subappend(a,y))))$

procedure irev u$
begin
  u:=simp u$
  return
    if cdaaar u=1 and cdaar u=cdr u and fixp cdar u then
        if cdr u=1 then
            if cdar u<0 then list('difference,
                                   caaaar u,
                                   -cdar u)
               else list('plus,
                         caaaar u,
                         cdar u)
           else if cdar u<0 then list('difference,
                                      caaaar u,
                                      list('quotient,
                                           -cdar u,
                                           cdr u))
                  else list('plus,
                            caaaar u,
                            list('quotient,
                                 cdar u,
                                 cdr u))
      else prepsq u
end$

unfluid '(coord unvars sunvars interpolp novars ncor nvar intp icor gvar
        hi hip1 hip2 him1 him2 lhs rhs lsun lun xxgrid resar)$


procedure gentempst$
list('gentemp,xread t)$

put('gentemp,'stat,'gentempst)$
put('gentemp,'formfn,'formgentran)$

put('outtemp,'stat,'endstat)$
flag('(outtemp),'eval)$


algebraic$

endmodule$

end$