No BSD License  

Highlights from
Advanced Mathematics and Mechanics Applications Using MATLAB, 3rd Edition

image thumbnail

Advanced Mathematics and Mechanics Applications Using MATLAB, 3rd Edition

by

 

14 Oct 2002 (Updated )

Companion Software (amamhlib)

[x,y,z,rmin,dmin,d]=frusdist ...
function [x,y,z,rmin,dmin,d]=frusdist ...
         (rb,rt,h,nb,nt,ns,nc,re,vecs,r0)
% 
% [x,y,z,rmin,dmin,d]=
%      frusdist(rb,rt,h,nb,nt,ns,nc,re,vecs,r0)
% ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
%
% This function determines the point on the 
% surface of a conical frustum which is closest 
% to another point outside the surface. A grid 
% of points on the surface is formed along with 
% a matrix of distances from the remote point to 
% the surface points. The closest point is 
% determined by sorting elements of the 
% distance matrix. The surface and the line from 
% the nearest point to the closest surface point 
% are drawn. A surface is also plotted to depict 
% the distance matrix as a function of the axial 
% and circumferencial variables used to 
% construct the surface. For clarity the 
% negative of the distance matrix is plotted so 
% the peak on the surface identifies the point 
% having minimum distance.
% 
% rb,rt - the base and top radii of the frustum.
%         The base radius can be larger or 
%         smaller than the top radius.
% h     - the height of the frustum
% nb,nt - the number of plot increments taken 
%         on the base and on the top
% ns,nc - the number of increments taken on 
%         the side and around the circumference.
% re    - a vector to the center of the base
% vecs  - a matrix having two columns which 
%         define the orientation of the conical 
%         frustum. vecs(:,1) gives the direction 
%         of the axis of the frustum. The cross 
%         product of vecs(:,1) and vecs(:,2) 
%         gives the direction of the third base 
%         vector defining a triad of local base
%         vectors centered at re on the base. 
%         (See function frustum for more detail.)
% r0    - coordinate vector of the remote point 
%         for which the closest point on the 
%         frustum is sought
%
% x,y,z - matrices of coordinate points on the 
%         surface of the frustum
% rmin  - the vector for the point closest to 
%         the remote point r0 
% dmin  - the shortest distance  
% d     - the matrix containing distances from 
%         the outside point to surface points. 
%
% User m functions called: 
%    frustum, cubrange, surfesy

% Default data case
if nargin==0
  rb=1; rt=2; Rb=2*[1,1,1]; Rt=2*[3,3,3]; 
  va=Rt-Rb; h=norm(va); re=Rb; 
  vecs=[va(:),[1;-1;0]]; r0=[4,0,8];
  nb=10; nt=10; ns=20; nc=35;
end

% Call function frustum to generate points on 
% the surface
r=[rb;rt]; n=[nb;ns;nt;nc]; 
[x,y,z]=frustum(r,h,n,vecs);
x=x+re(1); y=y+re(2); z=z+re(3);

% Form the matrix containing distances from 
% the outside point to surface points. 
d=sqrt((x-r0(1)).^2+(y-r0(2)).^2+(z-r0(3)).^2); 

