Elliptical math problems

[Ariphaos]

Okay, I've spent quite a bit of the day on this and I'm stuck. I'm trying to find two things:

1: Given the circumference P and major axis a of a symmetrical (foci equidistant from the origin) ellipse, find the minor axis.
- Approximations are fine, I'm just trying to find something better than
b = P/π - a
- I tried inverting Ramanujan's equation, but I got stuck at
b(6c - 8a - 6b) = 6aa - 6ac - c

Where c = P/π. I'm not even sure if it's worth the effort if there's a better equation to start from in the first place.

2: Given an arc length on the above ellipse no more than a quarter of its circumference, find the angle subtended by the length.
* Wikipedia would suggest it's one of these
* This method would be interesting but I can't figure a way to know the specific r at a given arc-length point around the ellipse without already knowing the angle.
* This discussion seems to apply but I'm not certain that's what I want?

[Howedar]

This sounds like a job for a new coordinate system, but that's just a random near-bed thought.

I was afraid you wanted help on elliptic PDE's...

[Kuroneko]

1: Given the circumference P and major axis a of a symmetrical (foci equidistant from the origin) ellipse, find the minor axis.
- Approximations are fine, I'm just trying to find something better than
b = P/π - a

Alright, how about P ≅ cπ√2, were c² = a²+b²?

- I tried inverting Ramanujan's equation, ...

When one is talking about someone as prolific as Ramanajuan, please be more specific--Ramanujan had at least two approximations for this. Based on your comments, I you're probably talking about P ≅ π[3(a+b) - sqrt{(a+3b)(3a+b)}], where a,b are semimajor and semiminor axes, respectively, this has solution Δ = 3P²+12πaP-20π²a², b = [3P-4πa+sqrt(Δ)]/(6π). On the other hand, if you intend to use his other approximation, good luck with that.

2: Given an arc length on the above ellipse no more than a quarter of its circumference, find the angle subtended by the length.
* Wikipedia would suggest it's [an elliptic integral]

It is, in fact: the arc length of an arc θ from the minor axis is given by the elliptic integral of the second kind, E(φ,e), where e is the eccentricity of the ellipse and tan φ = (a/b)tan θ. In general, you must make some sort of restriction such as the starting point of the arc; otherwise, the problem is ill-defined.

[Edit×2: corrected φ + factor of π]

* This method would be interesting but I can't figure a way to know the specific r at a given arc-length point around the ellipse without already knowing the angle.

That looks like too much work.

* This discussion seems to apply but I'm not certain that's what I want?

Mr. Mathar's series solution would be applicable, assuming it is correct, but please make note that he defined the elliptic integral of the second kind slightly differently than usual.

[Kuroneko]

This sounds like a job for a new coordinate system, but that's just a random near-bed thought.

It's a good idea for seeing just why the elliptic integral of the second kind is important in the first place, but it doesn't really help in getting a closed-form inversion.

If we take the ellipse with eccentricity ε<1 as having the foci along x = ±a and make the transformation z = x+iy = aε cos iw, where w = u+iv. Thus, x = aε cosh u sin u, y = aε sinh u sin v, so that the points on the ellipse share the same u-coordinate acosh(1/ε). The metric is ds² = dx² + dy² = a²dv²[1 - ε²cos²v] at constant u = acosh(1/ε). Integrating ds, we have the elliptic integral of the second kind, with sin replaced with cos (which means I need to correct my previous post--getting sin back would mean measuring angles from the minor axis). Now, we have tan θ = y/x = tanh u tan v, so at u = acosh(1/ε), tan θ = sqrt[1-ε²] tan v = [b/a] tan v.

[Ariphaos]

Alright, how about P ≅ cπ√2, were c² = a²+b²?

Well, I was gunning for a somewhat higher standard for 'better'. Namely, the goal was something that's not going to get noticed under double precision - Ramanajan qualifies.

