BTW, "elliptic sector" could also reasonably be used to describe a region
having a "vertex" at the center, rather than at a focus, of the ellipse.
But the area of such a region has already been dealt with adequately; see
the end of <http://mathforum.org/dr.math/problems/brady2.28.98.html>, for
example. This article deals with focal (rather than central) sectors, as
described in the previous paragraph.
Here's a response of mine, edited from a recent discussion in another math
newgroup:
-------------------------------------------
"Pal" <o...@nospampal.com> wrote:
> I am trying to calculate the sector area from a foci as in Kepler's 2nd
> Law. I was hoping to be able to do this in polar form.
Fine. With a focus at the pole, an equation for the ellipse is
r = a(1-e^2)/(1-e*cos(t))
where r denotes radius, a semimajor axis length, e eccentricity,
and t angle (in radians, with t=0 giving aphelion).
> So if a planet is at 20* longitude it has a radius of R1. Then it moves
> to 70* longitude with a radius of R2, what is the area swept by these 2
> locations. I would be using the orbital parameters of the planets as
> given by NASA.
Good. You'll be working with a specific planet and will know (or can
calculate from the NASA data) the semiaxes lengths a and b and the
eccentricity e. Since you already know the polar angles corresponding to
R1 and R2, rather than using those radii themselves, we can use the
angles (in radians), which I'll call A and B (with B > A), in the following
calculation.
Define the function F(t) =
t + e*Sqrt(1-e^2)*sin(t)/(1-e*cos(t)) + 2*Arctan(e*sin(t)/(1-e*cos(t)+Sqrt(1-e^2)))
where Arctan denotes the principal-valued inverse tangent function.
Then the area of the sector, swept out about the focus at the pole as the
angle t changes from A to B, is
area = 1/2*a*b [F(B) - F(A)].
I've never seen or derived this result before. I hope I've expressed it as
simply as possible. Surely it must be a standard result. Can any other
reader provide a reference for this or an equivalent result?
-------------------------------------------
Additional notes:
The area formula above is obviously a generalization of the formula
1/2*r*r [B - A] for the area of a sector of a circle.
It is nice that only the principal-valued inverse tangent function is used
above. Some other formulations might require care in dealing with different
branches of the inverse tangent.
Comments and suggestions are welcome.
David Cantrell
This is can be handled by unflattening the ellipse to a circle. See
http://www.whim.org/nebula/math/images/kepler.gif
Offset the ellipse so that the center is at the origin:
( ae + r cos(t), r sin(t) )
1-e^2 1-e^2
= a( e + --------- cos(t), --------- sin(t) )
1+ecos(t) 1+ecos(t)
a
= --------- ( e(1+ecos(t)) + (1-e^2)cos(t), (1-e^2)sin(t) )
1+ecos(t)
a
= --------- ( e+cos(t), (1-e^2)sin(t) ) [1]
1+ecos(t)
Unflatten the ellipse to a circle by stretching in the y direction by
1/sqrt(1-e^2):
a
--------- ( e+cos(t), sqrt(1-e^2)sin(t) ) [2]
1+ecos(t)
We can verify that [2] is a circle by computing its absolute value:
a 2 2 2
--------- sqrt( e + 2ecos(t) + cos (t) + (1-e^2)sin (t) )
1+ecos(t)
a 2 2
= --------- sqrt( 1 + 2ecos(t) + e cos (t) )
1+ecos(t)
a
= --------- ( 1+ecos(t) )
1+ecos(t)
= a [3]
The big trick is one that is used pretty commonly with integration
substitutions which use z = tan(s/2); compute the tangent of half the
angle of [2] using tan(s/2) = y/(a+x), where s is the angle relative to
the origin:
sqrt(1-e^2)sin(t)
tan(s/2) = --------------------
1+ecos(t) + e+cos(t)
sqrt(1-e^2) sin(t)
= ------------ --------
1+e 1+cos(t)
1-e
= sqrt( --- ) tan(t/2) [4]
1+e
Thus, the position on the circle is simply a(cos(s),sin(s)). The angle
s is called the eccentric anomaly. See
http://mathworld.wolfram.com/EccentricAnomaly.html
The area we are looking for is the area swept out by the radius of the
circle (a^2 s/2), minus the triangle formed by the center of the circle,
the focus of the ellipse, and the planet (a^2 e sin(s)/2), all scaled by
sqrt(1-e^2). That is, the area swept out by a line through the focus of
an ellipse from angle 0 to angle t is
a^2 sqrt(1-e^2)
A = --------------- (s - e sin(s)) [5]
2
where s is computed from t as in [4]. The quantity s - e sin(s) is
called the mean anomaly since it has a constant rate of change with
respect to time for planetary motion.
Rob Johnson
r...@whim.org
Sure.
[snip]
OK, let's see how your formulation works in a simple example. (Later we'll
see how my formulation works in the same example.) Let t = 2*pi. To use [5]
we must first find s from [4]. But, with t = 2*pi, [4] merely tells us that
tan(s/2) = 0. So we know that s = 2*arctan(0), where arctan denotes the
multivalued inverse tangent relation. But which of that infinity of
possible values of arctan(0) are we to choose? Of course, Rob, you and I
could easily say "We must choose the branch of arctan so that we get s =
2*pi.", but it would probably not be so easy for a naive person -- much
less, for a computer -- to choose the correct branch. I would guess that
many naive users would conclude, incorrectly, that s = 0, and would then be
puzzled why [5] gives them area = 0 for the full ellipse.
In short, your formulation does not give a solution in closed form, in the
sense that the choice of the branch of arctan is not specified. My
formulation, by contrast, is in closed form, always using the
principal-valued function Arctan. For the example, using my formulation
above with angles B = 2*pi and A = 0, we obtain
area = 1/2*a*b [F(2*pi) - F(0)]
= 1/2*a*b [ 2*pi - 0 ]
= pi*a*b
as required for the area of the entire ellipse.
David Cantrell
I was trying to supply in [5], a simpler looking formula. [4] is a
relation, so using it does allow multiple solutions. Since s goes from
0 to 2pi while t goes from 0 to 2pi, the solution to [4] that we want is
1-e t
s = 2 atan( sqrt( --- ) tan(t/2) ) + 2pi round( --- ) [6]
1+e 2pi
where round(x) = floor(x+1/2). The right summand in [6] is the closest
integer multiple of 2pi to t, which is 0 when t is in (-pi,pi).
Rob Johnson
r...@whim.org
Certainly [6] correctly addresses the concern I'd mentioned.
BTW, for anyone interested in implementing Rob's formulation or mine: It
might _appear_ that the two formulations give different results. In fact,
they give the same result. However, one must be aware that Rob's
formulation takes t=0 at perihelion, whereas mine takes t=0 at aphelion.
Thus, our two formulations are "out of phase" by pi.
Final comment: If t = pi in [6], there is a problem. What would be the
value of Arctan(Sqrt((1-e)/(1+e)) tan(pi/2))? Of course, IMO, there's no
problem assigning a value to tan(pi/2). But if we have naively programmed
[6] in, say, a computer algebra system, what will it give for the
arctangent's value? In Mathematica at least, we get "Indeterminate". OTOH,
my function F(t) avoids any such problems as long as the ellipse is
nondegenerate ( 0 <= e < 1 ).
Regards,
David Cantrell
http://groups.google.com/groups?selm=20030403223940.701$w...@newsreader.com
I only refer to it since it has become quite lengthy.
I see now what is going on. You derived from
1-e
tan(s/2) = sqrt( --- ) tan(t/2) [1]
1+e
that
t-s e sin(t)
tan( --- ) = -------------------------- [2]
2 1 + e cos(t) + sqrt(1-e^2)
so that computing s becomes
e sin(t)
s = t - 2 atan( -------------------------- ) [3]
1 + e cos(t) + sqrt(1-e^2)
And now, s can be used in the equation
a^2 sqrt(1-e^2)
A = --------------- (s - e sin(s)) [4]
2
I agree that [3] is a better formula for s than
1-e t
s = 2 atan( sqrt( --- ) tan(t/2) ) + 2pi round( --- ) [5]
1+e 2pi
since it avoids the problem when tan(t/2) blows up, which is where
round(t/2pi) is discontinuous (using 2pi as a single constant).
It should now be the case the [3] and [4] above are equivalent to your
formula, accounting for our phase difference.
Rob Johnson
r...@whim.org
I have collected a paper I wrote a long time ago about Kepler's Laws and
parts of this thread into a couple of web pages. If you are interested,
please see http://www.whim.org/nebula/math/kepler.html and its link to
"Planetary Orbits".
Rob Johnson
r...@whim.org