meshanellip.m

function [node,face,elem]=meshanellip(c0,rr,tsize,maxvol)
%
% [node,face,elem]=meshanellip(c0,rr,opt)
%
% create the surface and tetrahedral mesh of an ellipsoid
%
% author: Qianqian Fang, <q.fang at neu.edu>
%
% input: 
%   c0:  center coordinates (x0,y0,z0) of the ellipsoid
%   rr:  radii of an ellipsoid, 
%        if rr is a scalar, this is a sphere with radius rr
%        if rr is a 1x3 or 3x1 vector, it specifies the ellipsoid radii [a,b,c]
%        if rr is a 1x5 or 5x1 vector, it specifies [a,b,c,theta,phi]
%           where theta and phi are the rotation angles along z and x 
%           axes, respectively. Rotation is applied before translation.
%   tsize: maximum surface triangle size on the sphere
%   maxvol: maximu volume of the tetrahedral elements
%
% output:
%   node: node coordinates, 3 columns for x, y and z respectively
%   face: integer array with dimensions of NB x 3, each row represents
%         a surface mesh face element 
%   elem: integer array with dimensions of NE x 4, each row represents
%         a tetrahedron; if ignored, only produces the surface
%
% example:
%   [node,face,elem]=meshanellip([10,10,-10],[30,20,10,pi/4,pi/4],0.5,0.4);
%   plotmesh(node,elem,'x>10');axis equal;
%
% -- this function is part of iso2mesh toolbox (http://iso2mesh.sf.net)
% 

if(nargin<3)
    error('you must at least provide c0, rr and tsize, see help for details');
end

rr=rr(:)';
if(length(rr)==1)
    rr=[rr,rr,rr];
elseif(length(rr)==3)
    % do nothing
elseif(length(rr)~=5)
    error('the number of elements for rr is not correct. see help for details');
end

rmax=min(rr(1:3));
if(nargout==3)
    if(nargin==3)
	maxvol=tsize*tsize*tsize;
    end
    [node,face,elem]=meshunitsphere(tsize/rmax,maxvol/(rmax*rmax*rmax));
else
    [node,face]=meshunitsphere(tsize/rmax);
end

node=node*diag(rr(1:3));
if(length(rr)==5)
   theta=rr(4);
   phi=rr(5);
   Rz=[cos(theta) sin(theta) 0;-sin(theta) cos(theta) 0;0 0 1];
   Rx=[1 0 0;0 cos(phi) sin(phi);0 -sin(phi) cos(phi)];
   node=(Rz*Rx*(node'))';
end
node=node+repmat(c0(:)',size(node,1),1);