reduce-historical/r36/cslsrc/gnuplot.red

module plot; % device and driver independent plot services.

% Author: Herbert Melenk.
%    Adjusted by A C Norman to be compatible with CSL - the original
%    was written to be fairly PSL-specific.

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!*          %  definition of y-bigstep for smooth

      plotoptions!*        % list for collecting the options.

);

fluid '(

      plotderiv!*          %  derivative for 2d plot

);


!#if (member 'psl lispsystem!*)

macro procedure !@if u;
   <<if eval cadr u then caddr u else cadddr u>>;

!@if(errorp errorset('!*writingfaslfile,nil,nil)
         or not !*writingfaslfile
         or errorp errorset('(load fcomp),nil,nil),
       <<prin2t "no support for fast float!";
% ACN believes that the following line is just WRONG in that the
% parentheses written in (x) are ignored by the RLISP parser, to the detriment
% of the way that DM will handle things!
         dm(fdeclare,(x),nil); dm(thefloat,(x),cadr x);>>,
       nil);

!#else

symbolic macro procedure fdeclare u; nil;
symbolic macro procedure thefloat u; cadr u;

!#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,'! !.!.! );
>>;

!*msg := t;

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

global'(!*plotrefine);

fluid '(plotprecision!*);

plotprecision!* := 0.9995;

switch  plotrefine;        % if OFF no refinments are computed
                           % for 2D plot.

!*plotrefine:=t;

fluid '(!*show_grid);

switch show_grid;

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;
    if null plotdriver!* then
      rederr "no active device driver for PLOT";
    m:=plotrounded(nil);
    !*plotoverflow := nil;
    plotranges!* := plotfunctions!* := nil;
    plotstyle!* := 'lines;
    plotdriver(init);
    for each option in u do ploteval1 reval option;
    ploteval2();
    plotrounded(m);
  end;

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 eqcar(option,'equal) and (x:=get(cadr  option,do))
             then apply(x,list caddr option)
       else ploteval0 option;
   end;

symbolic procedure ploteval0 option;
  begin scalar l,r,opt;
    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!*);
        typerr(option,"plot option")
      >>;
    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>>;
    l:=reval cadr option;
    r:=reval caddr option;
    if l=0 or r=0 then return
      plotfunctions!*:=('implicit.if l=0 then r else l)
                       . plotfunctions!*;
    if not idp l then rederr "illegal option in PLOT";
    if l memq '(size terminal view) then
      <<plotoptions!*:=(l.r).plotoptions!*; return>>;
    if eqcar(r,'!*interval!*) then
      % must be a range
     <<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;
   for each u in plotfunctions!* do
     <<if car u eq 'implicit 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:=car dvar
   >>;
   if length ivars=1 then ploteval2x(car ivars,dvar) else
   if length ivars=2 then ploteval3xy(car ivars,cadr ivars,dvar)
    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 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 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 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) 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);

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!*;
     % compute algebraic derivative.
    plotderiv!*:=reval {'df,f,x};
    plotderiv!*:= if smemq('df,plotderiv!*) then nil
        else rdwrap plotderiv!* . plotderiv!*;
    ff:=rdwrap f;
    plot_xmesh := reval plot_xmesh;
    p:=float plot_xmesh;
    d:=(d0:=(xhi-xlo))/p;
    v:=xlo;
    for i:=0:plot_xmesh do
     <<vv:=if i=0 or i=plot_xmesh 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;
     >>;
     variation!* :=
     if yhi and not(yhi=plotmax!*) then (yhi-ylo) else
        ploteval2xvariation l;

    if !*plotrefine 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,8)
   >>;
   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 plotexp3; % Computing surface z=f(x,y) with regular grid.

% Author: Herbert Melenk, ZIB Berlin.

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,dx,dy,xx,yy,l,ff;
     ff := rdwrap f;
     xlo:=rdwrap xlo; xhi:=rdwrap xhi;
     ylo:=rdwrap ylo; yhi:=rdwrap yhi;
     plot_xmesh := reval plot_xmesh;
     plot_ymesh := reval plot_ymesh;
     dx:=float(xhi-xlo)/float(plot_xmesh);
     dy:=float(yhi-ylo)/float(plot_ymesh);
     l:=
     for j:=0:plot_ymesh collect
     <<
       for i:=0:plot_xmesh collect
       <<
        xx:=(xlo+i*dx); yy:=(ylo+j*dy);
        u:=plotevalform(ff,f,{x . xx,y . yy});
        if eqcar(u,'overflow) then u:=nil;
        if null u then  % look for an isolated singularity
          u:=ploteval3xysingular(ff,f,x,xx,dx,y,yy,dy,zhi,zlo);
        u := if zhi and u>zhi then zhi else
             if zlo and u<zlo then zlo else u;
        {xx,yy,u}
       >>
     >>;
     return l;
   end;

symbolic procedure ploteval3xysingular(ff,f,x,xx,dx,y,yy,dy,zhi,zlo);
 % set up an iteration approaching a cirital point.
   <<dx:=dx/4; dy:=dy/4;
     ploteval3xysingular1(ff,f,x,xx,dx,y,yy,dy,zhi,zlo,
       plotevalform(ff,f,{x . (xx+dx), y . (yy+dy)}),0)
    >>;

symbolic procedure ploteval3xysingular1(ff,f,x,xx,dx,y,yy,dy,zhi,zlo,w,c);
  if null w then 0 else if c>8 then w else
  if w>zhi then zhi else
  if w<zlo then zlo else
  begin scalar wnew;
    dx:=dx/2; dy:=dy/2;
    wnew := plotevalform(ff,f,{x . (xx+dx), y . (yy+dy)});
    return
     if null wnew then w else
     if abs(wnew-w) <abs wnew/20 then wnew else
       ploteval3xysingular1(ff,f,x,xx,dx,y,yy,dy,zhi,zlo,wnew,c+1);
  end;

endmodule;

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

% 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 '(p3!-triacount!* p3!-trias!* !*p3!-trace);

symbolic procedure ploteval2xyimpl(rx,ry,f,x,y);
  begin scalar l,form,g;
    form := plot!-form!-prep(cdr f,x,y);
    l:=p3!-plot(car rx,cadr rx, car ry,cadr ry,reval plot_xmesh,form);
    if !*show_grid then g:= p3!-show!-meshes();
    plotdriver(plot!-2imp,x,y,l,g,car rx,cadr rx,car ry,cadr ry);
  end;

symbolic procedure p3!-init();
   <<p3!-triacount!*:=0;
     p3!-trias!*:=mkxvect();
   >>;

symbolic procedure mk!-point(x0,y0,fcn);
         {x0,y0,apply2(fcn,x0,y0)};

symbolic procedure p3!-delete!-pt!-reference(i,p);
    cdr cddr p := delete(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 (p3!-triacount!* := p3!-triacount!* #+1);
    p:={i,p1,p2,p3,p3!-tria!-zerop!*(p1,p2,p3)};
    xputv(p3!-trias!*,i,p);
    aconc(p1,i); aconc(p2,i); aconc(p3,i);
      if !*p3!-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 p3!-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 p3!-tria!-zerop(w);
 % Fast access to stored zerop property.
     cadddr cdr w;

symbolic procedure p3!-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(p3!-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:=delete(nil,l);
   return for each w in l collect xgetv(p3!-trias!*,w);
  end;

symbolic procedure p3!-edge!-length(p1,p2);
  begin scalar dx,dy;
    fdeclare('p3!-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 p3!-tria!-surface w;
 begin scalar x1,x2,x3,y1,y2,y3,p1,p2,p3;
   fdeclare('p3!-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 p3!-tria!-length w;
 begin scalar p1,p2,p3;
   w:=cdr w;
   p1:=car w; p2:=cadr w; p3:=caddr w;
   return
    (0.3*(p3!-edge!-length(p1,p2)
            + p3!-edge!-length(p2,p3)
              + p3!-edge!-length(p3,p1)));
 end;

symbolic procedure p3!-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 p3!-tria!-goodpoint(w);
 % 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
       p3!-tria!-goodpoint1(p1,f1,p3,f3,p2,f2)
   else if s1=s3 then
       p3!-tria!-goodpoint1(p1,f1,p2,f2,p3,f3)
   else
       p3!-tria!-goodpoint1(p2,f2,p1,f1,p3,f3)
 end;

symbolic procedure p3!-tria!-goodpoint1(p1,f1,p2,f2,p3,f3);
 % 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 p3!-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(p3!-trias!*,w); % record.
     if !*p3!-trace then print!-tria("refine ",w);
   i:=car w; w :=cdr w;

     % delete reference to this triangle.
   p3!-delete!-pt!-reference(i,car w);
   p3!-delete!-pt!-reference(i,cadr w);
   p3!-delete!-pt!-reference(i,caddr w);

     % find longest edge
   d1:=p3!-edge!-length(car w,cadr w);
   d2:=p3!-edge!-length(cadr w,caddr w);
   d3:=p3!-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(p3!-trias!*,w); % record.
   i:=car w; w:=cdr w;
   p3!-delete!-pt!-reference(i,car w);
   p3!-delete!-pt!-reference(i,cadr w);
   p3!-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 p3!-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(p3!-trias!*,w); % record.
%     if !*p3!-trace then print!-tria("refine ",w);
%   i:=car w; w :=cdr w;
%
%     % delete reference to this triangle.
%   p3!-delete!-pt!-reference(i,car w);
%   p3!-delete!-pt!-reference(i,cadr w);
%   p3!-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 p3!-plot(xlo,xhi,ylo,yhi,pts,fcn);
  begin scalar p1,p2,p3,p4,sf,z;
   p3!-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:=p3!-plot!-refine(sf,fcn);
     if !*p3!-trace then
       for each w in z do print!-tria("zero triangle:",w);
   z:=p3!-plot!-collect(z);
   return p3!-plot!-draw(z);
   return nil;
 end;

symbolic procedure p3!-plot!-refine(sflimit,fcn);
  begin scalar z,zn,w,s,limit2,again;
    integer i;
    limit2 := 0.15*sflimit;

  % In the first stage I refine all triangles until they
  % are fine enough for a coarse overview.
    again := t;
phase1:
    z:=nil; again:=nil;
    for i:=p3!-triacount!* step -1 until 1 do
    << w := xgetv(p3!-trias!*,i);
       if p3!-tria!-length w > sflimit then
           <<again:=t; p3!-tri!-refine!-one!-tria (w,fcn)>>
        else if not again and p3!-tria!-zerop w
           then z:=car w.z;
    >>;
    if 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.
phase2:
    again:=nil;
    for each w in z do
      if (s:=p3!-tria!-length (w:=xgetv(p3!-trias!*,w))) >limit2 then
       <<for each q in p3!-tri!-refine!-one!-tria (w,fcn) do
           if p3!-tria!-zerop q and not memq(car q,zn)
              then zn:=car q.zn;
         again:=t>>
      else if p3!-tria!-zerop w and not memq(car w,zn)
          then zn:=car w.zn;
    z:=zn; zn:=nil;
    if again 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.
phase3:
    for i:=p3!-triacount!* step -1 until 1 do
      p3!-refine!-phase3(i,i,6,fcn,limit2*0.3);

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

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

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

    for each w1 in z do
       for each w2 in p3!-neighbours w1 do
          if car w2 > car w1 and p3!-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 p3!-plot!-draw(l);
  begin scalar r,s;
  return
   for each w in l collect
   <<r:=nil;
     for each q in w join
       <<s:=p3!-tria!-goodpoint q;
         if r neq s then {r:=s}>>
   >>;
  end;


symbolic procedure p3!-show!-meshes();
 % generate plot of meshes;
  begin scalar d,l,w,s,p1,p2; integer i;
    i:=1;
 loop:
    w:=xgetv(p3!-trias!*,i);
    if null w then 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 plotnum; % Numeric evaluation of algebraic expressions.

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

global '(!*plotinterrupts);

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 not r and not c 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;
      return o;
    >>
   else
    <<
      !*rounded:=cadr d;
      if not caddr d then
         <<!*roundbf:=nil; setdmode('rounded,nil)>>;
      if car(d) then setdmode(car d,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 rdwrap1 car f . rdwrap1 cdr f;

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 {'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 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 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;

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

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

endmodule;

module xvect; % Support for vectors with adaptive length.

% Author: Herbert Melenk, ZIB-Berlin.

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

symbolic procedure xputv(v,n,w);
  begin scalar i,j;
    i:=iquotient(n,1024); j:=iremainder(n,1024);
    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,1024); j:=iremainder(n,1024);
    while(i>0 and v) do
    <<v:=cdr v; i:=i #- 1>>;
    w:=if v then igetv(car v,j);
    return w
  end;

endmodule;

module gnuplot;  % REDUCE interace for gnuplot graphics.

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

global '(!*plotusepipe !*plotpause plotdta!* plotdta!*2 plotmax!*
         plotmin!* plotcmds!* plotcommand!* plotheader!*
         plotcleanup!* plottmp!* lispsystem!* !*plotinterrupts
         plotoptions!* plotdriver!*);

fluid '(plotstyle!*);

plotdriver!*:='gnuplot;

put('gnuplot,'do,'gp!-do);

put('gnuplot,'option,'gp!-option);


!#if (member 'psl lispsystem!*)

   !#if (member 'unix  lispsystem!*,

      in "$reduce/plot/rgpunx.red"$

   !#elif (intersection '(dos os2 winnt alphant) lispsystem!*)

      in "$reduce\plot\rgpwin.red"$

   !#elif (member 'vms lispsystem!*)

      in "$reduce/plot/rgpvms.red"$

   !#endif

!#elif (member 'csl lispsystem!*)

   in "../cslsrc/rgpcsl.red"$

!#endif

endmodule;

%==================================================================

module gnupldrv; % main GNUPLOT driver.

% Author: Herbert Melenk.

global '(!*plotusepipe !*trplot !*plotkeep !*plotrefine
         !*plotinterrupts);

switch  plotusepipe;       % use pipes

switch  trplot;            % list Gnuplot commands to REDUCE
                           % output (e.g. screen).

switch  plotkeep;          % if ON, the command and data files are
                           % not erased after calling Gnuplot.

global '(
        !*plotpause        % Gnuplot pause command at the end:
                           % nil: no pause
                           % -1: Gnuplot will ask the user for
                           %     a return key.
                           % number>0: Gnuplot will wait <number>
                           % seconds.


      plotcommand!*          % string: command to start gnuplot


      plotcmds!*             % file for collecting commands

      plotdta!*            % files for collecting data

      plotheader!*           % list of Gnuplot commands (strings)
                           % for initializing GNUPLOT

      plotcleanup!*          % list of system commands (strings)
                           % for cleaning up after gnuplot

);

if null plotcommand!* then rederr
      " no support of GNUPLOT for this installation";


fluid '(plot!-files!* plotpipe!*);

symbolic procedure gp!-init();
   <<
    plot!-files!* := plotdta!*;
    plotoptions!*:=  nil;
    PlotOpenDisplay();
  >>;

put('gnuplot,'init,'gp!-init);

symbolic procedure plot!-filename();
   <<plot!-files!* := cdr plot!-files!*; u>>
      where u=if null plot!-files!* then
          rederr "ran out of scratch files" else car plot!-files!*;

symbolic procedure gp!-reset();
   if !*plotusepipe and plotpipe!* then
    <<  plotprin2 "exit"; plotterpri();
        close plotpipe!*; plotpipe!*:=nil;>>;

put('gnuplot,'reset,'gp!-reset);

symbolic procedure PlotOpenDisplay();
   begin
    if null plotpipe!* then
    if not !*plotusepipe then plotpipe!* := open(plotcmds!*,'output)
        else <<plotpipe!* :=pipe!-open(plotcommand!*,'output)>>;
    if null plotheader!* then nil else
    if atom plotheader!* then <<plotprin2 plotheader!*; plotterpri()>>
     else if eqcar(plotheader!*,'list) then
      for each x in cdr plotheader!* do <<plotprin2 x; plotterpri()>>
     else typerr(plotheader!*,"gnuplot header");
   end;

symbolic procedure gp!-show();
   if !*plotusepipe and plotpipe!* then
     << channelflush  plotpipe!*; >>
    else
   <<if !*plotpause then plotprin2lt{"pause ",!*plotpause};
     close  plotpipe!*;
     plotpipe!* := nil;
     if plotcommand!* then
       <<plot!-exec plotcommand!*;
         if not !*plotkeep then
            for each u in plotcleanup!* do system u;
       >>;
    >>;

put('gnuplot,'show,'gp!-show);

symbolic procedure plot!-exec u; system u;

symbolic procedure plotprin2 u;
   <<prin2 u; wrs v;
     if !*trplot then prin2 u>> where v=wrs plotpipe!*,!*lower=t;

symbolic procedure plotterpri();
   <<terpri(); wrs v;
     if !*trplot then terpri() >> where v=wrs plotpipe!*;

symbolic procedure plotprin2lt l;
   <<for each x in l do plotprin2 x; plotterpri()>>;

fluid '(plotprinitms!*);

symbolic procedure plotprinexpr u;
   begin scalar plotprinitms!*,!*lower,v;
     !*lower:=t;
     v := wrs plotpipe!*;
     plotprinitms!* := 0;
     if eqcar(u,'file) then
        <<prin2 '!"; prin2 cadr u;prin2 '!"; prin2 " ">>
     else
        errorset(list('plotprinexpr1,mkquote u,nil),nil,nil);
     wrs v;
   end;

symbolic procedure plotprinexpr1(u,oldop);
   begin scalar op;
     if plotprinitms!* > 5 then
        <<prin2 "\"; terpri(); plotprinitms!*:=0>>;
     if atom u then
        <<prin2 if u='e then 2.718281 else
                if u='pi then 3.14159 else u;
          plotprinitms!* := plotprinitms!*+1>>
          else
     if eqcar(u,'!:rd!:) then
         plotprinexpr1 (if atom cdr u then cdr u else
                           bf2flr u,nil)
          else
     if (op:=car u) memq '(plus times difference quotient expt) then
           plotprinexpr2(cdr u,get(car u,'PRTCH),
               oldop and not (op memq(oldop memq
                      '(difference plus times quotient expt)))
               ,op)
          else
     if op='MINUS then
          <<prin2 "(-";
            plotprinexpr1(cadr u,t);
            prin2 ")">>
          else
     if get(car u,'!:RD!:) then
         <<prin2 car u; plotprinexpr2(cdr u,'!, ,t,nil)>>
          else
        typerr(u," expression for printing")
   end;

symbolic procedure plotprinexpr2(u,sep,br,op);
   <<if br then prin2 " (";
     while u do
     <<plotprinexpr1(car u,op);
       u := cdr u;
       if u then prin2 sep>>;
     if br then prin2 ") "
   >>;

symbolic procedure gnuploteval u;
 % Support of explicit calls to GNUPLOT in algebraic mode.
  begin scalar m,evallhseqp!*;
    evallhseqp!* := t;
    m:=plotrounded(nil);
    PlotOpenDisplay();
    for each v in u do
    <<plotprinexpr reval v; plotprin2 " ">>;
    plotterpri();
    plotrounded(m);
  end;

put('gnuplot,'psopfn,'gnuploteval);

% Declare options which are supported by GNUPLOT:

flag ('(

          % keyword options
    contour nocontour logscale nologscale surface nosurface

          % equation type options
    hidden3d xlabel ylabel zlabel title

),'gp!-option);

put('gnuplot,'option,'gp!-option);

symbolic procedure plotpoints u;
  begin scalar f,fn,of,dim,w;
     fn := plot!-filename();
     f := open(fn,'output);
     of := wrs f;
     w:={'plotpoints0,mkquote(nil.u)};
     dim:=errorset(w,t,nil);
     wrs of;
     close f;
     if ploterrorp dim then
        rederr "failure during plotting point set";
     return if car dim=2 then {'file,fn,'x} else {'file,fn,'x,'y};
  end;

symbolic procedure plotpoints0 u;
  begin scalar z,bool;
    integer n;
   for each x in cdr u do
    if not bool and eqcar(x,'list) then n:=plotpoints0 x
      else
     <<bool:=t; n:=n#+1;
       z:=rdwrap reval x;
       if not numberp z then <<wrs nil; typerr(x,"number")>>;
       prin2 z; prin2 " ";
     >>;
   terpri();
   return n;
  end;

symbolic procedure plotpoints1 u;
  begin scalar f,fn,of,y;
     fn :=  plot!-filename();
     f := open(fn,'output);
     of := wrs f;
     for each x in u do
     <<for each y in x do
       << if null y or nil member y then t else
           for each z in y do <<plotprin2number z; prin2 " ">>;
         terpri()
       >>;
       terpri();
     >>;
     wrs of;
     close f;
     return fn;
  end;

symbolic procedure plotprin2number u;
  prin2 if floatp u and abs u < plotmin!* then "0.0" else u;


flag ('(xlabel ylabel zlabel),'plotstring);

symbolic procedure gp!-plotoptions();
  <<if not('polar memq plotoptions!*)then
      plotoptions!* := 'nopolar . plotoptions!*;
    if not('contour memq plotoptions!*)then
      plotoptions!* := 'nocontour . plotoptions!*;
  for each x in plotoptions!* do
    <<plotprin2 "set ";
      if idp x then plotprin2 x else
      <<plotprin2 car x;
        plotprin2 " ";
        if flagp(car x,'plotstring) then plotprin2 '!";
        plotprin2 cdr x;
        if flagp(car x,'plotstring) then plotprin2 '!">>;
      plotterpri()
    >>;
  >>;

symbolic procedure plotstyle1();
   if plotstyle!* then
    <<plotprin2 " \";
     plotterpri();
     plotprin2 " with ";
     plotprin2 plotstyle!*;
     plotprin2 " ";
   >>;

symbolic procedure plottitle option;
  <<plotprin2 "set title ";
    plotprin2 '!";
    plotprin2 option;
    plotprin2 '!";
    plotterpri()>>;

put('title,'gp!-do,'plottitle);

symbolic procedure plotstyle option;
  if option memq'(lines points linespoints impulses dots errorbars)
     then plotstyle!* := option
  else typerr(caddr option, "plot style option");

put('style,'gp!-do,'plotstyle);


% Drivers for different picture types.

symbolic procedure gp!-2exp(x,y,pts,fp);
  % x:   name of x coordinate,
  % y:   name of y coordinate,
  % pts: list of computed point sets,
  % fp:  list of user supplied point sets.
  begin scalar cm;
     plotoptions!* := 'noparametric .  plotoptions!*;
     plotprin2lt{"set size 1,1"};
     plotprin2lt{"set xlabel ",'!",x,'!"};
     plotprin2lt{"set ylabel ",'!",y,'!"};
     gp!-plotoptions();
     plotprin2lt{"set nokey"};
     if pts or fp then plotprin2 "plot ";

     for each f in reversip pts do
     << if cm then <<plotprin2 ",\"; plotterpri()>>;
        plotprin2 "'"; plotprin2 plotpoints1 f; plotprin2 "'";
        plotprin2 " with lines"; cm:=t;
     >>;
     if fp then
     << if cm then <<plotprin2 ",\"; plotterpri()>>;
        plotprin2 "'"; plotprin2 fp; plotprin2 "'";
        plotprin2 if cm then " with points" else " with lines";
     >>;
     plotterpri();
  end;

put('gnuplot,'plot!-2exp,'gp!-2exp);

symbolic procedure gp!-3exp(x,y,z,f);
 % x:   name of x coordinate,
 % y:   name of y coordinate,
 % z:   name of z coordinate,
 % f:   orthogonal list of point lists.
  begin scalar h;
     h:=member('hidden3d,plotoptions!*);
     if h then f:=for each l in f collect
       for each p in l collect {caddr p};
     f:=plotpoints1 f;
     plotprin2lt{"set nohidden3d"};
     if not h then plotoptions!* := 'parametric .
           delete('noparametric,plotoptions!*)
        else
            plotoptions!* := 'noparametric .
           delete('parametric,plotoptions!*);
     plotprin2lt{"set view 60,30,1,1"};
     plotprin2lt{"set size 1,1"};
     if h then plotprin2lt{"set format xy ",'!",'!"};
     plotprin2lt{"set xlabel ",'!",x,'!"};
     plotprin2lt{"set ylabel ",'!",y,'!"};
     plotprin2lt{"set zlabel ",'!",z,'!"};
     gp!-plotoptions();
     plotprin2lt{"set nokey"};
     plotprin2 "splot ";
     plotprin2 "'"; plotprin2 f; plotprin2 "'";
     plotprin2 " with lines ";
     plotterpri();
     plotprin2lt{"set nohidden3d"};
     plotprin2lt{"set format xy"};
  end;

put('gnuplot,'plot!-3exp!-reg,'gp!-3exp);

symbolic procedure gp!-2imp(x,y,l,g,xmin,xmax,ymin,ymax);
 % x,y:   names of coordinates,
 % l:     point lists for funtion,
 % g:     nil or point lists for grid,
 % xmin..ymax: minimum and maximum coordinate values.
  begin scalar f,q;
    q:={{xmin,ymin},nil,{xmin,ymax},nil,
        {xmax,ymin},nil,{xmax,ymax}};
    plotoptions!* := 'noparametric .  plotoptions!*;
    f:=plotpoints1 (q.l);
    plotprin2lt{"set size 1,1"};
    plotprin2lt{"set xlabel ",'!",x,'!"};
    plotprin2lt{"set ylabel ",'!",y,'!"};
    gp!-plotoptions();
    plotprin2lt{"set nokey"};
    plotprin2 "plot "; plotprin2  "'"; plotprin2 f; plotprin2 "'";
    plotprin2 " with lines";
    if g then
    <<plotprin2 ", '"; plotprin2 plotpoints1 g;
      plotprin2 "' with lines";
    >>;
    plotterpri();
  end;

put('gnuplot,'plot!-2imp,'gp!-2imp);

endmodule;

end;