Ellipse arc length

Computes the length of an ellipse arc given the two radiuses and the starting and ending angles.

Description

I had to compute the length of arcs of Earth's meridian. I noticed that there was no analytical solution for this so I wrote this very simple function.

It computes the arc length of an ellipse centered on (0,0) with radius a (along OX) and radius b (along OY)

x(t) = a.cos(t)
y(t) = b.sin(t)

with angle t (in radians) between t1 and t2.

The solution is obtained numerically by dividing the arc in small straight segments. In addition, the second output is the length of the complete ellipse compyted with an approximated formula (Ramanujan).

Luc

• MATLAB release: MATLAB 7.5 (R2007b)

[Dana Smith]

It would seem to me that before you can consider a question of the arc length of an ellipse, you would first ask, what is the finite number of arcs whose assembly comprises an ellipse? In other words, dividing an ellipse along the major and/or axes, the resultant curved line is composed of how many arcs and how do you determine the end points of each?

[Cameron]

Try:

dt = (t2-t1)/nseg;
tt = t1+.5*dt:dt:t2;
dx = -a*dt*sin(tt);
dy = b*dt*cos(tt);
L = sum(sqrt(dx.^2+dy.^2));

It's faster this way, but the result is slightly different (proportional and more accurate). I'm not sure why.

[John D'Errico]

This code rather blithely generates 100,000 linear segments (the same number for all problems) then uses a piecewise linear approximation. Far more logical is to just call a tool designed to solve numerical integration problems.

I'll use the case provided by the author as an example. Here they are computing an approximate perimeter around a half-great circle of the earth that passes through the poles. It looks like the units are kilometers.

>> Req=6378.137;
>> Rpo=6356.752;
>> ellipsearc(Req,Rpo,-pi/2,pi/2)
ans =
          20003.9309621345

This seems reasonable. How correct is it? As I said, we can use a tool better designed to solve the problem than piecewise linear arcs. Here I'll do the entire computation in two lines. Define the perimeter of the ellipse parametrically in terms of angular displacement, and form differentials. Then just call quadgk.

>> arcfun = @(theta,a,b) sqrt((a.*sin(theta)).^2 + (b.*cos(theta)).^2);
>> quadgk(@(t) arcfun(t,Req,Rpo),-pi/2,pi/2)
ans =
          20003.9309654251

As we see, the two results compare reasonably well, but quadgk is returning a result that differs in the last 5 decimal places. Is the quadgk call more accurate, or less so? This is easily seen by computing the arc length of a semi-circular arc. We know that to be simply pi.

>> quadgk(@(t) arcfun(t,1,1),-pi/2,pi/2)
ans =
          3.14159265358979

>> ellipsearc(1,1,-pi/2,pi/2)
ans =
          3.14159265307115

And for those who cannot recall the digits of pi...

>> pi
ans =
          3.14159265358979

So it appears that ellipsearc is only returning about 10 correct digits, whereas quadgk got all of them correct. If you are going to make the effort to generate 100k points on that ellipse, the result should be nearly perfect! To me, this is a major flaw. Use a better algorithm.

Which approach is faster? I'll concede that ellisearc is reasonably coded. It is vectorized internally, in the sense that no explicit loops are used. (IT is NOT externally vectorized, so that if you had multiple arcs to integrate, the code will bomb if you tried to pass in a vector of parameters.) I'll compare the times using Steve Eddins timeit tool, which will give an accurate prediction.

>> timeit(@() quadgk(@(t) arcfun(t,2,3),-pi/2,pi/2))
ans =
            0.000734411552

>> timeit(@() ellipsearc(2,3,-pi/2,pi/2))
ans =
            0.001485160487

So quadgk is twice as fast, at least on this test problem. So if you really want to compute the arc length of an elliptical arc, just use quadgk as I've shown.

How about the code itself? Well, it does have some error checks. It substitutes defaults for two of the parameters if they are not supplied.

Mlint flags a couple of things. For example, there is no need to have a return statement after an error is generated, since the error automatically drops you out. As well, the code uses & in an if clause, instead of &&, which is more efficient. In neither case are these major problems, but both cases are a bit sloppy. Good code has these sloppy things cleaned up.

A more sloppy point that I don't like at all is the internal comment where 100,000 points are generated on the ellipse. Here it is:

%number of straight segments (complete ellipse)
nseg=100000; % you may adjust this parameter

Good code does not have you edit the code to change a potentially important parameter.

How about the help? I'd give the help a fair grade. It is certainly not good, but also not terribly bad. There was no explicit statement of the parameters for the code, what they are, what they mean. In fact, it is necessary to edit the code to learn what the parameters are and their order. This part I very much dislike. help (or doc) should tell you everything you need to know.

I did find several examples of use in the help. This is good. There was even a wiki-link about the ellipse.

I did find an H1 line, so when I try out lookfor, I get this response:

>> lookfor ellipse
ellipsearc - computes the arc length of an ellipse of equation:

Overall, I want to give the code a fair rating, somewhere in the middle of the scale. It is better than a 1 star code, and certainly not worth 5 stars. Even 4 stars would be too much credit for code that makes you edit the code to learn to use it. So it must be either 2 or 3 stars. There are a couple of mlint flags, but only minor things. You should not need to edit code to change a parameter. And the use of what I'd consider a poor numerical scheme for the integration is another demerit in my eyes. All of this is pushing me to rate it as 2 stars. On the other hand, the author did make an effort in the help, did have an H1 line, did have internal comments, and did use a vectorized scheme in the code. It is even reasonably fast, moderately competitive in speed with my own solution, though 5 digits less accurate. As much as I want to rate it 2.5 stars or so, it does not quite rise up to 3 in my eyes. I'd be willing to revisit the code when changes are made though, and new ratings by a reviewer supersede old ones.