Artifact 3b1ccd40424e88cd499942da2fc93cfb920b8ea2e03d302e054a028297eaadc2:
- Executable file
r37/packages/plot/plotexp2.red
— part of check-in
[f2fda60abd]
at
2011-09-02 18:13:33
on branch master
— Some historical releases purely for archival purposes
git-svn-id: https://svn.code.sf.net/p/reduce-algebra/code/trunk/historical@1375 2bfe0521-f11c-4a00-b80e-6202646ff360 (user: arthurcnorman@users.sourceforge.net, size: 8730) [annotate] [blame] [check-ins using] [more...]
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; 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;