eclipse.mws
Area of an Eclipse
by Jason Schattman, Waterloo Maple, jschattman@maplesoft.com
During an eclipse of one planet by another, how does the visible area of the eclipsed planet change, as viewed from earth? In this Maple exploration, we compute the area in terms of the apparent radii of the planets,
r
and
R
, and the apparent distance
d
between their centers.
> |
restart; r := 1: R := 2:
|
> |
for i from 1 to 50 do
c1 := plottools[disk]([0,0], R, color=green):
c2 := plottools[disk]([-.9*(R-r)-2*r*i/50,0], r, color=red):
frame[i] := plots[display](c2,c1):
end do:
plots[display](seq(frame[i],i=1..50), insequence=true, scaling=constrained, axes=none);
|
In terms of the picture above, we seek the area of the visible green region as a function of
d
,
R
and
r
.
We first tell Maple our assumptions explicitly. In particular,
,
, and
d
is in
> |
assume( r > 0, R >= r, d <= R+r, d >= R-r ): interface(showassumed=0);
|
Let's rotate the diagram above by 90 degrees so that we can treat the intersecting arcs as functions of
x
and find the area by integration. We consider separately the cases that
d
is greater or less than
. These cases determine whether we use the upper or lower arc of the red circle for the integration.
In Case 1, shown below,
. We integrate between the green arc and the lower arc of the red circle. Then the area of the eclipse is
> |
a1 := plottools[circle]([0,0], 1, color=red, thickness=2):
a2 := plottools[arc]([0,1.5], 2, 5*Pi/4..7*Pi/4, color=green, thickness=2):
plots[display](a1,a2, scaling=constrained, title="Case 1");
|
In the second case, shown below,
. We integrate between the green arc and the upper arc of the red circle. Then the area of the eclipse is then
> |
a1 := plottools[circle]([0,0], 1, color=red, thickness=2):
a2 := plottools[arc]([0,1.9], 2, 5*Pi/4..7*Pi/4, color=green, thickness=2):
plots[display](a1,a2, scaling=constrained, axes=normal, title="Case 2");
|
We'll work out case 2 first. We can express the lower arc of the green circle as
and the upper arc of the red circle as
.
> |
y1 := d - sqrt(R^2-x^2);
y2 := sqrt(r^2-x^2);
|
Now we let Maple compute the intersection points of the circles by solving
y1 = y2
for
x
. These are the end points of the integration.
> |
x1,x2 := solve( y1=y2, x );
|
The antiderivative of the difference between the curves.
> |
Int( y2 - y1, x ) = int( y2 - y1, x );
|
For Maple-convenience, we make an explicit function
g(x)
out of the antiderivative using Maple's
unapply
command..
> |
g := unapply(rhs(%), x):
|
By symmetry about the y-axis, the area is double the definite integral from 0 to the right-most intersection point. Therefore, the formula for the visible area in Case 2 is
.
> |
Pi*R^2 - 2*(g(x1)-g(0));
|
So that we can easily plot this formula, we again make the area formula an explicit function of
R, r
and
d
using the
unapply
command.
> |
VA1 := unapply( %, R, r, d ):
|
Now we go back to case 1. Here we represent the lower arc of the circle by the function
Antiderivative of the difference between the curves.
> |
Int( y1 - y2, x ) = int( y1 - y2, x );
|
Again, we make an explicit function
g(x)
out of the antiderivative above using Maple's
unapply
command..
> |
g := unapply(rhs(%), x):
|
By symmetry about the y-axis, the area is double the definite integral from 0 to the right-most intersection point. Therefore, the formula for the visible area in Case 1 is
.
> |
Pi*R^2 - (Pi*r^2 - 2*(g(x1)-g(0)));
|
> |
VA2 := unapply(%, R, r, d):
|
Finally, we make the visible area a two-piece function of d, where the breakpoint occurs at
> |
VisibleArea := (R,r,d)->piecewise( d >= sqrt(R^2-r^2),
VA1(R,r,d),
d < sqrt(R^2-r^2),
VA2(R,r,d)):
|
The area when
R = r
, that is, the planets have the same apparent size.
Two quick reality checks. When
R=r
and
d=2*R,
the red circle no longer obscures the green one, so the area of the eclipse should just be the area of green circle.
When
d = R-r
, the red circle is completely contained in the green one, so the area of the eclipse should just be the difference of the circles' areas.
> |
simplify( VA2(R, r, R-r) );
|
Plot of the area when
R = r = 1
> |
plot( VisibleArea(1, 1, d), d=0..2, thickness=2, color=blue);
|
The area when R > r.
> |
plot( [VisibleArea(1, 1, d), VisibleArea(1, .95, d)], d=0..2, thickness=3, color=[blue,orange]);
|
The cases
r = 1, 1.125, 1.25, 1.375
and
1.5
.
> |
plot( [seq(VisibleArea(1+i/8, 1, d), i=0..4)], d=0..3, thickness=2);
|