reduce-historical/r36/plot.red

From hearn@rand.orgSat Feb  3 09:55:22 1996
Date: Sat, 03 Feb 96 01:00:05 -0800
From: Tony Hearn <hearn@rand.org>
To: shar-list@rand.org
Subject: Shar File


# This is a shell archive.  Remove anything before this line, then
# unpack it by saving it in a file and typing "sh file".  (Files
# unpacked will be owned by you and have default permissions.)
#
# This archive contains:
# plot/plot.red

echo x - plot/plot.red
if [ -f plot/plot.red ]
   then
    mv plot/plot.red \
       plot/plot.red.old
  else
    echo  "*** New file plot/plot.red created"
fi
cat > "plot/plot.red" \
    << '//E*O*F plot/plot.red//'
module plot; % device and driver independent plot services.

% Author: Herbert Melenk.

% Minor corrections by Winfried Neun (October 1995)

create!-package('(plot),nil);

global '(

      plotdriver!*         % modulename of the actual driver.

      plotmax!*            % maximal floating point number which
                           % gnuplot supports on the machine
                           % (mostly IEEE double or single precision).

      plotmin!*            % lower bound (=1/plotmax!*)

      variation!*          %  defintion of y-bigstep for smooth

      plotoptions!*        % list for collecting the options.

);

fluid '(

      plotderiv!*          %  derivative for 2d plot

);