% Compute the minimum and the related 
% surface point
[dmin,J]=min(min(d)); [dmin,I]=min(min(d'));
rmin=[x(I,J);y(I,J);z(I,J)]; R=[rmin,r0(:)];

% Generate points to plot the line from the 
% outside point to the nearest surface point. 
R=[linspace(rmin(1),r0(1),50);
   linspace(rmin(2),r0(2),50);
   linspace(rmin(3),r0(3),50)];

% Create a window range for undistorted 
% plotting
v=cubrange([[x(:),y(:),z(:)];r0(:)']);

% Draw the line and then the conical frustum
hold off, close, colormap('default');
[M,N]=size(R);
plot3(R(1,:),R(2,:),R(3,:),'-', ...
      R(1,N),R(2,N),R(3,N),'o');
hold on; 
surfesy(x,y,z,'x axis','y axis','z-axis',...
        'Closest Point on a Surface',v);
axis(v); figure(gcf); hold off;
print -deps closept
input('Press [Enter] to continue','s'); 

% Draw a surface showing -d on which the 
% highest point identifies the minimum distance 
% point
[naxial,mcircum]=size(d); N=(1:naxial)/naxial;
M=(1:mcircum)/mcircum; surf(M,N,-d); 
title('Inverted Minimum Distance Surface');
xlabel('circumferential coordinate');
ylabel('radial coordinate');
zlabel('negative of d matrix'); figure(gcf);
print -deps minsurf

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

function [x,y,z]=frustum(r,h,n,vecs)
%
% [x,y,z]=frustum(r,h,n,vecs)
% ~~~~~~~~~~~~~~~~~~~~~~~~~~~
%
% This function creates points defining a 
% frustum of a cone having its axis in 
% the z direction oriented relative to different 
% axes depending on parameters specified in 
% argument vecs.
%
% r    - a vector containing [rb,rt] where rb 
%        is the base radius and rt is the top 
%        radius. If only one number is input 
%        then rt is set equal to rb.
% h    - height of the frustum
% n    - a vector defining the number of radial 
%        increments taken on the base, the 
%        side, the end and the circumference. 
%        n has the form [nb,ns,nt,nc]. 
%        Using [1,1,1,4] generates a cube.
% vecs - a matrix having three rows and either 
%        one or two columns which determines 
%        the axis orientation of the frustum. 
%        When vecs is not input the surface has 
%        its base on the the xy plane and its
%        longitudinal axis along the z 
%        direction. If vecs is present, the 
%        longitudinal axis is in the direction 
%        of vecs(:,1) and the shifted y axis 
%        is in the direction of vecs(:,1) 
%        crossed into vecs(:,2).  
% 
% x,y,z- matrices of coordinate points on the
%        surface of the frustum
%
% User m functions called:  none
%----------------------------------------------

if nargin<4, vecs=[]; end

% The default data creates a cube of unit 
% side length
if nargin==0 
  r=1/sqrt(2)*[1,1]; h=1; n=[1,1,1,4]; 
end 

rb=r(1); nb=n(1); ns=n(2); rt=rb;
if length(r)==1, rt=rb; else rt=r(2); end
if length(n)>2, nt=n(3); else, nt=nb; end
if length(n)>3, nc=n(4); else, nc=36; end

% Generate radius values for rotation about 
% the z axis
R=[linspace(0,rb*(nb-1)/nb,nb), ...
   linspace(rb,rt,ns+1),...
   linspace(rt*(nt-1)/nt,0,nt)]';
z=[zeros(1,nb),linspace(0,h,ns+1),h*ones(1,nt)]';
z=z(:,ones(1,nc+1));

% Make a surface of revolution by rotation 
% about the z axis
th=linspace(pi/nc,pi/nc+2*pi,nc+1);
x=R*cos(th); y=R*sin(th);
if nargin<4 | isempty(vecs), return, end

% If vecs is present shift the axes of the 
% frustum appropriately.
[N,M]=size(x); e3=vecs(:,1); e3=e3/norm(e3); 
if size(vecs,2)==1
  u=null(e3'); u=[u(:,1),cross(e3,u(:,1)),e3];
else
  e2=cross(e3,vecs(:,2)); e2=e2/norm(e2); 
  u=[cross(e2,e3),e2,e3];
end
w=[x(:),y(:),z(:)]*u'; x=reshape(w(:,1),N,M);
y=reshape(w(:,2),N,M); z=reshape(w(:,3),N,M);

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

function surfesy(x,y,z,xlab,ylab,zlab,titl,v)
%
% surfesy(x,y,z,xlab,ylab,zlab,titl,v)
% ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
%
% This function provides an easy input 
% interface to function surf.
%
% x,y,z          - data for surface plotting
% xlab,ylab,zlab - labels for the coordinate 
%                  axes
% titl           - a title for the plot
% v              - a vector to set the axis 
%                  range. If no value is input, 
%                  a range is found to make an 
%                  undistorted plot. If a 
%                  single number is input, the 
%                  default scaling is used 
%
% User m functions called:  cubrange
%----------------------------------------------

if nargin<8, v=cubrange(x,y,z); end 
if nargin<7, titl=[]; end
if nargin<6, xlab=''; ylab=''; zlab=''; end
surf(x,y,z); xlabel(xlab); ylabel(ylab);
zlabel(zlab); title(titl);
if length(v)>2, axis(v); end
grid on; figure(gcf);

Contact us