orthdisk

 node=orthdisk(c0,c1,r,ndiv)

 Defining a 3D disk that is orthogonal to the vector c1-c0 

 author: Qianqian Fang (fangq <at> nmr.mgh.harvard.edu)

 input:
     c0: a 1x3 vector for the origin
     c1: a 1x3 vector to define a direction vector c1-c0
     r: the radius of the disk that is orthogonal to c1-c0, passing through c0
     ndiv: division count to approximate a circle by a polygon, if ignored, ndiv=20
     v1: a 1x3 vector specifying the first point on the output 3D disk. if
         v1 is not perpendicular to c1-c0, the disk rotation axis is
         changed to cross(v1,cross(c1-c0,v1)); 
     angle0: angle0 represents the angle (in radian) of the 1st point in 
         the 3D disk if placed on the x-y plane.

 output:
     node: the 3D vertices of the disk

 -- this function is part of iso2mesh toolbox (http://iso2mesh.sf.net)

0001 function node=orthdisk(c0,c1,r,ndiv,v1,angle0)
0002 %
0003 % node=orthdisk(c0,c1,r,ndiv)
0004 %
0005 % Defining a 3D disk that is orthogonal to the vector c1-c0
0006 %
0007 % author: Qianqian Fang (fangq <at> nmr.mgh.harvard.edu)
0008 %
0009 % input:
0010 %     c0: a 1x3 vector for the origin
0011 %     c1: a 1x3 vector to define a direction vector c1-c0
0012 %     r: the radius of the disk that is orthogonal to c1-c0, passing through c0
0013 %     ndiv: division count to approximate a circle by a polygon, if ignored, ndiv=20
0014 %     v1: a 1x3 vector specifying the first point on the output 3D disk. if
0015 %         v1 is not perpendicular to c1-c0, the disk rotation axis is
0016 %         changed to cross(v1,cross(c1-c0,v1));
0017 %     angle0: angle0 represents the angle (in radian) of the 1st point in
0018 %         the 3D disk if placed on the x-y plane.
0019 %
0020 % output:
0021 %     node: the 3D vertices of the disk
0022 %
0023 % -- this function is part of iso2mesh toolbox (http://iso2mesh.sf.net)
0024 %
0025 
0026 len=sqrt(sum((c0-c1).*(c0-c1)));
0027 v0=c1-c0;
0028 
0029 if(nargin>=5)
0030     vt=cross(v0,v1);
0031     if(abs(dot(v0(:),v1(:)))>1e-5) % input is not orthogonal
0032         v0=cross(v1,vt);
0033     end
0034 end
0035 
0036 if(nargin<6)
0037     angle0=0;
0038 end
0039 
0040 if(nargin<=2)
0041     r=1;
0042 end
0043 if(nargin<=3)
0044     ndiv=20;
0045 end
0046 
0047 dt=2*pi/ndiv;
0048 theta=angle0+dt:dt:2*pi+angle0;
0049 cx=r*cos(theta);
0050 cy=r*sin(theta);
0051 pp=[cx(:) cy(:) zeros(ndiv,1)];
0052 node=rotatevec3d(pp,v0)+repmat(c0(:)',size(pp,1),1);