When one is talking about someone as prolific as Ramanajuan, please be more specific--Ramanujan had at least two approximations for this. Based on your comments, I you're probably talking about P ≅ π[3(a+b) - sqrt{(a+3b)(3a+b)}], where a,b are semimajor and semiminor axes, respectively, this has solution Δ = 3P²+12πaP-20π²a², b = [3P-4πa+sqrt(Δ)]/6. On the other hand, if you intend to use his other approximation, good luck with that.

Wow, yes, and sorry. Thank you very much - I'll have to try working that out, I feel rusty now.

It is, in fact: the arc length of an arc θ from the minor axis is given by the elliptic integral of the second kind, E(φ,e), where e is the eccentricity of the ellipse and tan φ = (b/a)tan θ. In general, you must make some sort of restriction such as the starting point of the arc; otherwise, the problem is ill-defined.

Err, isn't that finding the length then, and not the angle?

Say I begin with an ellipse centered on the origin with semimajor axis a and semiminor axis b, with the foci located at +/- c along the x-axis (probably not important). I have an arc of known length starting at either [0,b] (or [a,0], as long as it's consistent) - but I want to find the angle. No sort of θ is known, unless I'm being wrongheaded and am supposed to build a unit-ellipse of some sort? Since I know the circumference I can make that equal 2π.

Ultimately, I can also work with calculating the arc length given theta (from the same points), but that's far from the optimal solution.

That looks like too much work.

That too.

Mr. Mathar's series solution would be applicable, assuming it is correct, but please make note that he defined the elliptic integral of the second kind slightly differently than usual.

Slightly differently?

[Ariphaos]

solution Δ = 3P²+12πaP-20π²a², b = [3P-4πa+sqrt(Δ)]/6.

This appears to be off by a factor of π? It seems it should be
b = [3P-4πa+sqrt(Δ)]/(6π)

[Kuroneko]

This appears to be off by a factor of π? It seems it should be b = [3P-4πa+sqrt(Δ)]/(6π)

It's a typo; yes, b = [3P-4πa+sqrt(Δ)]/(6π).

Err, isn't that finding the length then, and not the angle?

Yes--the problem is then inverting the elliptic integral of the second kind. That's the series Mr. Mathar is trying to build in that post.

Say I begin with an ellipse centered on the origin with semimajor axis a and semiminor axis b, with the foci located at +/- c along the x-axis (probably not important). I have an arc of known length starting at either [0,b] (or [a,0], as long as it's consistent) - but I want to find the angle.

Right. Since tan θ = (b/a)tan v, this problem is then finding v such that the given arc length is s = E(v,ε) for constant eccentricity ε = sqrt[1-b²/a²], θ measured from (0,b). Sorry for the confusion about θ and v/φ; I shouldn't have relied on memory when posting the first time. My second post contains the correct derivation.

Ultimately, I can also work with calculating the arc length given theta (from the same points), but that's far from the optimal solution.

If you're willing to do that, Newton's method would be a good idea, since dE/dv = sqrt[1-ε²sin²v] is comparatively simple.

Slightly differently?

His is E(φ,m) = Int{0,φ}[ sqrt{1-m sin²α}dα ], so that m is the square of the elliptic modulus, i.e., the square of the eccentricity rather than the eccentricity itself.

[Ariphaos]

Urgh.

Taking this a step at a time (namely so I can check my work before I try the inversion), I worked with the expansion from
here.

Which declares an infinite series of integrals of:

(-1)^n C(½,n) e^2n sin (2na)

For n = 0, 1, 2 ...

Which is very nice but since I have no idea what a combination of ½ is supposed to be, I'm plugging in numbers half-blind. The following loop still has a ~.1% error or so:

Code: Select all

long double arclength (long double ldTheta, long double ldEccentricity, unsigned short iterations)
  {
    long double t = 0, s = 0;
    unsigned short n = 0, k = 0;

    long double ldResult = 0;

    for (n = 0; n <= iterations; n++)
      {
        t = (long double) combination (2*n, n);
        for (s = 0, k = 1; k <= n; k++)
          {
            s += pow (-1.0, k) * combination (2*n, n - k) * sin (2*k * ldTheta) / k;
          }
        t = ((ldTheta * t * t / pow (2.0, 2*n) + s) / pow (2.0, 2*n)) * pow (ldEccentricity, 2*n) / (1 - 2*n);
        ldResult += t;
      }

    return ldResult;
  }

It functions perfectly with Theta at pi/2 or 0, of course, but I'm obviously not treating the sine summation correctly. :-/

[Kuroneko]

Which declares an infinite series of integrals of:
(-1)^n C(½,n) e^2n sin (2na) For n = 0, 1, 2 ...
Which is very nice but since I have no idea what a combination of ½ is supposed to be, I'm plugging in numbers half-blind.

For non-negative integer n, C(1/2,n) = Γ(3/2)/[Γ(3/2-n)Γ(n+1)] = (-1)^{n+1} (2n-3)!!/[2^n n!] = (-1)^{n+1} (2n-2)!/[2^{2n-1} n! (n-1)!]. You can implement the general case of this rather simply:

Code: Select all

typedef long double     ldouble;
ldouble _ldChoose(ldouble x, unsigned int n) {
  ldouble res = 1.0; unsigned int k = 0;
  while ( k < n ) { res *= (x-(ldouble)(k))/((ldouble)(k+1)); k++; }
  return res;
}

Edit: Oh yes, you can see from the last form of C(1/2,n) that it is also equal to (-1)^{n+1} C(2n,n) / [ 4^n (2n-1) ]. That's probably the easiest one to adapt without making a choose function for real arguments.

[Ariphaos]

...choose function for real arguments.

Is there a good explanation for that? What that means, I mean...

Anyway, thanks, that makes my error obvious:

Code: Select all

long double arclength (long double ldTheta, long double ldEccentricity)
  {
    long double t = 0, sines = 0, c = 0, pow4n = 0;
    long double parameter = ldEccentricity * ldEccentricity;
    unsigned short n = 0, n2 = 0, k = 0;
    short s = 0;

    long double ldResult = 0;

    for (n = 0; n < 10; n++)
      {
        n2 = n << 1;
        pow4n = pow (4.0, n);
        for (sines = 0, k = 1; k <= n; k++)
          {
            s = (k % 2)?-1:1;
            sines += ((long double) s) * combination (n2, n - k) * sin ((k << 1) * ldTheta) / ((long double) k);
          }
        c = -1.0 * combination (n2, n) / (pow4n * (n2-1));
        t = ldTheta * combination (n2, n) / ((long double) pow4n);
        ldResult += c * (t + sines / pow4n) * pow (parameter, n);
      }

    return ldResult;
  }

The above function still has an error of about 3 parts per million compared to the results from functions.wolfram... one part per quadrillion with theta = pi/2. The latter is simply because Visual C doesn't really use the long double type. It might be the same reason the sine summation gets thrown off too, but I find it hard to believe that the error is so large.

[Kuroneko]

...choose function for real arguments. Is there a good explanation for that? What that means, I mean...

Combinatorially, not in general, although C(x,k) is the coefficient of z^k in the Maclaurin expansion of (1+z)^x. Specific values of x can have some appropriate combinatorial interpretation, but I don't think the general case does. However, |C(n,k)| for integer n<0 can have the interpretation of picking k objects out of |n| boxes, with no restrictions as to how many objects can be picked from each box (other than being non-negative, of course).

[Ariphaos]

Combinatorially, not in general, although C(x,k) is the coefficient of z^k in the Maclaurin expansion of (1+z)^x. Specific values of x can have some appropriate combinatorial interpretation, but I don't think the general case does.

Something special about sqrt (1 - sin^2) - integrals and 1/2 then, or does it get used elsewhere?

I honestly don't even remember how to construct series integrals (as if that wasn't obvious). :-/

However, |C(n,k)| for integer n<0 can have the interpretation of picking k objects out of |n| boxes, with no restrictions as to how many objects can be picked from each box (other than being non-negative, of course).

That sounds like math abuse.

[Kuroneko]

Something special about sqrt (1 - sin^2) - integrals and 1/2 then, or does it get used elsewhere?

It's the binomial theorem, and it's valid for complex numbers, if we write C(x,y) = x!/[y!(x-y)!] and x! = Int_0^\infty[ t^x exp(-t) dt ]. One can verify by integration by parts that x! is the standard factorial for non-negative integers and (x+1)! = (x+1)x! for all complex numbers, if x! is defined. We actually have (x-1)! = Γ(x) as the so-called "gamma function" for this. Unfortunately, it's not defined at nonpositive integers, so we can't use it for those cases (although since we're generally taking ratios, the limits might be well-defined).

That sounds like math abuse.

Not really. We can define the symbol C(n,k) as with the above code, which works for most complex numbers as well. Or if we treat the binomial theorem formally as defining the symbol C(n,k):
(1+z)^x = Sum_{k=0}^\infty C(x,k) z^k,
then C(x,k) exists even for negative integers x and actually complex numbers as well, since the right-hand side is just the MacLaurin series of the left-hand side, and it is well-defined. Now, we have
(1-z)^{-1} = (1 + z + z² + z³ + ... )
as the generating function of the number of ways to pick k objects out of one box, i.e., the coefficient of z^k is that number (which is 1, as it should be). Hence if we have n boxes, the coefficient of z^k in
(1-z)^{-n} = (1 + z + z² + z³ + ... )^n
gives us the number of ways to pick k objects out of n boxes with no restrictions as to how many we can pick from each box. For example, for n = 5 boxes and k = 2 objects, we can pick both out of one box (C(5,1) = 5 choices) or two out of different boxes (C(5,2) = 10 choices); thus the total number is 15. Sure enough,
(1+z)^{-5} = 1+5z+15z²+35z³+...
has the correct coefficient of z². Hence C(-5,2) = (-1)²15 = 15. Similarly, we would have C(-5,3) = -35.

Alternatively, we can say that |C(-n,k)| is the number of ways we can pick k things out of n choices if we're allowed to pick something more than once, but without order. (If with order, the answer is obviously n^k.) We don't have to do this, by the way--|C(-n,k)| = C(n+k-1,k). That's the form that's usually used.

[Ariphaos]

It's the binomial theorem, and it's valid for complex numbers, if we write C(x,y) = x!/[y!(x-y)!] and x! = Int_0^\infty[ t^x exp(-t) dt ]. One can verify by integration by parts that x! is the standard factorial for non-negative integers and (x+1)! = (x+1)x! for all complex numbers, if x! is defined. We actually have (x-1)! = Γ(x) as the so-called "gamma function" for this. Unfortunately, it's not defined at nonpositive integers, so we can't use it for those cases (although since we're generally taking ratios, the limits might be well-defined).

That is insanely cool and yet none of my Calc or higher math teachers ever mentioned that.

Not really. We can define the symbol C(n,k) as with the above code, which works for most complex numbers as well. Or if we treat the binomial theorem formally as defining the symbol C(n,k):
(1+z)^x = Sum_{k=0}^\infty C(x,k) z^k,
then C(x,k) exists even for negative integers x and actually complex numbers as well, since the right-hand side is just the MacLaurin series of the left-hand side, and it is well-defined. Now, we have
(1-z)^{-1} = (1 + z + z² + z³ + ... )
as the generating function of the number of ways to pick k objects out of one box, i.e., the coefficient of z^k is that number (which is 1, as it should be). Hence if we have n boxes, the coefficient of z^k in
(1-z)^{-n} = (1 + z + z² + z³ + ... )^n
gives us the number of ways to pick k objects out of n boxes with no restrictions as to how many we can pick from each box. For example, for n = 5 boxes and k = 2 objects, we can pick both out of one box (C(5,1) = 5 choices) or two out of different boxes (C(5,2) = 10 choices); thus the total number is 15. Sure enough,
(1+z)^{-5} = 1+5z+15z²+35z³+...
has the correct coefficient of z². Hence C(-5,2) = (-1)²15 = 15. Similarly, we would have C(-5,3) = -35.

Alternatively, we can say that |C(-n,k)| is the number of ways we can pick k things out of n choices if we're allowed to pick something more than once, but without order. (If with order, the answer is obviously n^k.) We don't have to do this, by the way--|C(-n,k)| = C(n+k-1,k). That's the form that's usually used.

Much understood, danke.

Back on topic, I tried Mathar's solution and couldn't make it work (I probably simply didn't understand it). Regardless, now that I have a somewhat usable (if not exactly perfect) length function...

Am I to understand that, given your post (and setting a=1): ds = sqrt (1 - ε²cos²v)dv

...that's actually the inversion? I'm not sure I follow. Namely:

transformation z = x+iy = aε cos iw, where w = u+iv. Thus, x = aε cosh u sin u, y = aε sinh u sin v, so that the points on the ellipse share the same u-coordinate acosh(1/ε)

I assume aε is the angular eccentricity, but the iw / iv / u transformations I don't quite understand. I'll have to look at it more when more free time is had later in the week.

[Kuroneko]

Back on topic, I tried Mathar's solution and couldn't make it work (I probably simply didn't understand it). Regardless, now that I have a somewhat usable (if not exactly perfect) length function...

You could simply try Newton's method instead.

I assume aε is the angular eccentricity, ...

No; here aε is just the product of the semimajor axis and the eccentricity, which gives the positions of the foci. (Come to think of it, there's a small typo in that post in this regard--the ε is missing in one place. Apologies.)

Am I to understand that, given your post (and setting a=1): ds = sqrt (1 - ε²cos²v)dv

[img=right]http://home.att.net/~numericana/answer/elliptic.gif[/img] Yes. This is the line element, and its integral gives the arc length. Note that if the angles are measured from the y-axis instead, the cosine turns into sine, and we get the elliptic integral of the second kind. The v-parameter is Mr. Zviman's θ on that site you referenced, but it is not the Euclidean angle. The actual Euclidean angle θ' is related by tan θ' = (b/a)tan v.

... but the iw / iv / u transformations I don't quite understand. I'll have to look at it more when more free time is had later in the week.

Here, it's just an excuse to transform into elliptic coordinates (u,v). I put it like that simply because I'm used to thinking about it in those terms (it's a nice conformal transformation).

[Ariphaos]

You could simply try Newton's method instead.

Mrf, I suppose I should. The error becomes extreme with high eccentricities, although I'm really not dealing with those much yet.

Come to think of it, would be a heck of a lot faster, too. Will work on that next.

Anyway, I don't suppose this:
s = sqrt[1-ε²sin²v]
s^2 = 1 - ε²sin²v
s^2 - ε²sin²v = 1
ε²sin²v = 1 - s^2
sin²v = (1 - s^2)/ε²
sin v = sqrt (1 - s^2)/ε
v = asin (sqrt (1 - s^2)/ε)

...and going from there has a prayer of working?

Yes. This is the line element, and its integral gives the arc length. Note that if the angles are measured from the y-axis instead, the cosine turns into sine, and we get the elliptic integral of the second kind. The v-parameter is Mr. Zviman's θ on that site you referenced, but it is not the Euclidean angle. The actual Euclidean angle θ' is related by tan θ' = (b/a)tan v.

Oh, yes, I got that, the input to the function is:
atan (tan (θ) / b)

Which is different than the wolfram site and I should probably just merge it into the function.

Here, it's just an excuse to transform into elliptic coordinates (u,v). I put it like that simply because I'm used to thinking about it in those terms (it's a nice conformal transformation).

Ahh, figured it was something like that. Thanks.

[Kuroneko]

Anyway, I don't suppose this:
s = sqrt[1-ε²sin²v]
...
v = asin (sqrt (1 - s^2)/ε)
...and going from there has a prayer of working?

Er, not quite: what that gives you is s = Int[ sqrt[1-ε²sin²v] dv ], except v now measured from the minor axis, since what was given was ds in terms of dv. This is exactly what you're doing when you're programming the elliptic integral function above; it's just a derivation of why that integral gives you the correct answer. Also, it allows one to explicitly relate the true Euclidean angle and to the v-parameter (or θ-parameter, that previous site you had called it).

Oh, yes, I got that, the input to the function is:
atan (tan (θ) / b)
Which is different than the wolfram site and I should probably just merge it into the function.

Eh? Wolfram's equation (58) gives t = atan([a/b]tan θ), calling the t what v is if measured from the minor axis. This is exactly what you have, except for a missing 'a'.

[Ariphaos]

Er, not quite: what that gives you is s = Int[ sqrt[1-ε²sin²v] dv ], except v now measured from the minor axis, since what was given was ds in terms of dv. This is exactly what you're doing when you're programming the elliptic integral function above; it's just a derivation of why that integral gives you the correct answer. Also, it allows one to explicitly relate the true Euclidean angle and to the v-parameter (or θ-parameter, that previous site you had called it).

I was hoping I could cheat and get:
dv = asin (sqrt (1 - s^2)/ε)ds
Or something like that and just integrate.

Eh? Wolfram's equation (58) gives t = atan([a/b]tan θ), calling the t what v is if measured from the minor axis. This is exactly what you have, except for a missing 'a'.

I'm referring to this page which I was using to test the validity of my function. My function takes the modified angle and not the Euclidean one.

[Ariphaos]

*snip*

Don't mind my previous post, I must have been tired, and the Newton's Method inversion works perfectly, thanks.

I've expanded the 'infinite sums' part of the series on my own, which is much faster since it's recursive. I should probably look into doing it for the main integral as well.

My error from before was due to the fact that I simply missed this. Is it possible to apply that to the original series?

[Kuroneko]

My error from before was due to the fact that I simply missed this. Is it possible to apply that to the original series?

Hmm... that site seems to be mistaken in its expansion, although as far as I can tell the summation form is correct. It is indeed faster, but I can't see any obvious way to adapt it for the arc length case. Below is a hastily written implementation of both the complete elliptic integral (giving 1/4 arc length of an ellipse with major axis 1) and the Gauss-Kummer series.

The complete EllipticE is very badly behaved for elliptic modulus ~ 1; I haven't got around to making an alternate method near there.

Code: Select all

#define _LD_PI          (3.1415926535897932384626433832795)
#define _LD_PI2         (1.5707963267948966192313216916398)
typedef long double     ldouble;
typedef unsigned int    uint;

/* Rising factorial */
ldouble _ldPochhammerR(ldouble x, uint n) {
  ldouble res = 1.0; uint k = 0;
  while ( k < n ) { res *= x+(ldouble)(k++); }
  return res;
}

/***
 * Hypergeometric series
 *   z      -- argument to evaluate
 *   p, q   -- sizes of a, b arrays, respectively
 *   n1, n2 -- lower and upper limits of sum
 * TODO: save last term to quickly restart computation
 ***/
ldouble _ldHyperF(ldouble z, uint p, uint q, ldouble *a, ldouble *b,
		  uint n1, uint n2) {
  ldouble res = 0.0, term = 1.0;
  uint n, k;

  if ( n2 < n1 ) return 0.0;
  /* Get (n1)st term, or (n1-1)st if last != NULL */
  for ( k = 0; k < p; term *= _ldPochhammerR(a[k++],n1) );
  for ( k = 0; k < q; term /= _ldPochhammerR(b[k++],n1) );
  term /= _ldPochhammerR(1.0,n1); /* (n1)!, but as long double */
  term *= pow(z,(ldouble)(n1)); res = term;
  /* Summation loop */
  for ( n = n1+1; n <= n2; n++ ) {
    term *= z/(ldouble)(n);
    for ( k = 0; k < p; term *= ((ldouble)(n-1)+a[k++]) );
    for ( k = 0; k < q; term /= ((ldouble)(n-1)+b[k++]) );
    res += term;
  }
  return res;
}

/***
 * Complete Elliptic Integral of the Second Kind
 *   (bad convergence for k ~ 1; needs refinement)
 ***/
ldouble _ldEllipticEC(ldouble k, uint n) {
  ldouble a[] = {-0.5,0.5}, b[] = {1.0};
  return _LD_PI2*_ldHyperF(k*k,2,1,a,b,0,n);
}

/***
 * Gauss-Kummer Series
 *  not so good at ellipticity ~ 1; needs valid argument checking
 ***/
ldouble _ldGaussKummer(ldouble a, ldouble b, uint n) {
  ldouble A[] = {-0.5,-0.5}, B[] = {1.0};
  ldouble h = (a-b)/(a+b); h = h*h;
  return _LD_PI*(a+b)*_ldHyperF(h,2,1,A,B,0,n);
}

[Ariphaos]

My own code. Renamed things as appropriate and added the ldouble typedef, since someone might want to plug in a high-precision type in the future.

Visual C (and some other compilers) do not use the long double type since it was intended to be internal to the FPU and people consider it iffy, so my comments are based on double precision. Anyway...

Code: Select all

typedef long double ldouble;

// Simple combination function.  Breaks down for n > ~1026
ldouble combination (unsigned short n, unsigned short k)
  {
    ldouble ldResult = 1;

    for (int i = 0; i < k; i++)
      {
        ldResult = ldResult * (n - i) / (k - i);
      }

    return ldResult;
  }

// Gauss-Kummer expansion for finding the circumference of an ellipse.
// Takes the semiminor axis and assumes semimajor axis = 1.
// Begins to break down for modulus > .96 
// Double precision cannot converge for modulus > ~.999944
ldouble ldGaussKummer (ldouble ldMinor, unsigned short usIterations)
  {
    ldouble ldResult = M_PI*(1 + ldMinor);
    ldouble c = 0, pow4n = 1.0;
    ldouble ldH = ((1.0 - ldMinor)*(1.0 - ldMinor)) / ((1.0 + ldMinor)*(1.0 + ldMinor));
    ldouble ldHp = 1.0;
    unsigned short n = 0, n2 = 0;

    for (n = 1; n < usIterations; n++)
      {
        n2 = n << 1;
        ldHp *= ldH;
        pow4n = pow4n * 4.0;
        c = combination (n2, n) / ((1 - n2)*pow(-4.0, n));
        ldResult += M_PI*(1 + ldMinor) * ldHp * c*c;
      }
    return ldResult;
  }

// Computes an arc length of a given elliptical angle using the incomplete
// integral of the second kind.  ldTheta is the elliptic angle, not the Euclidean
// one.  Passed in via:
// ldTheta = atan (tan (ldAngle) * ldMajor / ldMinor);
// Begins to break down for modulus > .5
// Double precision does not converge for modulus > .97
ldouble ldEllipticE (ldouble ldTheta, ldouble ldModulus, unsigned short usIterations)
  {
    ldouble s = 1;
    ldouble t = 0;
    ldouble c = 0;
    ldouble pow4n = 1.0;
    ldouble ldCosTheta = cos (ldTheta);
    ldouble ldSinTheta = sin (ldTheta);
    ldouble ldResult = ldTheta;
    unsigned short n = 0, n2 = 0;

    for (n = 1; n < usIterations; n++)
      {
        n2 = n << 1;
        pow4n = pow4n * 4.0;
        t = combination (n2, n) / (long double) pow4n;
        c = -1.0 * (t / (n2-1)) * pow (ldModulus, n);
        s = -1.0 * pow (ldSinTheta, n2 - 1) * ldCosTheta / (long double) n2 + ((long double) (n2-1) / (long double) n2) * s;
        ldResult += c * (ldTheta * t + (s - t));
      }

    return ldResult;
  }

// Inverse of the Elliptic function using Newton's Method.
// Limited by the EllipticE function in derives from.  Generally converges within 6 iterations
// for double precision, sometimes faster.
ldouble ldInverseEllipticE (ldouble ldLength, ldouble ldModulus, unsigned short usIterations)
  {
    unsigned short i = 0;
    ldouble ldResult = ldLength;
    ldouble ldSin = 0;

    for (i = 0; i < usIterations; i++)
      {
        ldSin = sin (ldResult);
        ldResult = ldResult - (ldEllipticE (ldResult, ldModulus, 513) - ldLength) /
                   sqrt (1 - ldModulus * ldSin * ldSin );
      }

    return ldResult;
  }

[Ariphaos]

Err forgot to convert the casting to the typedef. Oh well.

Hmm... that site seems to be mistaken in its expansion, although as far as I can tell the summation form is correct. It is indeed faster, but I can't see any obvious way to adapt it for the arc length case. Below is a hastily written implementation of both the complete elliptic integral (giving 1/4 arc length of an ellipse with major axis 1) and the Gauss-Kummer series.

The expansion is correct, I thought it was off too initially.

h = ((1.0 - sqrt (1 - m))*(1.0 - sqrt (1 - m))) / ((1.0 + sqrt (1 - m))*(1.0 + sqrt (1 - m)))

Where m is the modulus... It just seems simple enough that I figured it might be worth trying to see what it did with the individual terms and factoring the changes into the recursive sine integral... if that's at all possible.

Alternatively, I could expand / toy with the recursion and see what patterns can be abused.

[Kuroneko]

You can speed up the calculation of the Gauss-Kummer series a significant amount if you don't re-calculate coefficient each time. The binomial and power functions are relatively costly to call on each iteration, and your particular implementation will be roughly O(n²) time for n as the number of terms. In the Gauss-Kummer series, the initial term is a_0 = 1 and the subseqent terms are given by a_{k+1} = a_k[(n-1/2)/(n+1)]². That will bring it back to O(n).

For example:

Code: Select all

ldouble ldGaussKummer2(ldouble ldMinor, unsigned short usIterations) {
  ldouble ldResult = 1.0, ldNterm = 1.0;
  ldouble ldH = (1.0 - ldMinor)/(1.0 + ldMinor);
  unsigned short n;

  for (n = 1; n < usIterations; n++) {
    ldouble t = ldH*((ldouble)(n)-1.5)/((ldouble)(n));
    ldNterm *= t*t;
    ldResult += ldNterm;
  }
  return M_PI*(1.0+ldMinor)*ldResult;
} 

In real-time, this implementation seems to be about 12:30:50 times faster on my system for n = 16,64,128 iterations, respectively (based on the number of ticks for a loop of 2^16 identical calls with b = 0.99). The ldEllipticE function can be made significantly faster in a similar manner. Try calculating the next term from the previous without combination() or pow().

[Ariphaos]

Dur, yeah, was planning on doing that but it was late and the method wasn't obvious to me at the time. Here's the 'cleaned' EllipticE function:

Code: Select all

ldouble ldEllipticE (ldouble ldTheta, ldouble ldModulus)
  {
    ldouble s = 1.0;
    ldouble ldCosTheta = cos (ldTheta);
    ldouble ldSinTheta = sin (ldTheta);
    ldouble ldSinTerm = ldSinTheta;
    ldouble ldNTerm = 1.0;
    ldouble ldMTerm = 1.0;
    ldouble ldResult = ldTheta, ldPrevResult = 0;
    unsigned short n = 0, n2 = 0, n4 = 0;

    for (n = 1; ldResult != ldPrevResult; n++)
      {
        ldPrevResult = ldResult;
        n2 = n << 1;
        n4 = n2 << 1;
        ldNTerm *= ((ldouble) n4-2.0) / ((ldouble) n4);
        ldMTerm *= ldModulus;
        s = -1.0 * ldSinTerm * ldCosTheta / (ldouble) n2 + ((ldouble) (n2-1) / (ldouble) n2) * s;
        ldResult += -1.0 * (ldNTerm / (n2-1)) * ldMTerm * (ldTheta * ldNTerm + (s - ldNTerm));
        ldSinTerm *= ldSinTheta * ldSinTheta;
      }

    return ldResult;
  }

Still somewhat messy due to having to basically factor s back out each time, but it's ten times faster for Iterations=32 and three hundred times faster for Iterations=128 (measuring 65536 runs by the millisecond) at b=.1 and Euclidean angle=1 (before I had it quit on success).

Thanks.

Still, getting something that would converge on the order of the Gauss-Kummer series would be nicer, if that's even possible. Unfortunately I really have no idea where to start. Or I could just use GK since computing arclengths on elongated ellipses will probably not need to be frequent.