[xi,eta,x,y]=evengrid(a,b,nxi,nu)

function [xi,eta,x,y]=evengrid(a,b,nxi,nu)
% [xi,eta,x,y]=evengrid(a,b,nxi,nu)
% This function generates an approximately 
% uniform grid of points on an ellipse
% a,b    - ellipse semi-diameters
% nxi    - number of evenly spaced points
%          generated on the y-xias
% nu     - number of boundary points used
%          to compute arc length (default
%          value is 1001)
% xi,eta - curvilinear coordinate vectors
% x,y    - Cartesian coordinate arrays

%          HBW, 12/01/04

if nargin<4, nu=1001; end
if nargin==0, a=1; b=.5; nxi=8; end

% Compute xi to make the points on the y-axis
% evenly spaced
h=sqrt(a^2-b^2); xi=asinh(linspace(0,b/h,nxi));

% Compute the arc length
etau=linspace(0,2*pi,nu)';
z=a*cos(etau)+i*b*sin(etau);
s=[0;cumsum(abs(diff(z)))];

%Make point spacing even on the boundary
neta=ceil(s(end)/b*nxi);
if mod(neta,2)==0, neta=neta+1; end
eta=interp1(s,etau,linspace(0,s(end),neta)');

% Generate the point grid
zeta=xi(ones(neta,1),:)+i*eta(:,ones(1,nxi));
z=h*cosh(zeta); x=real(z); y=imag(z);