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 . fp else % WN 25.9.98
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;
%% extremely careful here !! see : plot(e^(x/(2-x)),x);
%% WN 29. 9.99
w := errorset({'times2,dy32,dy32},nil,nil);
if ploterrorp w then goto disc else w:=car w;
d := errorset({'times2,dy12,dy12},nil,nil);
if ploterrorp d then goto disc else d:=car d;
w := (w + dx32**2);
d := (d + dx12**2);
d:= errorset({'times2,w,d},nil,nil);
if ploterrorp d then goto disc else d:=car d;
d := sqrt d;
%% original 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;
g := rat := lev := ml := nil;
%%%%% 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
end;