ellipse points regularly spaced, or inverse of the incomplete elliptic integral of the second (not first) kind
Felipe G. Nievinski
second (not first) kind? I'm familiar with the Jacobi amplitude as the
inverse of the incomplete elliptic integral of the first kind, but
can't seem to find in the literature an inverse for the second kind.
Alternatively, it'd help to know how to obtain the first kind given
the second one (both incomplete).
I found this statement at <http://dlmf.nist.gov/19.25#v>:
"[elliptic] integrals of the second kind, ... are not invertible
in terms of single-valued functions"
though I find the treatment given, in terms of the so-called elliptic
symmetric integrals, a little hard to digest.
I'm trying to discretize the perimeter of an ellipse at regularly
spaced points.Trivial discretizations end up causing an undesirable
high/low density of points at maximum/minimum distance from the
center. (A brute force approach would be to discretize at a
sufficiently high density everywhere then discard excessive points
where necessary.)
I am given the ellipse semi-major and semi-minor axes, a and b,
respectively; eccentricity is obtained as e=sqrt(1-(b/a)^2). To define
regularly spaced points, I compute the ellipse perimeter or
circumference length C and, given a step length s, define intermediary
arc lengths as len = 0:s:C. Now the problem is to obtain the x,y
coordinates of the ellipse at each such arc length values. (The
circumference is given by C = 4 * a * E(k), where E(k) is the complete
elliptical integral of the second kind, and k=e^2).
A brief review of the background theory. The ellipse arc length can be
calculated as:
len = a * E(phi,k),
in terms of the incomplete elliptical integral of the _second_ kind
v = E(phi,k).
I'm not aware of an explicit formulation for its inverse,
phi = Einv(v,k).
In contrast, for the incomplete elliptical integral of the _first_
kind,
u = F(phi,k)
its inverse is the so-called Jacobi's amplitude function:
phi = am(u,k),
Sometimes the arc length is expressed as:
len = a * Eps(u,k),
in terms of the so-called Jacobi's epsilon function:
Eps(u,k) = E(am(u,k),k).
I don't see the benefit of this formulation, unless we can relate the
output of the first and second elliptic integrals, u and v. And if we
do, Jacobi's amplitude function usually involves the so-called Jacobi
elliptic functions:
sn(u) = sin(phi),
cn(u) = cos(phi),
which can be related directly to the point position by
x = a * sn(u,k),
y = b * cn(u,k).
Putting it all together, we could proceed in one of the following
ways:
v = len / a
invert for phi = Einv(v,k) (how?)
x = a * cos(phi)
y = b * sin(phi)
or
v = len / a
relate u to v (how?)
get sn(u,k) and cn(u,k) given u via routine ellipj.m
x = a * sn(u,k)
y = b * cn(u,k)
I'm inclined to obtain phi = Einv(v,k) via a semi-brute force
approach, i.e., sampling densely v = E(phi,k) then keeping only the
phi values yielding more or less regularly spaced v, e.g.:
% (using <http://elliptic.googlecode.com/svn/trunk/elliptic12.m>)
s = 0.1;
a = 1;
b = 1/2;
esq = 1 - (b/a)^2;
phi = linspace(0, 2*pi)';
[F, E] = elliptic12 (phi, esq);
len = a * E;
dlen = [0; diff(len)];
figure, plotyy(phi,len, phi,dlen)
[ignore, temp] = ellipke (esq);
C = 4 * a * temp;
C - len(phi==2*pi) % DEBUG
len2 = (0:s:C)'; len2(end+1) = C;
phi2 = interp1(len, phi, len2, 'nearest');
dlen2 = [0; diff(len2)];
figure, plotyy(phi2,len2, phi2,dlen2)
figure, hold on, plot(phi,len,'.-k'), plot(phi2,len2,'ok')
figure, hold on, plot(phi,dlen,'.-k'), plot(phi2,dlen2,'ok')
%x = a * cos(phi);
%y = b * sin(phi);
varphi = atan2(a*sin(phi), b*cos(phi));
x = a * cos(varphi);
y = b * sin(varphi);
ind = ismember(phi, phi2);
figure, hold on, plot(x, y, '.-k'), plot(x(ind), y(ind), 'ok'), axis
equal
figure, hold on, plot(x(ind), y(ind), 'o-k'), axis equal
Though if I'm going this way I might just as well avoid integrals
altogether and do simply a pure brute force approach:
s = 0.1;
a = 1;
b = 1/2;
theta = linspace(0, 2*pi)';
x = a * cos(theta);
y = b * sin(theta);
figure, plot(x, y, '.-k'), axis equal
dx = [0; diff(x)];
dy = [0; diff(y)];
dlen = sqrt(dx.^2 + dy.^2);
len = cumsum(dlen);
figure, plotyy(theta,len, theta,dlen), title('len')
C = len(phi==2*pi);
len2 = (0:s:C)'; len2(end+1) = C;
theta2 = interp1(len, theta, len2, 'nearest');
dlen2 = [0; diff(len2)];
figure, plotyy(theta2,len2, theta2,dlen2), title('E')
figure, hold on, plot(theta,dlen,'.-k'), plot(theta2,dlen2,'ok')
ind = ismember(theta, theta2);
figure, hold on, plot(x, y, '.-k'), plot(x(ind), y(ind), 'ok'), axis
equal
figure, hold on, plot(x(ind), y(ind), 'o-k'), axis equal
References:
<http://mathworld.wolfram.com/EllipticIntegraloftheFirstKind.html>
<http://mathworld.wolfram.com/EllipticIntegraloftheSecondKind.html>
<http://mathworld.wolfram.com/JacobiAmplitude.html>
<http://mathworld.wolfram.com/JacobiEllipticFunctions.html>
<http://mathworld.wolfram.com/Ellipse.html>
<http://en.wikipedia.org/wiki/Ellipse>
<http://en.wikipedia.org/wiki/Jacobi's_elliptic_functions>
<http://en.wikipedia.org/wiki/Elliptic_integral>
<http://reference.wolfram.com/mathematica/tutorial/
EllipticIntegralsAndEllipticFunctions.html>
<http://www.mhtl.uwaterloo.ca/courses/me755/web_chap3.pdf>
<http://dlmf.nist.gov/19>
<http://dlmf.nist.gov/22>
moiseev.igor
http://www.mathworks.de/matlabcentral/newsreader/view_thread/283132
Found some discussion on inverting elliptic integrals:
http://mathforum.org/kb/message.jspa?messageID=1672065&tstart=0
and final expansion procedure
http://mathforum.org/kb/message.jspa?messageID=1672067&tstart=0
From you refs (http://dlmf.nist.gov/19.25#v) got that it seems
possible to invert the incomplete elliptic integral of the second kind
using the Carlson symmetric integrals and related functions, but the
numeric algorithm to perform it fast as is done by Gauss-Legendre's
algorithm for the first kind integrals unfortunately does not exist.
Hope it helps.