function rub(n,d,clipvalue,islog,numpoints)
% function rub(n,d,clipvalue,islog,numpoints)
%
% This function draws a 3d rubber sheet surface
%    of the magnitude of a transfer function with
%    numerator n and denominator d.
% clipvalue is how high to limit the infinit peaks
% islog is a boolian 1=log of height
%                    0= linear
% numpoints is how many points you want on
%    the x grid and on the y grid.
%
% This function uses the surf comand and as such you may
%   use the VIEW comand to change view angle.
%
% try n=[1 4]
%      d=[1 12 49]
% rub(n,d,1000,1,40)
% 
% or
% n=poly([-1+2.001j  -1-2.001j]);
% d=poly([-1 -5 -14]);
% rub(n,d,1000,1,40)
%
%   See also MESH, VIEW, FIGURE, AXES.
%
%   Copyright (c) 1998 by Ulta Tech Instuments.
%   $Revision: 2.01 $  $Date: 2003/02/10 $ for Octave
%   $Revision: 1.01 $  $Date: 2002/02/04 $
%   $Revision: 1.0 $  $Date: 1998/01/24 $

oneplot()
gset term windows


np=numpoints-1;
z(numpoints,numpoints)=0;
rn=roots(n);
if isempty(rn);
   rn=0;
end

rd=roots(d);
m(1)=max(abs(real(rn)));
m1(1)=max(abs(imag(rn)));
m(2)=max(abs(real(rd)));
m1(2)=max(abs(imag(rd)));
m1rn=max(m);
m3rn=max(m1);
if m3rn==0
   m3rn=m1rn/4;
end

for k1=0:np
   for k=0:np
      s=-((4+m1rn)*k1/np)+(m3rn*2*(-1+k/(np/2)))*j;
      nv=polyval(n,s);
      dv=polyval(d,s);
      if islog==1
         v=log(abs(nv/dv));
      else
         v=abs(nv/dv);
      end
      
      if v>clipvalue
         v=clipvalue;
      end   
      if v<-clipvalue
         v=-clipvalue;
      end
      
      z(k1+1,k+1)=v;
   end
end
w1=abs(imag(s));
sd=real(s);
w=-(w1):w1/(np/2):w1;
ss=0:sd/np:sd;
mesh(w,-ss,z)
view(-10,20)