@@ -1,1650 +1,1650 @@
-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
+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