!#if(or (errorp (errorset '!*writingfaslfile nil nil)) 
        (not !*writingfaslfile) 
        (errorp (errorset '(load fcomp) nil nil)))
   prin2t "no support for fast float!"; 
   eval'(dm fdeclare (x) nil); 
   eval'(dm thefloat (x)(cadr x));
!#endif

endmodule;

module plotsynt; % Support for the syntax of the plot command.

% Author: Herbert Melenk.

% Create .. as the infix operator if not yet done.

!*msg := nil;  % prevent message  ".. redefined" during load

newtok '( (!. !.) !*interval!*);  

if not(gettype '!*interval!* = 'operator) then
<< 
   precedence .., or;
   algebraic operator ..;
   put('!*interval!*,'PRTCH,'! !.!.! );
>>;

mkop 'point;

!*msg := t;

fluid '(plot!-points!* plot!-refine!* plot!-contour!*);

global '(plot_xrange plot_yrange plot_zrange);
share plot_xmesh,plot_ymesh,plot_xrange,plot_yrange,plot_zrange;

fluid '(plotprecision!*);

plotprecision!* := 0.9995;

fluid '(!*show_grid test_plot);

switch show_grid;
switch test_plot; % for test printouts

if null plotmax!* then
<< 
   load!-package 'arith;
   if not !!plumax then roundconstants();
   plotmax!* := !!plumax;     % IEEE double precision 
>>;

plotmin!*:= 1.0/plotmax!*;


fluid '(plotranges!* plotfunctions!*  plotstyle!* !*plotoverflow
        !*roundbf);

put('plot,'psopfn,'ploteval);

% I need the following definition only at compile time.
macro procedure plotdriver u;
  {'apply,{'get,'plotdriver!*,mkquote cadr u},'list.cddr u};

symbolic procedure ploteval u;
  begin scalar m,!*exp;
    if null plotdriver!* then 
      rederr "no active device driver for PLOT";
    m:=plotrounded(nil);
    plot!-points!* := {20};
    plot!-refine!* := 8;
    !*plotoverflow := nil;
    plotranges!* := plotfunctions!* := nil;
    plotstyle!* := 'lines;
    plotdriver(init);
    for each option in u do ploteval1 plot!-reval option;
    errorset('(ploteval2),t,nil);
    plotrounded(m);
  end;
 
symbolic procedure plot!-reval u;
 % Protected call reval: simplify u, but don't call any 
 % algebraic procedure.
   begin scalar w;
     w:={nil};
     u:=plot!-reval1(u,w);
     return car w and u or reval u;
   end; 

symbolic procedure plot!-reval1(u,w);
  if idp u then reval u else
  if atom u or eqcar(u,'!:dn!:) or get(car u,'dname) then u else   %WN
  if eq(car u,'!*sq) then plot!-reval1 (reval u,w) else            %WN
  <<if flagp(car u,'opfn) and
         memq(car u,'(first second rest rhs lhs)) then
     << u := reval u;                                              %WN
        plot!-reval1(u,w)>> else
   << if flagp(car u,'opfn) then car w:=t;
    car u . for each q in cdr u collect plot!-reval1(q,w) >> >>;

symbolic procedure ploteval1 option;
   begin scalar x,do;
     do := get(plotdriver!*,'do);
     if pairp option and (x:=get(car option,do)) 
             then apply(x,list option) else
     if pairp option and (x:=get(car option,'plot!-do)) 
             then apply(x,list option) else
     if eqcar(option,'equal) and (x:=get(cadr option,do))
             then apply(x,list caddr option) else
     if eqcar(option,'equal) and (x:=get(cadr option,'plot!-do))
             then apply(x,list caddr option) 
       else ploteval0 option;
   end;
  
symbolic procedure ploteval0 option;
  begin scalar l,r,opt,w;
    opt:=get(plotdriver!*,'option);
    if flagp(option,opt) then
      <<plotoptions!*:=option . plotoptions!*; return>>;
    if eqcar(option,'list) then
      <<option := cdr option;
        if option and eqcar(car option,'list) then
          return (plotfunctions!*:= 
             ('points.plotpoints option).plotfunctions!*);
        for each o in option do ploteval0 o; return;
      >>;
    if eqcar(option,'equal) and flagp(cadr option,opt) then
      <<plotoptions!*:=(cadr option.caddr option). plotoptions!*; 
       return>>;
    if not eqcar(option,'equal) then 
      <<plotfunctions!*:= (nil.option) . plotfunctions!*; return>>;

      % Handle equations.
    l:=plot!-reval cadr option;
    r:=plot!-reval caddr option;
    if plot!-checkcontour(l,r) then return
      plotfunctions!*:=('implicit.l) . plotfunctions!*;
    if not idp l then typerr(option,"illegal option in PLOT");

    if l memq '(size terminal view) then
      <<plotoptions!*:=(l.r).plotoptions!*; return>>;

       % iteration over a range?
    if eqcar(r,'times) and eqcar(caddr r,'!*interval!*) 
        and evalnumberp(w:=cadr r) and evalgreaterp(w,0) and
        not evalgreaterp(w,1) 
     then <<plot!-points!*:=append(plot!-points!*,
              {l.reval{'floor,{'quotient,1,w}}});
            r:=caddr r>>;

    if eqcar(r,'quotient) and eqcar(cadr r,'!*interval!*) 
       and fixp caddr r and caddr r > 0
     then <<plot!-points!*:=append(plot!-points!*,{l.caddr r});
            r:=cadr r>>;

       % range?
    if eqcar(r,'!*interval!*) then
     <<r:='!*interval!* . revalnuminterval(r,t);
       plotranges!* := (l . r) . plotranges!*>>
      else
       plotfunctions!* := (l . r) . plotfunctions!*;
  end;

symbolic procedure ploteval2 ();
   % all options are collected now;
  begin scalar dvar,ivars,para,impl;
   for each u in plotfunctions!* do 
     <<impl:=impl or car u eq 'implicit;
       para:=eqcar(cdr u,'point);
       if impl and dvar and dvar neq car u then
          rederr "mixture of implicit and regular plot not supported";
       dvar:=car u or dvar;
       ivars := plotindepvars(cdr u,ivars)>>;
      % classify
   if null dvar then
   <<dvar:='(x y z);
     for each x in ivars do dvar:=delete(x,dvar);
     if dvar then dvar:=if 'y memq dvar then 'y else car dvar;
   >>;
   if para and length ivars=1 then plotevalpara1(car ivars) else
   if para and length ivars=2 then plotevalpara2(car ivars,cadr ivars)
    else if length ivars=1 then ploteval2x(car ivars,dvar) else
   if length ivars=2 then ploteval3xy(car ivars,cadr ivars,dvar) else
   if length ivars=3 and impl then 
            ploteval3impl(car ivars,cadr ivars,caddr ivars)
    else typerr('list . for each p in plotfunctions!* collect 
                         if null car p then cdr p else
                         {'equal,car p,cdr p},
                " plot option or function");
  plotdriver(show);
  end;

symbolic procedure plot!-checkcontour(l,r);
  % true if the job is a contour expression.
   if length plotindepvars(l,nil)=2 then
     if r=0 then <<plot!-contour!*:={0};t>>
       else eqcar(r,'list) and
     <<plot!-contour!*:= for each x in cdr r collect   
        <<x:=plot!-reval x; l:=l and adomainp x; x>>;
       l>>;

symbolic procedure plotrange(x,d);
   begin scalar y;
     y:=assoc(x,plotranges!*);
     y:=if y then cdr y else d;
     if y=0 or null y then % return nil;
     y:={'!*INTERVAL!*, - plotmax!*, plotmax!*};
     if not eqcar(y,'!*INTERVAL!*) then
        typerr(y,"plot range");
     return {plotevalform0(rdwrap cadr y,nil) , 
             plotevalform0(rdwrap caddr y,nil)};
   end;

symbolic procedure plot!-points(x);
  (if w then cdr w else car plot!-points!*)
    where w=assoc(x,cdr plot!-points!*);

symbolic procedure plotseteq(u,v);
    null u and null v or car u member v 
       and plotseteq(cdr u,delete(car u,v));

symbolic procedure plotindepvars(u,v);
    if idp u then 
       if member(u,v) or member(u,'(e pi)) 
                      or u eq 'i and !*complex then v 
               else  u . v 
    else if eqcar(u,'file) then cddr u 
    else if pairp u then
      if eqcar(u,'!:dn!:) or get(car u,'dname) then v else
      if member(car u,'(plus minus difference times quotient expt)) or
	 get(car u,'!:RD!:) or get(car u,'simpfn)
	    or eqcar(getd(car u),'expr)
       then <<for each x in cdr u do v:=plotindepvars(x,v); v>>
     else typerr(u,"expression in function to plot")
    else v;

remprop('plotshow,'stat);
 
symbolic procedure plotshow();
   plotdriver(show);

put('plotshow,'stat,'endstat);

remprop('plotreset,'stat);

symbolic procedure plotreset();
   plotdriver(reset);

put('plotreset,'stat,'endstat);

put('points,'plot!-do,
    function(lambda(x);car plot!-points!*:=ieval x));

put('refine,'plot!-do,
    function(lambda(x);plot!-refine!*:=ieval x));

endmodule; % plotsynt.

module plotexp2; % Compute explicit 2-dim graph y=f(x);

symbolic procedure ploteval2x(x,y);
  begin scalar xlo,xhi,ylo,yhi,rx,ry,fp,fcn,fcns,pts;
     if y='implicit then 
        rederr "no implicit plot in one dimension";  
     rx:=plotrange(x,
       reval(plot_xrange or '(!*interval!* -10 10)));
     xlo:=car rx;
     xhi:=cadr rx;
     fcns:= reverse plotfunctions!*;
     ry:=plotrange(y, reval(plot_yrange or nil));
     if ry then <<ylo:=car ry; yhi:=cadr ry>>;
     while fcns do
     <<fcn := car fcns; fcns := cdr fcns;
       if eqcar(fcn,'points) then fp:=caddr fcn else
       pts:=ploteval2x1(cdr fcn,x,xlo,xhi,ylo,yhi).pts;
     >>;
     plotdriver(plot!-2exp,x,y,pts,fp);
  end;

symbolic procedure ploteval2x1(f,x,xlo,xhi,ylo,yhi);
   begin scalar plotsynerr!*,l,d,d0,u,v,vv,p,mx,mn,ff;
        scalar plotderiv!*;
        integer nx;
     % compute algebraic derivative.
    ff:= errorset({'reval,mkquote {'df,f,x}},nil,nil); 
    if not errorp ff and not smemq('df,car ff) then 
% Hier irrte Goethe.  % This comment is for Herbert, please keep it
%      plotderiv!*:= rdwrap plotderiv!* . plotderiv!*;
       plotderiv!*:= rdwrap car ff . car ff;
    ff:=rdwrap f;
    p:=float (nx:=plot!-points(x));
    d:=(d0:=(xhi-xlo))/p;
    v:=xlo;
    for i:=0:nx  do
     <<vv:=if i=0 or i=nx  then v 
           else v+d*(random(100)-50)*0.001;
       u:= plotevalform(ff,f,{x.vv});
       if plotsynerr!* then typerr(f,"function to plot");
       if eqcar(u,'overflow) then u:=nil;
       if u then
       << 
          if yhi and u>yhi then u:=yhi else
          if ylo and u<ylo then u:=ylo;;
          if null mx or u>mx then mx:=u;
          if null mn or u<mn then mn:=u
       >>;
       l:=(vv.u).l;
       v:=v+d;
     >>;
     if null mx or null mn then rederr "plot, sampling failed";
     variation!* :=
     if yhi and not(yhi=plotmax!*) then (yhi-ylo) else
        (mx-mn); %  ploteval2xvariation l;

    if plot!-refine!*>0 then
        l:=smooth(reversip l,ff,f,x,mx,mn,ylo,yhi,d);
    return {for each x in l collect {car x,cdr x}};
   end;


symbolic procedure ploteval2xvariation l;
  begin scalar u;
   % compute geometric mean value.
    integer m;
    u:=1.0;
    for each p in l do
     <<m:=m+1; p:=cdr p;
       if p and p neq 0.0 then u:=u*abs p;
     >>;
    return expt(u,1/float m);
  end;
     
symbolic procedure smooth(l,ff,f,x,maxv,minv,ylo,yhi,d);
  begin scalar rat,grain,cutmax,cutmin,z,z0;
   z:=l;
   cutmax :=  yhi or (- 2*minv + 3*maxv);
   cutmin :=  ylo or (3*minv - 2*maxv);
   grain  :=  variation!* *0.05;
   rat := d / if numberp maxv and maxv > minv then (maxv - minv) 
               else 1;
    % Delete empty points in front of the point list.
   while z and null cdar z and cdr z and null cdadr z do z:=cdr z;
    % Find left starting point if there is no one yet.
   if z and null cdar z and cdr z then
   <<z0:= findsing(ff,f,x,ylo,yhi,cadr z,car z);
     if z0 then l:=z:=z0.cdr z; 
   >>;
   while z and cdr z do
   <<z0:=z; z:=cdr z;
       smooth1(z0,ff,f,x,cutmin,cutmax,grain,rat,0,plot!-refine!*)
   >>;
   return l;
  end;

symbolic procedure smooth1(l,ff,f,x,minv,maxv,g,rat,lev,ml);
    smooth2(l,ff,f,x,minv,maxv,g,rat,lev,ml);


symbolic procedure smooth2(l,ff,f,x,minv,maxv,g,rat,lev,ml);
  if lev >= ml then 
      smooth3(l,ff,f,x,minv,maxv,g,rat,lev,ml)
    else
  if null cdar l then t else
 begin scalar p0,p1,p2,p3,x1,x2,x3,y1,y2,y3,d;
       scalar dy12,dy32,dx12,dx32,z,a,w;
%%%%%    fdeclare(x1,x2,x3,y1,y2,y3,rat,d,dx12,dx32,dy12,dy32);
    lev:= add1 lev;
    p1:=car l; p3:=cadr l;
    x1:=car p1; y1:=cdr p1;
    x3:=car p3; y3:=cdr p3;
  nopoint:
    if null y3 then
    <<if null cddr l then return(cdr l:=nil);
      x2:=x3; y2:=y3; cdr l:=cddr l;
      p3:=cadr l; x3:=car p3; y3:=cdr p3;
      if y3 then goto outside else goto nopoint;
    >>;
    % Generate a new point
    x2:=(x1+x3)*0.5;
    if null (y2 := plotevalform(ff,f,{x.x2})) 
       or eqcar(y2,'overflow) then goto outside;
    if y2 < minv or y2 > maxv then goto outside;

    dy32 := (y3 - y2) * rat; dx32 := x3 - x2;
    dy12 := (y1 - y2) * rat; dx12 := x1 - x2;
    d :=  sqrt((dy32**2 + dx32**2) * (dy12**2 + dx12**2));
    if zerop d then return t;
    w := (dy32*dy12 + dx32*dx12);
    d:= errorset({'quotient,w,d},nil,nil);
     % d is cos of angle between the vectors p2p1 and p2p3.
     % d near to 1 means that the angle is almost 0 (needle).
    if ploterrorp d then goto disc else d:=car d;
    a:=abs(y3-y1)<g;
    if d > plotprecision!* and a and lev=ml then goto disc;

      % I have a good point, insert it with destructive replacement.
    cdr l := (x2 . y2) . cdr l;
      % When I have an almost 180 degree angle, I compare the
      % derivative at the midpoint with the difference quotient.
      % If these are close to each other, I stop the refinement.
    if -d  > plotprecision!* 
       and (null plotderiv!* 
             or
           ((w:=plotevalform(car plotderiv!*,cdr plotderiv!*,{x.x2}))
                and abs (w - (y3-y1)/(x3-x1))*rat <0.1))
      then return t;
    smooth2(cdr l,ff,f,x,minv,maxv,g,rat,lev,ml);
    smooth2(l,ff,f,x,minv,maxv,g,rat,lev,ml);
    return t;

    % Function has exceeded the boundaries. I try to locate the screen 
    % crossing point.
  outside:
    cdr l := (x2 . nil) . cdr l;
    z := cdr l;    % insert a break
    p2:= (x2 . y2);  % artificial midpoint

    p0:=findsing(ff,f,x, minv, maxv, p1, p2);
    if p0 then
      << cdr l := p0 . z;
         smooth2(l,ff,f,x,minv,maxv,g,rat,lev,ml) >>;
    p0 := findsing(ff,f,x, minv, maxv, p3, p2);
    if p0 then
      << cdr z := p0 . cdr z;
         smooth2(cdr z,ff,f,x,minv,maxv,g,rat,lev,ml) >>;
    return;
  
  disc:  % look for discontinuities.
    return smooth3(l,ff,f,x,minv,maxv,g,rat,lev,ml);
  end;

symbolic procedure smooth3(l,ff,f,x,minv,maxv,g,rat,lev,ml);
 % detect discontinuities.
 begin scalar p1,p2,p3,x1,x2,x3,y1,y2,y3,d;
       scalar x2hi,y2hi,x2lo,y2lo,dir,hidir,chi,clo;
       scalar lobound,hibound;
       integer n;
 %%%%%      fdeclare(x1,x2,x3,y1,y2,y3,d,hidir);
    p1:=car l; p3:=cadr l;
    x1:=car p1; y1:=cdr p1;
    x3:=car p3; y3:=cdr p3;
    if abs(y3-y1)<variation!* then return t;

  % Bigstep found. 
    hibound:=variation!**10.0; lobound:=-hibound;
    if y1>y3 then
    <<x2hi:=x1; y2hi:=y1; x2lo:= x3; y2lo:=y3; hidir:=-1.0>>
       else
    <<x2hi:=x3; y2hi:=y3; x2lo:= x1; y2lo:=y1; hidir:=1.0>>;
    n:=0; d:= (x3-x1)*0.5; x2:=x1+d;
  % intersection loop. Cut remaining interval into two parts.
  % If there is again a high increase in one direction we assume
  % a pole.
  next_point:
    if null (y2 := plotevalform(ff,f,{x.x2}))
       or eqcar(y2,'overflow) then goto outside;
    if y2 < minv then 
      <<x2lo:=x2;y2lo:=minv;dir:=hidir>> 
    else if y2 < y2lo then 
      <<if y2lo<0.0 and y2<y2lo+y2lo and y2<lobound then clo:=t;
        x2lo:=x2;y2lo:=y2;dir:=hidir;>> 
    else if y2 > maxv then 
     <<x2hi:=x2;y2hi:=maxv;dir:=-hidir>> 
    else if y2 > y2hi then 
     <<if y2hi>0.0 and y2>y2hi+y2hi and y2>hibound then chi:=t;
       x2hi:=x2;y2hi:=y2;dir:=-hidir;>> else
      goto no_disc;
    if dir and (n:=n+1)<20 and (not clo or not chi) then
    <<d:=d/2; x2:=x2+d*dir; goto next_point>>;
  no_disc:
    if null dir then return t;
    if clo then y2lo:=minv;
    if chi then y2hi:=maxv;
  outside:
    p1:=(x2hi.y2hi); p3:=(x2lo.y2lo);
    if hidir=1.0 then <<p2:=p3;p3:=p1;p1:=p2>>;
    cdr l := p1 . (car p1.nil) . p3 . cdr l;
    return;
  brk:
    cdr l := (car p1.nil)  . cdr l;
    return;
  end;
    
symbolic procedure findsing(ff,f,x,minv,maxv,p1,p3);
% P3 is a point with a singularity.
% Try to locate the point where we exceed minv/maxv by subdivision.
begin scalar p1, p2, p3, x1, x2, x3, y1, y2, y3, d, x2n;
    x1:=car p1; y1:=cdr p1;
    x3:=car p3; y3:=cdr p3;
    d := (x3-x1)*0.5; x2n:=x1+d;
    for i:=1:5 do
    << d:=d*0.5; x2:= x2n;
       if null(y2 := plotevalform(ff,f,{x.x2}))
          or eqcar(y2,'overflow)
          or y2 < minv or y2 > maxv
        then x2n := x2n - d
        else << p2 := (x2 . y2); x2n := x2n + d>>
    >>;
   if null p2 then return nil;

    % generate uniform maxima/minima
    x2 := car p2; y2 := cdr p2;
    y2 := if (eqcar(y2,'overflow) and cdr y2<0) or y2<minv 
          then minv 
     else if eqcar(y2,'overflow) or y2>maxv then maxv else y2;
    return (x2 . y2)
end;

endmodule;   % plotexp2

module pltpara; % Computing parametric curve.
                 % (x,y) = (x(t),y(t))
                 %     or
                 % (x,y,z) = (x(t),y(t),z(t))
                 %     or
                 % (x,y,z) = (x(t,u),y(t,u),z(t,u))

% Author: Herbert Melenk, ZIB Berlin.

symbolic procedure plotevalpara1(x);
  begin scalar xlo,xhi,ylo,yhi,rx,ry,fcn,fcns,pts;
     rx:=plotrange(x,
       reval(plot_xrange or '(!*interval!* -10 10)));
     xlo:=car rx;
     xhi:=cadr rx;
     fcns:= reverse plotfunctions!*;
     % ry:=plotrange(y, reval(plot_yrange or nil));
     if ry then <<ylo:=car ry; yhi:=cadr ry>>;
     while fcns do
     <<fcn := cddar fcns; fcns := cdr fcns;
       pts:=plotevalpara11(fcn,x,xlo,xhi).pts;
     >>;
     if length fcn=2 then  plotdriver(plot!-2exp,'x,'y,list pts,nil)
        else plotdriver(plot!-3exp!-reg,'x,'y,'z,pts)
   end;

symbolic procedure plotevalpara11(fm,x,xlo,xhi);
   begin scalar plotsynerr!*,l,d,d0,u,v,p,fl;
        scalar plotderiv!*;
        integer nx;
    fl:= for each f in fm collect rdwrap f.f;
    p:=float (nx:=plot!-points(x));
    d:=(d0:=(xhi-xlo))/p;
    v:=xlo;
    for i:=0:nx  do
     <<u:= for each f in fl collect plotevalform(car f,cdr f,{x.v});
       if plotsynerr!* then typerr(fm,"function to plot");
       if smemq('overflow,u) then u:=nil;
       l:=u.l;
       v:=v+d;
     >>;
    return reversip l;
   end;

symbolic procedure plotevalpara2(x,y);
  begin scalar xlo,xhi,ylo,yhi,rx,ry,fcn,fcns,pts;
     rx:=plotrange(x,
       reval(plot_xrange or '(!*interval!* -10 10)));
     xlo:=car rx; xhi:=cadr rx;
     fcns:= reverse plotfunctions!*;
     ry:=plotrange(y, 
       reval(plot_yrange or '(!*interval!* -10 10)));
     ylo:=car ry; yhi:=cadr ry;
     fcn := cddar fcns; fcns := cdr fcns;
     if length fcn neq 3 then typerr(cdar fcns,"function to plot");
     pts:=plotevalpara21(fcn,x,xlo,xhi,y,ylo,yhi);
     plotdriver(plot!-3exp!-reg,'x,'y,'z,pts)
   end;

symbolic procedure plotevalpara21(fm,x,xlo,xhi,y,ylo,yhi);
   begin scalar plotsynerr!*,l,ll,dx,dy,u,v,p,fl,w;
        scalar plotderiv!*;
        integer nx,ny;
    fl:= for each f in fm collect rdwrap f.f;
    p:=float(nx:=plot!-points(x));
    dx:=(xhi-xlo)/p; 
    p:=float(ny:=plot!-points(y));
    dy:=(yhi-ylo)/p;
    v:=xlo;
    for i:=0:nx  do
     <<w:= ylo; l:=nil;
       for j:=0:ny do 
       <<u:= for each f in fl collect 
              plotevalform(car f,cdr f,{x.v,y.w});
         if plotsynerr!* then typerr(fm,"function to plot");
         if smemq('overflow,u) then u:=nil;
         l:=u.l;
         w:=w+dy
       >>;
       v:=v+dx;
       ll:=l.ll;
     >>;
    return ll;
   end;

endmodule;

module plotexp3; % Computing surface z=f(x,y) with regular grid.

% Author: Herbert Melenk, ZIB Berlin.

% A rectangular grid is encoded as list of rows.
% A row is a list of points.
% A point is a list {v,h,x,y,z} where
%   z,y are the coordinates and z is the function value.
% The boolean values v ("vertical" = y direction) and 
%    h ("horizontal" = x direction) are true, 
%   if the connection to the neighbour point in that
%   direction is valid, nil if the connection is broken.

symbolic procedure ploteval3xy(x,y,z);
  begin scalar rx,ry,rz,f,fcn;     
     rx:=plotrange(x,
      reval(plot_xrange or '(!*interval!* -10 10)));
     ry:=plotrange(y,
      reval(plot_yrange or '(!*interval!* -10 10)));
     rz:=plotrange(z, reval(plot_zrange or nil));
     fcn := car reverse plotfunctions!*;
     if z='implicit then 
       return ploteval2xyimpl(rx,ry,fcn,x,y);
     f:= if eqcar(fcn,'points) then caddr fcn else
        ploteval3xy1(cdar plotfunctions!*,
         z,if rz then car rz, if rz then cadr rz,
         x,car rx,cadr rx,
         y,car ry,cadr ry);
     plotdriver(plot!-3exp!-reg,x,y,z,f);
  end;

symbolic procedure ploteval3xy1(f,z,zlo,zhi,x,xlo,xhi,y,ylo,yhi);
   begin scalar u,l,ll,ff,p,r,w;
%        scalar x1,x2,y1,y2,xx,yy,p1,p2,p3,p4,f1,f2,f3,f4;
        integer nx,ny;
     ff := rdwrap f;
     xlo:=rdwrap xlo; xhi:=rdwrap xhi;
     ylo:=rdwrap ylo; yhi:=rdwrap yhi;
     nx:=plot!-points(x); ny:=plot!-points(y);
       % compute the basic grid.
     r:=ploteval3xy1pts(f,ff,z,zlo,zhi,x,xlo,xhi,nx,y,ylo,yhi,ny);
     l:=cdr r;
     w:=car r;
     r:={l};
        % create refined environments for the bad points
     for each q in w do
      r:=cdr 
       ploteval3xy1pts(f,ff,z,zlo,zhi,x,car q,cadr q,4,y,caddr q,
		       cadddr q,4)
         .r;
%       % look for singularities or points of big variation
%     ll:=l;
%     while ll and cdr ll do
%     <<p := pair(car ll,cadr ll); ll:=cdr ll;
%       while p and cdr p do 
%       <<p1:=caar p; p2:=caadr p; p3:=cdar p; p4:=cdadr p; p:=cdr p;
%         if (f1:=car cdddr p1) and (f2:=car cdddr p2)
%                 and (f3:=car cdddr p3) and (f4:=car cdddr p4) then
%         <<xx:=((x1:=caddr p1) + (x2:=caddr p2))*0.5;
%           yy:=((y1:=cadddr p1) + (y2:=cadddr p3))*0.5;
%           u:=plotevalform(ff,f,{x . xx,y . yy});
%           if null u or eqcar(u,'overflow) or
%             numberp u and abs u > (abs f1+abs f2+abs f3+abs f4)*0.5
%              then 
%              <<r:=cdr
%                    ploteval3xy1pts(f,ff,z,zlo,zhi,x,x1,x2,3,y,y1,y2,3)
%                 .r;
%                   % cut connections
%                %car p1:= nil; cadr p1:= nil; 
%              >>;
%         >>;
%       >>;
%     >>;
     return ploteval3xy3 r;
   end;

symbolic procedure ploteval3xy1pts
          (f,ff,z,zlo,zhi,x,xlo,xhi,nx,y,ylo,yhi,ny);
  % Compute orthogonal graph over ordinary grid. Return additionally
  % a list of inner rectangles with singular points.  
   begin scalar u,dx,dy,xx,yy,l,w;
     dx:=float(xhi-xlo)/float(nx);
     dy:=float(yhi-ylo)/float(ny);
     l:=
     for j:=0:nx collect
     <<
       for i:=0:ny collect
       <<
        xx:=(xlo+i*dx); yy:=(ylo+j*dy);
        u:=plotevalform(ff,f,{x . xx,y . yy});
        if null u or eqcar(u,'overflow) 
         or numberp u and
             (zhi and u>zhi or zlo and u<zlo) then
        <<u:=nil;
          if 0<j and j<nx and 0<i and i<ny then
             w:={xx-dx,xx+dx,yy-dy,yy+dy}.w;
        >>;
        {t,t,xx,yy,u}
       >>
     >>;
     return w.l;
   end;

symbolic procedure ploteval3xy2 l;
    ploteval3xy3 {l};

symbolic procedure ploteval3xy3 ll;
  % Decompose ll into maximal regular grids.
 begin scalar l,lr,c,w,y,r,nw;
       integer n,m;
   while ll do
   <<l:=car ll; ll:=cdr ll;
       % find stripe with maximal lower left rectangle.
     w:={car l,cadr l}; l:=cdr l; 
     n:=min(ploteval3xybkpt(car w,0,nil),     % hor. only
            ploteval3xybkpt(cadr w,0,t));     %  hor. and vert
     c := t;
     while c and l and cdr l do
     <<m:=ploteval3xybkpt(cadr l,0,t);
       if izerop n and izerop m or n #>0 and not(n #> m) then
          <<l:= cdr l; w:=nconc(w,{car l})>>
       else c:=nil;
     >>;
     if cdr l then ll:=l . ll;
       % cut off subnet
     l:=nil; r:=nil; nw:=nil;
     for each row in w do
     <<if izerop n then row := cdr row
        else
       r:=(for i:=1:n collect <<y:=cddar row; row:=cdr row; y>>).r;
       if row then nw:=row.nw;
     >>;
     nw:= reversip nw;
        %print n;print{"streifen",length w,cddr caar w,
        %      cddr lastcar lastcar w};
        %print "gut";print r;print "rest";print nw;
        %if yesp "kill" then rederr "schluss";
     if nw then ll:=nw.ll;
     if r and cdr r then lr:=r.lr;
    >>;
    return lr;
  end;
          
symbolic procedure ploteval3xybkpt(w,n,m);
  % m=t: look for horizontal and vertical interrupts. Otherwise
  % look only for horizontal interrupts.
   if null w then n else
   if nil memq cddar w then n    % undefined point
    else
   if null cadar w               % x - break.
        or (m and null caar w)   % y - break
      then n+1 else
   ploteval3xybkpt(cdr w,n#+1,m);

endmodule;

module plotimp2; % Implicit plot: compute the varity {x,y|f(x,y)=c}.

% Author: Herbert Melenk, ZIB Berlin.

% data structure: 
%
%    point = {x,y,f(x,y),t1,t2,t3,...}
%               where ti are the numbers of the triangles which
%               have this point as a node. 
%               The point list is unique - adjacent triangles share
%               the list for equal nodes. The node numbers are
%               updated in place.
%
%    triangle = {t#,p1,p2,p3}
%       triangles are stored in triangle vector by number
%

fluid '(imp2!-triacount!* imp2!-trias!* !*imp2!-trace);

imp2!-triacount!*:=0;

symbolic procedure ploteval2xyimpl(rx,ry,f,x,y);
  begin scalar ll,l,form,g;
   for each c in plot!-contour!* do
   << form := plot!-form!-prep({'difference,cdr f,c},x,y);
      l:=imp2!-plot(car rx,cadr rx, car ry,cadr ry,
                   plot!-points(nil),form);
      ll:=append(ll,l);
    >>;
    if !*show_grid and null cdr plot!-contour!* 
         then g:= imp2!-show!-meshes();
    plotdriver(plot!-2imp,x,y,ll,g,car rx,cadr rx,car ry,cadr ry);   
  end;           

symbolic procedure imp2!-init();
  << imp2!-finit();
     if null imp2!-trias!* then imp2!-trias!* :=mkxvect()>>;

symbolic procedure imp2!-finit();
  <<if imp2!-trias!* then
      for i:=0:imp2!-triacount!* do xputv(imp2!-trias!*,i,nil);
    imp2!-triacount!*:=0;
  >>;

symbolic procedure mk!-point(x0,y0,fcn); 
         {x0,y0,apply2(fcn,x0,y0)};
 
symbolic procedure imp2!-delete!-pt!-reference(i,p);
    cdr cddr p := deletip(i,cdddr p);
 
symbolic procedure mk!-tria(i,p1,p2,p3);
  % make a triangle from 3 points. If the number is given,
  % reuse it. Otherwise generate a fresh number.
  begin scalar p; integer i;
    i := i or (imp2!-triacount!* := imp2!-triacount!* #+1);
    p:={i,p1,p2,p3,imp2!-tria!-zerop!*(p1,p2,p3)};
    xputv(imp2!-trias!*,i,p);
    aconc(p1,i); aconc(p2,i); aconc(p3,i);
      if !*imp2!-trace then print!-tria("new triangle ",p);
    return p;
  end;

symbolic procedure print!-tria(tx,w);
  <<prin2 tx; prin2 car w; w:=cdr w;
    prin2l {{car car w,cadr car w,{caddr car w}},
            {car cadr w,cadr cadr w,{caddr cadr w}},
            {car caddr w,cadr caddr w,{caddr caddr w}}};
      terpri();
  >>;

symbolic procedure imp2!-tria!-zerop!*(p1,p2,p3);
 % Here I test whether the triangle contains a zero line.
  begin scalar f1,f2,f3;
    f1 := caddr p1; 
    f2 := caddr p2;
    f3 := caddr p3;
    return f1*f2 <= 0.0 or f1*f3 <= 0.0;
  end;
    
symbolic procedure imp2!-tria!-zerop(w);
 % Fast access to stored zerop property.
     cadddr cdr w;

symbolic procedure imp2!-neighbours p;
 % Compute the direct neighbours of p - the traingles which have
 % two nodes respectively one complete edge in common with p.
  begin scalar l,p1,p2,p3; integer n;
   if fixp p then p:=xgetv(imp2!-trias!*,p);
   n:= car p; p:=cdr p;
   p1:=delete(n,cdddr car p);
   p2:=delete(n,cdddr cadr p);
   p3:=delete(n,cdddr caddr p);
   l:={find!-one!-common(p1,p2),
       find!-one!-common(p2,p3),
       find!-one!-common(p3,p1)};
   while nil memq l do l:=deletip(nil,l);
   return for each w in l collect xgetv(imp2!-trias!*,w);
  end;

symbolic procedure imp2!-edge!-length(p1,p2);
  begin scalar dx,dy;
    fdeclare('imp2!-edge!-length,dx,dy);
    dx:=thefloat car p1 - thefloat car p2;
    dy:=thefloat cadr p1 - thefloat cadr p2;
    return sqrt(dx*dx + dy*dy)
  end;

symbolic procedure imp2!-tria!-surface w;
 begin scalar x1,x2,x3,y1,y2,y3,p1,p2,p3;
   fdeclare('imp2!-tria!-surface,x1,x2,x3,y1,y2,y3);
   w:=cdr w;
   x1:=car (p1:=car w); y1:=cadr p1;
   x2:=car (p2:=cadr w); y2:=cadr p2;
   x3:=car (p3:=caddr w); y3:=cadr p3;
   return abs ((0.5*(x1*(y2-y3) + x2*(y3-y1) + x3*(y1-y2))));
 end;
    
symbolic procedure imp2!-tria!-length w;
 begin scalar p1,p2,p3;
   w:=cdr w;
   p1:=car w; p2:=cadr w; p3:=caddr w;
   return 
    (0.3*(imp2!-edge!-length(p1,p2) 
            + imp2!-edge!-length(p2,p3)
              + imp2!-edge!-length(p3,p1)));
 end;

symbolic procedure imp2!-tria!-midpoint(w);
   <<w:= cdr w;
     {(thefloat car car w 
        + thefloat car cadr w 
           + thefloat car caddr w)/3.0, 
      (thefloat cadr car w 
        + thefloat cadr cadr w 
           + thefloat cadr caddr w)/3.0}
   >>;

symbolic procedure imp2!-tria!-goodpoint(w,fn);
 % Here I assume that there is a zero in the triangle. Compute
 % a point inside the triangle which is as close as possible
 % to the desired line, but without local recomputation of
 % function values.
 begin scalar p1,p2,p3,f1,f2,f3,s1,s2,s3;
   w:=cdr w;
   p1:=car w; p2:=cadr w; p3:=caddr w;
   if (f1:=caddr p1)=0.0 then return {car p1,cadr p1} else
   if (f2:=caddr p2)=0.0 then return {car p2,cadr p2} else
   if (f3:=caddr p3)=0.0 then return {car p3,cadr p3};
   s1:=f1<0.0; s2:=f2<0.0; s3:=f3<0.0;
   return if s1=s2 then 
       imp2!-tria!-goodpoint1(p1,f1,p3,f3,p2,f2,fn)
   else if s1=s3 then
       imp2!-tria!-goodpoint1(p1,f1,p2,f2,p3,f3,fn)
   else
       imp2!-tria!-goodpoint1(p2,f2,p1,f1,p3,f3,fn)
 end;

%symbolic procedure imp2!-tria!-goodpoint1(p1,f1,p2,f2,p3,f3,fn);
% % Now I know that f2 has the opposite sign to f1 and f3.
% % I take the linearly interpolated zero of f on p1-p2 and p2-p3
% % return their arithmetic mean value which lies inside the 
% % triangle.
%  begin scalar x1,x2,y1,y2,s;
%    fdeclare (x1,x2,y1,y2,s,f1,f2,f3);
%    s:=f1-f2; 
%    x1:=(f1*thefloat car p2  - f2*thefloat car p1)/s;
%    y1:=(f1*thefloat cadr p2 - f2*thefloat cadr p1)/s;
%    s:=f3-f2; 
%    x2:=(f3*thefloat car p2  - f2*thefloat car p3)/s;
%    y2:=(f3*thefloat cadr p2 - f2*thefloat cadr p3)/s;
%    return {(x1+x2)*0.5, (y1+y2)*0.5};
%  end;

symbolic procedure imp2!-tria!-goodpoint1(p1,f1,p2,f2,p3,f3,fn);
 % Now I know that f2 has the opposite sign to f1 and f3.
 % F1 and f3 are non-zero.
 % I use the midpoint of the p1-p3 edge and look for an
 % approximation of a zero on the connection of the midpoint 
 % and p2 using regula falsi.
  begin scalar x1,x2,y1,y2,x3,y3,s;
    fdeclare (x1,x2,x3,y1,y2,y3,s,f1,f2,f3);
    f1:=(f1+f3)*0.5; 
    x1:=(thefloat car p1  + thefloat car p3)*0.5;
    y1:=(thefloat cadr p1  + thefloat cadr p3)*0.5;
    x2:=car p2; y2:=cadr p2;
    s:=f2-f1;
    x3:=(x1*f2 - x2*f1)/s;
    y3:=(y1*f2 - y2*f1)/s;
    f3:=apply2(fn,x3,y3);
    if f2*f3>=0 then
    <<s:=f1-f3; x3:=(x3*f1-x1*f3)/s; y3:=(y3*f1-y1*f3)/s>>
      else
    <<s:=f2-f3; x3:=(x3*f2-x2*f3)/s; y3:=(y3*f2-y2*f3)/s>>;
 done:
    return{x3,y3};
  end;

symbolic procedure imp2!-tri!-refine!-one!-tria (w,fn);
 % Refine one triangle by inserting a new point in the mid
 % of the longest edge. Also, refine the triangle adjacent
 % to that edge with the same point.
  begin scalar p1,p2,p3,pn,xn,yn,new,nb,y; integer i;
    scalar x1,x2,y1,y2,d1,d2,d3,s;
    fdeclare (x1,x2,y1,y2,s,d1,d2,d3);
   if fixp w then w :=xgetv(imp2!-trias!*,w); % record.
     if !*imp2!-trace then print!-tria("refine ",w);
   i:=car w; w :=cdr w;

     % delete reference to this triangle.
   imp2!-delete!-pt!-reference(i,car w);
   imp2!-delete!-pt!-reference(i,cadr w);
   imp2!-delete!-pt!-reference(i,caddr w);
  
     % find longest edge
   d1:=imp2!-edge!-length(car w,cadr w);
   d2:=imp2!-edge!-length(cadr w,caddr w);
   d3:=imp2!-edge!-length(car w,caddr w);
   if d1>=d2 and d1>=d3 then 
     <<p1:=car w; p2:=cadr w; p3:=caddr w>>
   else if d2>=d1 and d2>=d3 then 
     <<p1:=cadr w; p2:=caddr w; p3:=car w>>
   else 
     <<p1:=car w; p2:=caddr w, p3:=cadr w>>;
    % identify the neighbour triangle.
   nb := find!-one!-common(cdddr p1,cdddr p2);
    % compute new point almost at the mid.
   s:=0.45+0.01*random(10);  
   x1:=car p1; y1:=cadr p1;
   x2:=car p2; y2:=cadr p2;
   xn:=x1*s+x2*(1.0-s);
   yn:=y1*s+y2*(1.0-s);
   pn:=mk!-point(xn,yn,fn);
construct:
    % construct new triangles
   new:=mk!-tria(i,p1,pn,p3).new;
   new:=mk!-tria(nil,p2,pn,p3).new;
    % update the neighbour, if there is one
   if null nb then return new;
   w:=nb; nb:=nil;
   if fixp w then w :=xgetv(imp2!-trias!*,w); % record.
   i:=car w; w:=cdr w;
   imp2!-delete!-pt!-reference(i,car w);
   imp2!-delete!-pt!-reference(i,cadr w);
   imp2!-delete!-pt!-reference(i,caddr w);
   
    % find the far point
   p3 := if not((y:=car w) eq p1 or y eq p2) then y else
         if not((y:=cadr w) eq p1 or y eq p2) then y else
         caddr w;
   goto construct;
  end;

%symbolic procedure imp2!-tri!-refine!-one!-tria!-center (w,fn);
% % Refine one triangle by inserting a new point in the center.
%  begin scalar p1,p2,p3,pn,xn,yn,new,nb,y; integer i;
%    scalar x1,x2,x3,y1,y2,y3;
%   fdeclare (x1,x2,x3,y1,y2,y3);
%   if fixp w then w :=xgetv(imp2!-trias!*,w); % record.
%     if !*imp2!-trace then print!-tria("refine ",w);
%   i:=car w; w :=cdr w;
%
%     % delete reference to this triangle.
%   imp2!-delete!-pt!-reference(i,car w);
%   imp2!-delete!-pt!-reference(i,cadr w);
%   imp2!-delete!-pt!-reference(i,caddr w);
%  
%    % compute center.
%   p1:=car w; p2:=cadr w; p3:=caddr w;
%   x1:=car p1; y1:=cadr p1;
%   x2:=car p2; y2:=cadr p2;
%   x3:=car p3; y3:=cadr p3;
%   xn:=(x1+x2+x3)*0.33333;
%   yn:=(y1+y2+y3)*0.33333;
%   pn:=mk!-point(xn,yn,fn);
%construct:
%    % construct new triangles
%   new:=mk!-tria(i,p1,p2,pn).new;
%   new:=mk!-tria(nil,p2,p3,pn).new;
%   new:=mk!-tria(nil,p1,p3,pn).new;
%   return new;
%  end;

symbolic procedure find!-one!-common(u,v);
  % fast equivalent to "car intersection(u,v)".
  if null u then nil else
  if memq(car u,v) then car u else
     find!-one!-common(cdr u,v);

%%%%%% Main program for implicit plot.

symbolic procedure imp2!-plot(xlo,xhi,ylo,yhi,pts,fcn);
  begin scalar p1,p2,p3,p4,sf,z;
   imp2!-init();
    % setup initial triangularization.
   p1:=mk!-point(xlo,ylo,fcn); 
   p2:=mk!-point(xhi,ylo,fcn); 
   p3:=mk!-point(xhi,yhi,fcn); 
   p4:=mk!-point(xlo,yhi,fcn); 
   mk!-tria(nil,p1,p2,p3);
   mk!-tria(nil,p1,p3,p4);
   sf:=((xhi-xlo)+(yhi-ylo))*1.5/float pts;
   z:=imp2!-plot!-refine(sf,fcn);
     if !*imp2!-trace then 
       for each w in z do print!-tria("zero triangle:",w);
if !*test_plot then print "collect";
   z:=imp2!-plot!-collect(z);
if !*test_plot then print "draw";
   z:=imp2!-plot!-draw(z,fcn);
   if not !*show_grid then imp2!-finit();
   return z;
 end;
   
symbolic procedure imp2!-plot!-refine(sflimit,fcn);
  begin scalar z,zn,w,s,limit2,again;
    integer i,j;
    limit2 := 0.15*sflimit;
   
  % In the first stage I refine all triangles until they
  % are fine enough for a coarse overview.
    again := t;

if !*test_plot then print "phase1";
phase1:
    z:=nil; again:=nil;
    for i:=imp2!-triacount!* step -1 until 1 do
    << w := xgetv(imp2!-trias!*,i);
       if imp2!-tria!-length w > sflimit then
           <<again:=t; imp2!-tri!-refine!-one!-tria (w,fcn)>>
        else if not again and imp2!-tria!-zerop w 
           then z:=car w.z;
    >>;
    j:=j #+ 1;
    if j+j #< plot!-refine!* or again then goto phase1;
 
 % Next I refine further only the triangles which contain a zero.
 % I must store in z only the numbers of triangles because these
 % may be modified as side effects by copying.

if !*test_plot then print "phase2";
phase2:
    for each w in z do
      if (s:=imp2!-tria!-length (w:=xgetv(imp2!-trias!*,w))) >limit2
	 then <<for each q in imp2!-tri!-refine!-one!-tria (w,fcn) do
		if imp2!-tria!-zerop q and not memq(car q,zn)
		  then zn:=car q.zn>>;
    z:=zn; zn:=nil;
    if z then goto phase2;

 % In the third phase I refine those critical triangles futher:
 % non-zeros with two zero neighbours. These may be turning points
 % or bifurcations.

if !*test_plot then print "phase3";
phase3:
    for i:=imp2!-triacount!* step -1 until 1 do
      imp2!-refine!-phase3(i,i,plot!-refine!*/2,fcn,limit2*0.3);

 % Return the final list of zeros in ascending order.
    z:=for i:=1:imp2!-triacount!* join
      if imp2!-tria!-zerop(w:=xgetv(imp2!-trias!*,i)) then {w};
   return z;
  end;

symbolic procedure imp2!-refine!-phase3(i,low,lev,fn,lth);
  begin scalar w; integer c; 
    if lev=0 then return nil;
    w:=if numberp i then xgetv(imp2!-trias!*,i) else i;
    if car w<low or imp2!-tria!-length w < lth then return nil;
    lev:=lev-1;
    for each q in imp2!-neighbours w do
      if imp2!-tria!-zerop q then c:=c+1;
    if imp2!-tria!-zerop w  
       then (if c eq 2 then return nil)
       else (if c #< 2 then return nil);
    for each q in imp2!-tri!-refine!-one!-tria (w,fn) do
        imp2!-refine!-phase3(q,low,lev,fn,lth);
  end;

symbolic procedure imp2!-plot!-collect(z);
   begin scalar lines,l,q,p,s;

    for each w1 in z do
       for each w2 in imp2!-neighbours w1 do
          if car w2 > car w1 and imp2!-tria!-zerop w2 then
            q:=(w1.w2) . q;
 
 lineloop:
     if null q then return lines;
     l:={caar q, (p:=cdar q)}; q:= cdr q;
      % first I extend the back end.
     while q and p do
     <<
       if(s:= atsoc(p,q)) then l:=nconc(l,{p:=cdr s}) else
       if(s:=rassoc(p,q)) then l:=nconc(l,{p:=car s});
       if s then q:=delete(s,q) else p:=nil;
     >>;

      % next I extend the front end.
     p:=car l;
     while q and p do 
     <<
       if(s:=rassoc(p,q)) then l:=(p:=car s).l else
       if(s:= atsoc(p,q)) then l:=(p:=cdr s).l;
       if s then q:=delete(s,q) else p:=nil;
     >>;

     lines := l.lines;
     goto lineloop;
   end;

symbolic procedure imp2!-plot!-draw(l,fn);
  begin scalar r,s,q;
  q:=for each w in l collect
   <<r:=nil; 
     for each q in w join
       <<s:=imp2!-tria!-goodpoint(q,fn);
         if r neq s then {r:=s}>>
   >>;
  return q;
  end;


symbolic procedure imp2!-show!-meshes();
 % generate plot of meshes;
  begin scalar d,l,w,s,p1,p2; integer i;
    i:=1;
 loop:
    w:=xgetv(imp2!-trias!*,i);
    if null w then 
    <<imp2!-finit(); return l>>;
    w:=cdr w;
    d:= {{car car w, cadr car w},
         {car cadr w,cadr cadr w},
         {car caddr w,cadr caddr w},
         {car car w, cadr car w}} ;
    while d and cdr d do
    <<p1:=car d; p2:=cadr d; d:=cdr d;
      if car p1 > car p2 then <<w:=p2;p2:=p1;p1:=w>>;
      s:={p1,p2};
      if not member(s,l) then l:=s.l
    >>;
    i:=i+1;
    goto loop;
   end;
              
endmodule; % plotimp2

module plotimp3; % Implicit plot: compute the varity {x,y,z|f(x,y,z)=0}.

% Author: Herbert Melenk, ZIB Berlin.

% data structure: cubes. 

symbolic procedure ploteval3impl(x,y,z);
  begin scalar rx,ry,rz,f,fcn;     
     rx:=plotrange(x,
      reval(plot_xrange or '(!*interval!* -10 10)));
     ry:=plotrange(y,
      reval(plot_yrange or '(!*interval!* -10 10)));
     rz:=plotrange(z, 
      reval(plot_zrange or '(!*interval!* -10 10)));
     fcn := car reverse plotfunctions!*;
     f:= ploteval3impl1(cdar plotfunctions!*,
         x,car rx,cadr rx,
         y,car ry,cadr ry,
         z,car rz,cadr rz);
     plotdriver(plot!-3exp!-reg,x,y,z,f);
  end;

symbolic procedure ploteval3impl1(f,x,xlo,xhi,y,ylo,yhi,z,zlo,zhi);
   begin scalar u,dx,dy,dz,xx,yy,zz,l,ff,pts,val,w,q,qq,pt,found,done;
        integer nx,ny,nz;
     ff := rdwrap f;
     xlo:=rdwrap xlo; xhi:=rdwrap xhi;
     ylo:=rdwrap ylo; yhi:=rdwrap yhi;
     dx:=float(xhi-xlo)/float(nx:=plot!-points(x));
     dy:=float(yhi-ylo)/float(ny:=plot!-points(y));
     dz:=float(zhi-zlo)/float(nz:=plot!-points(z));
     pts := mk!-p!-array3(nx,ny,nz);
     val:= mk!-p!-array3(nx,ny,nz);

        % Step 1: compute a complete grid in 3d.
     for i:=0:nx do
     <<
       xx:=(xlo+i*dx);
       for j:=0:ny do
       <<
        yy:=(ylo+j*dy);
        for k:=0:nz do
        <<
          zz:=(zlo+k*dz);
          p!-put3(pts,i,j,k,{xx,yy,zz});
          u:=plotevalform(ff,f,{x . xx,y . yy,z . zz});
          if eqcar(u,'overflow) then u:=1.0;
          p!-put3(val,i,j,k,u);
        >>;
       >>
     >>;
 
        % Step 2: extract zero points.
   done := t;
   while done do
   <<done := nil; w:=
     for i:=0:nx #-1 collect
      for j:=0:ny #-1 collect
      <<q:=nil; found:=nil;
       pt := p!-get3(pts,i,j,1);
       for k:=nz step -1 until 0 do
        if null found then
        <<if null q then q:=p!-get3(val,i,j,k);
          qq:=p!-get3(val,i,j,k); 
          if q and qq and q*qq <= 0.0 then
            found := if q=0.0 then caddr p!-get3(pts,i,j,k)
                else ploteval3impl3(caddr p!-get3(pts,i,j,k),qq,
                                    caddr p!-get3(pts,i,j,k#+1),q);
          if q=0.0 or qq=0.0 or not found then p!-put3(val,i,j,k,nil);
          done:=done or found;
          q:=qq;
        >>;
        {t,t,car pt,cadr pt,found}
       >>;
     if done then l:=w.l;
    >>;
    return ploteval3xy3 l;
   end;

symbolic procedure ploteval3impl3(p1,f1,p2,f2);
  % linear interpolation
  <<
    fdeclare(f1,f2,p1,p2);
    (f1*p2 - f2*p1)/(f1-f2)>>; 

endmodule;

module plotnum; % Numeric evaluation of algebraic expressions.

fluid '(plotsynerr!* ploteval!-alist2!*);

global '(!*plotinterrupts);

flag('(plus plus2 difference times times2 quotient exp expt expt!-int
       minus sin cos tan cot asin acos acot atan log),'plotevallisp);

symbolic procedure plotrounded d;
 % Switching rounded mode safely on and off.
  begin scalar o,!*protfg,c,r,!*msg;
   !*protfg:=t;
   if null d then
    <<c:=!*complex; r:=!*rounded;
      if c then <<setdmode('complex,nil); !*complex := nil>>;
      if not r and dmode!* then
          <<o:=get(dmode!*,'dname); setdmode(o,nil)>>;
      o:={o,r,!*roundbf,c,precision 0};
      if not r then <<!*rounded:=t; setdmode('rounded,t)>>;
      !*roundbf:=nil;
      if c then <<setdmode('complex,t); !*complex := t>>;
      return o;
    >>
   else
    <<
       % reconstruct the previous configuration.
      if !*complex then setdmode('complex,nil);
      if null(!*rounded:=cadr d) then setdmode('rounded,nil);
      !*roundbf:=caddr d; 
      if car(d) then setdmode(car d,t);
      if !*complex then setdmode('complex,t);
      precision car cddddr d;
    >>;
   end;

symbolic procedure adomainp u;
 % numberp test in an algebraic form.
   numberp u or (pairp u and idp car u and get(car u,'dname))
             or eqcar(u,'minus) and adomainp cadr u;

symbolic procedure revalnuminterval(u,num);
 % Evaluate u as interval; numeric bounds required if num=T.
  begin scalar l;
    if not eqcar(u,'!*interval!*) then typerr(u,"interval");
    l:={reval cadr u,reval caddr u};
    if null num or(adomainp car l and adomainp cadr l)then return l;
    typerr(u,"numeric interval");
  end;

ploteval!-alist2!*:={{'x . 1},{'y . 2}};

symbolic procedure plot!-form!-prep(f,x,y);
  % f is a REDUCE function expression depending of x and y. 
  % Compute a form which allows me to evaluate "f(x,y)" as
  % a LISP expr.
  {'lambda,'(!&1 !&2),
     {'plot!-form!-call2,mkquote rdwrap f,mkquote f,
                mkquote x,'!&1,
                mkquote y,'!&2}};

symbolic procedure plot!-form!-call2(ff,f,x,x0,y,y0);
 % Here I hack into the statically allocated a-list in order
 % to save cons cells.
  begin scalar a;
     a:=car ploteval!-alist2!*;
     car a := x; cdr a := x0;
     a:=cadr ploteval!-alist2!*;   
     car a:= y; cdr a:= y0;
     return plotevalform(ff,f,ploteval!-alist2!*);
  end;

% The following routines are used to transform a REDUCE algebraic
% expression into a form which can be evaluated directly in LISP.

symbolic procedure rdwrap f;
  begin scalar r,!*msg,!*protfg;
    !*protfg:=t;
    r:=errorset({'rdwrap1,mkquote f},nil,nil);
    return if errorp r then 'failed else car r;
  end;

symbolic procedure rdwrap1 f;
if numberp f then float f
  else if f='pi then 3.141592653589793238462643
  else if f='e then 2.7182818284590452353602987
  else if f='i and !*complex then rederr "i not LISP evaluable" 
  else if atom f then f
  else if get(car f,'dname) then rdwrap!-dm f
  else if eqcar(f, 'MINUS) then
    begin scalar x;
       x := rdwrap1 cadr f;
       return if numberp x then minus float x else {'MINUS, x}
    end
  else if eqcar(f,'expt) then rdwrap!-expt f
  else
   begin scalar a,w;
     if null getd car f or not flagp(car f,'plotevallisp)
                then typerr(car f,"Lisp arithmetic function (warning)");
     a:=for each c in cdr f collect rdwrap1 c;
     if (w:=atsoc(car f,'((plus.plus2)(times.times2)))) 
             then return rdwrapn(cdr w,a);
     return car f . a;
   end;

symbolic procedure rdwrapn(f,a);
  % Unfold a n-ary arithmetic call.
   if null cdr a then car a else {f,car a,rdwrapn(f,cdr a)};

symbolic procedure rdwrap!-dm f;
 % f is a domain element.
  if car f eq '!:RD!: then 
          if atom cdr f then cdr f else bf2flr f
  else if car f eq '!:CR!: then rdwrap!-cr f
  else if car f eq '!:GI!:
   then rdwrap!-cmplex(f,float cadr f,float cddr f)
  else if car f eq '!:DN!: then rdwrap2 cdr f
  else << plotsynerr!*:=t;
       rerror(PLOTPACKAGE, 32, {get(car f, 'DNAME),
                                "illegal domain for PLOT"})
    >>;
  
symbolic procedure rdwrap!-cr f;
  begin scalar re,im;
    re := cadr f;
    if not atom re then re := bf2flr re;
    im := cddr f;
    if not atom im then im := bf2flr im;
    return rdwrap!-cmplx(f,re,im);
  end;

symbolic procedure rdwrap!-cmplx(f,re,im);
   if abs im * 1000.0 > abs re then typerr(f,"real number") else re;

symbolic procedure plotrepart u;
   if car u eq '!:rd!: then u else
   if car u memq '(!:gi!: !:cr!:) then '!:rd!: . cadr u;

symbolic procedure rdwrap!-expt f;
  % preserve integer second argument.
  if fixp caddr f then 
     if caddr f>0 then {'expt!-int,rdwrap1 cadr f,caddr f}
       else {'quotient,1,{'expt!-int,rdwrap1 cadr f,-caddr f}}
    else {'expt,rdwrap1 cadr f, rdwrap caddr f};

symbolic procedure rdwrap2 f;
% Convert from domain to LISP evaluable value.
  if atom f then f else float car f * 10^cdr f;

symbolic procedure rdwrap!* f;
  % convert a domain element to float.
  if null f then 0.0 else rdwrap f;

symbolic procedure rdunwrap f; 
    if f=0.0 then 0 else if f=1.0 then 1 else '!:rd!: . f;

symbolic procedure expt!-int(a,b); expt(a,fix b);

symbolic procedure plotevalform(ff,f,a);
  % ff is LISP evaluable,f REDUCE evaluable.
  begin scalar u,w,!*protfg,mn,r,!*msg;
    !*protfg := t;
    if ff='failed then goto non_lisp;
    u:= errorset({'plotevalform1,mkquote ff,mkquote a},nil,nil);
    if not ploterrorp u and numberp (u:=car u) then 
    <<if abs u > plotmax!* then return plotoverflow plotmax!* else
      return u;
    >>;
 non_lisp:
    w := {'simp, mkquote
     sublis(
       for each p in a collect (car p.rdunwrap cdr p),
       f)};
    u := errorset(w,nil,nil);
    if ploterrorp u or (u:=car u) eq 'failed then return nil;
    if denr u neq 1 then return nil;
    u:=numr u;
    if null u then return 0; % Wn
    if numberp u then <<r:=float u; goto x>>;
    if not domainp u or not memq(car u,'(!:rd!: !:gi!: !:cr!:)) 
          then return nil;
    if evalgreaterp(plotrepart u, fl2rd plotmax!*) then
     return plotoverflow plotmax!* else
    if evalgreaterp(fl2rd (-plotmax!*),plotrepart u) then
     return plotoverflow (-plotmax!*);
    r:=errorset({'rdwrap!-dm,mkquote u},nil,nil);
    if errorp r or(r:=car r) eq 'failed then return nil;
  x: return if mn then -r else r;
  end;

symbolic procedure plotoverflow u;
   <<if not !*plotoverflow then 
      lprim "plot number range exceeded";
     !*plotoverflow := t;
     'overflow . u
   >> where !*protfg = nil;

symbolic procedure plotevalform0(f,a);
  (if ploterrorp u 
       then typerr(f,"expression for plot") else car u)
     where u=
   errorset({'plotevalform1,mkquote f,mkquote a},nil,nil);

%symbolic procedure plotevalform1(f,a);
% begin scalar x,w;
%  if numberp f then return float f;
%  if (x:=assoc(f,a)) then return plotevalform1(cdr x,a);
%  if not atom f and flagp(car f,'plotevallisp) then
%    return eval 
%      (car f . for each y in cdr f collect plotevalform1(y,a));
%  if not atom f then f :=
%    car f . for each y in cdr f collect prepf plotevalform2(y,a);
%  w:=simp f;
%  if denr w neq 1 or not domainp numr w then goto err;
%  return rdwrap!* numr w;
% err:
%  plotsynerr!*:=t;
%  typerr(prepsq w,"numeric value");
%end;

symbolic procedure plotevalform1(f,a);
 begin scalar x;
  if numberp f then return float f;
  if (x:=assoc(f,a)) then return plotevalform1(cdr x,a);
  if atom f then go to err;
  return if cddr f then 
     idapply(car f,{plotevalform1(cadr f,a),plotevalform1(caddr f,a)})
    else
     idapply(car f,{plotevalform1(cadr f,a)}); 
 err:
  plotsynerr!*:=t;
  typerr(prepsq f,"numeric value");
end;


%symbolic procedure plotevalform2(f,a);
% begin scalar x,w;
%  if fixp f then return f;
%  if floatp f then return rdunwrap f;
%  if (x:=assoc(f,a)) then return plotevalform2(cdr x,a);
%  if not atom f and flagp(car f,'plotevallisp) then
%    return rdunwrap eval 
%      (car f . for each y in cdr f collect plotevalform1(y,a));
%  if not atom f then f :=
%    car f . for each y in cdr f collect prepf plotevalform2(y,a);
%  w:=simp f;
%  if denr w neq 1 or not domainp numr w then goto err;
%  return numr w;
% err:
%  plotsynerr!*:=t;
%  typerr(prepsq w,"numeric value");
%end;

symbolic procedure ploterrorp u;
   if u member !*plotinterrupts 
      then rederr prin2 "***** plot killed"
    else atom u or cdr u;

endmodule;

module parray; % multidimensional arrays.

symbolic procedure mk!-p!-array3(nx,ny,nz);
  <<for i:=0:nx do iputv(w,i,mk!-p!-array2(ny,nz)); w>>
       where w=mkvect(nx#+1);

symbolic procedure mk!-p!-array2(ny,nz);
  <<for i:=0:ny do iputv(w,i,mkvect(nz#+1)); w>>
       where w=mkvect(ny#+1);

symbolic procedure p!-get3(v,i,j,k);
  igetv(igetv(igetv(v,i),j),k);

symbolic procedure p!-put3(v,i,j,k,w);
  iputv(igetv(igetv(v,i),j),k,w);

endmodule;

module xvect; % Support for vectors with adaptive length.

% Author: Herbert Melenk, ZIB-Berlin.

symbolic procedure mkxvect(); {mkvect(128)};

symbolic procedure xputv(v,n,w);
  begin scalar i,j;
    i:=iquotient(n,128); j:=iremainder(n,128);
    while(i>0) do
    <<if null cdr v then cdr v:= mkxvect();
      v:=cdr v; i:=i #- 1>>;
    iputv(car v,j,w);
    return w;
  end;

symbolic procedure xgetv(v,n);
  begin scalar i,j,w;
    i:=iquotient(n,128); j:=iremainder(n,128);
    while(i>0 and v) do
    <<v:=cdr v; i:=i #- 1>>;
    w:=if v then igetv(car v,j);
    return w
  end;

endmodule;

end;
//E*O*F plot/plot.red//

exit 0