Calculating ellipse arc length via elliptic integral

arc lengthelliptic integrals

I'm looking for an exact equation to calculate the arc length of an ellipse. I've found https://keisan.casio.com/exec/system/1343722259 to correctly calculate the length, but cannot calculate the same result when implementing in Julia via Elliptic.jl or ArbNumerics.

The MathJax equation for the arc length on the Keisan calculator is the same as I've found elsewhere; in Julia:

import Elliptic
function ellipticArcLength(a,b, th0, th1)
  # b<a
  # 0 <= th <= pi/2

  k2 = 1-b^2/a^2
  k = sqrt(k2)
  r2 = th -> a^2*b^2/( b^2*cos(th)^2 + a^2*sin(th)^2 ) #r(th)^2
  r = th -> sqrt(r2(th))
  x = th -> r(th) * cos(th)
  L = a*Elliptic.E( x(th0)/a, k) - a*Elliptic.E( x(th1)/a, k)
  return L
end
 
println("36.874554322338 ?= $(ellipticArcLength(  50, 20, deg2rad(10), deg2rad(60) ))")
println("99.072689284541 ?= $(ellipticArcLength( 500, 20, deg2rad(10), deg2rad(60) ))")

yielding

36.874554322338 ?= 29.051564462163682

99.072689284541 ?= 98.17285637594554

That is, the 50×20 ellipse is significantly off while the 500×20 is close, but I expect the elliptic integral method to be exact.

Can anyone point out my error(s)?

(DLMF provides a 2-part formula which hasn't been close in my attempts.)

Thank you

Best Answer

using ArbNumerics

function ellipticArcLength(a, b, angle )
  phi = atan( a/b*tan(angle))
  m = 1 - (a/b)^2
  return b*elliptic_e( ArbReal(phi), ArbReal(m) )
end
function ellipticArcLength(a, b, start, stop)
  lStart = ellipticArcLength(a,b, start)
  lStop = ellipticArcLength(a,b, stop)
  return lStop - lStart
end

println("36.874554322338 ?= $(ellipticArcLength(  50, 20, deg2rad(10), deg2rad(60) ))")
println("99.072689284541 ?= $(ellipticArcLength( 500, 20, deg2rad(10), deg2rad(60) ))")

Multiplying by the major axis instead of the minor seems to be the root error, see https://math.stackexchange.com/a/1123737/974011 .

Related Question