Application Center - Maplesoft

App Preview:

Area of an Eclipse

You can switch back to the summary page by clicking here.

Learn about Maple
Download Application


 

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);

[Maple Plot]

In terms of the picture above, we seek the area of the visible green region as a function of d , R  and r .

>    restart;

We first tell Maple our assumptions explicitly.  In particular, 0 < r , r <= R , and d  is in [R-r, R+r]

>    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   sqrt(R^2-r^2) .  These cases determine whether we use the upper or lower arc of the red circle for the integration.

In Case 1, shown below, d < sqrt(R^2-r^2) .  We integrate between the green arc and the lower arc of the red circle.  Then the area of the eclipse is Pi*R^2-(Pi*r^2-integrated_area)

>    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");

[Maple Plot]

In the second case, shown below, sqrt(R^2-r^2) <= d .  We integrate between the green arc and the upper arc of the red circle.  Then the area of the eclipse is then Pi*R^2-integrated_area

>    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");

[Maple Plot]

We'll work out case 2 first.  We can express the lower arc of the green circle as y1 = d-sqrt(R^2-x^2)  and the upper arc of the red circle as y2 = sqrt(r^2-x^2) .  

>    y1 := d - sqrt(R^2-x^2);
y2 := sqrt(r^2-x^2);

y1 := d-(R^2-x^2)^(1/2)

y2 := (r^2-x^2)^(1/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 );

x1, x2 := 1/2*(2*R^2*d^2-d^4-r^4+2*r^2*d^2+2*r^2*R^2-R^4)^(1/2)/d, -1/2*(2*R^2*d^2-d^4-r^4+2*r^2*d^2+2*r^2*R^2-R^4)^(1/2)/d

The antiderivative of  the difference between the curves.

>    Int( y2 - y1, x ) = int( y2 - y1, x );

Int((r^2-x^2)^(1/2)-d+(R^2-x^2)^(1/2),x) = 1/2*x*(r^2-x^2)^(1/2)+1/2*r^2*arcsin(1/r*x)-d*x+1/2*x*(R^2-x^2)^(1/2)+1/2*R^2*arcsin(1/R*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)) .  

>    Pi*R^2 - 2*(g(x1)-g(0));

Pi*R^2-1/2*(2*R^2*d^2-d^4-r^4+2*r^2*d^2+2*r^2*R^2-R^4)^(1/2)/d*(r^2-1/4*(2*R^2*d^2-d^4-r^4+2*r^2*d^2+2*r^2*R^2-R^4)/d^2)^(1/2)-r^2*arcsin(1/2/r*(2*R^2*d^2-d^4-r^4+2*r^2*d^2+2*r^2*R^2-R^4)^(1/2)/d)+(2*R...
Pi*R^2-1/2*(2*R^2*d^2-d^4-r^4+2*r^2*d^2+2*r^2*R^2-R^4)^(1/2)/d*(r^2-1/4*(2*R^2*d^2-d^4-r^4+2*r^2*d^2+2*r^2*R^2-R^4)/d^2)^(1/2)-r^2*arcsin(1/2/r*(2*R^2*d^2-d^4-r^4+2*r^2*d^2+2*r^2*R^2-R^4)^(1/2)/d)+(2*R...
Pi*R^2-1/2*(2*R^2*d^2-d^4-r^4+2*r^2*d^2+2*r^2*R^2-R^4)^(1/2)/d*(r^2-1/4*(2*R^2*d^2-d^4-r^4+2*r^2*d^2+2*r^2*R^2-R^4)/d^2)^(1/2)-r^2*arcsin(1/2/r*(2*R^2*d^2-d^4-r^4+2*r^2*d^2+2*r^2*R^2-R^4)^(1/2)/d)+(2*R...
Pi*R^2-1/2*(2*R^2*d^2-d^4-r^4+2*r^2*d^2+2*r^2*R^2-R^4)^(1/2)/d*(r^2-1/4*(2*R^2*d^2-d^4-r^4+2*r^2*d^2+2*r^2*R^2-R^4)/d^2)^(1/2)-r^2*arcsin(1/2/r*(2*R^2*d^2-d^4-r^4+2*r^2*d^2+2*r^2*R^2-R^4)^(1/2)/d)+(2*R...

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 y2 = -sqrt(r^2-x^2)

>    y2 := -sqrt(r^2-x^2);

y2 := -(r^2-x^2)^(1/2)

Antiderivative of the difference between the curves.

>    Int( y1 - y2, x ) = int( y1 - y2, x );

Int(d-(R^2-x^2)^(1/2)+(r^2-x^2)^(1/2),x) = d*x-1/2*x*(R^2-x^2)^(1/2)-1/2*R^2*arcsin(1/R*x)+1/2*x*(r^2-x^2)^(1/2)+1/2*r^2*arcsin(1/r*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))) .  

>    Pi*R^2 - (Pi*r^2 - 2*(g(x1)-g(0)));

Pi*R^2-Pi*r^2+(2*R^2*d^2-d^4-r^4+2*r^2*d^2+2*r^2*R^2-R^4)^(1/2)-1/2*(2*R^2*d^2-d^4-r^4+2*r^2*d^2+2*r^2*R^2-R^4)^(1/2)/d*(R^2-1/4*(2*R^2*d^2-d^4-r^4+2*r^2*d^2+2*r^2*R^2-R^4)/d^2)^(1/2)-R^2*arcsin(1/2/R*...
Pi*R^2-Pi*r^2+(2*R^2*d^2-d^4-r^4+2*r^2*d^2+2*r^2*R^2-R^4)^(1/2)-1/2*(2*R^2*d^2-d^4-r^4+2*r^2*d^2+2*r^2*R^2-R^4)^(1/2)/d*(R^2-1/4*(2*R^2*d^2-d^4-r^4+2*r^2*d^2+2*r^2*R^2-R^4)/d^2)^(1/2)-R^2*arcsin(1/2/R*...
Pi*R^2-Pi*r^2+(2*R^2*d^2-d^4-r^4+2*r^2*d^2+2*r^2*R^2-R^4)^(1/2)-1/2*(2*R^2*d^2-d^4-r^4+2*r^2*d^2+2*r^2*R^2-R^4)^(1/2)/d*(R^2-1/4*(2*R^2*d^2-d^4-r^4+2*r^2*d^2+2*r^2*R^2-R^4)/d^2)^(1/2)-R^2*arcsin(1/2/R*...
Pi*R^2-Pi*r^2+(2*R^2*d^2-d^4-r^4+2*r^2*d^2+2*r^2*R^2-R^4)^(1/2)-1/2*(2*R^2*d^2-d^4-r^4+2*r^2*d^2+2*r^2*R^2-R^4)^(1/2)/d*(R^2-1/4*(2*R^2*d^2-d^4-r^4+2*r^2*d^2+2*r^2*R^2-R^4)/d^2)^(1/2)-R^2*arcsin(1/2/R*...
Pi*R^2-Pi*r^2+(2*R^2*d^2-d^4-r^4+2*r^2*d^2+2*r^2*R^2-R^4)^(1/2)-1/2*(2*R^2*d^2-d^4-r^4+2*r^2*d^2+2*r^2*R^2-R^4)^(1/2)/d*(R^2-1/4*(2*R^2*d^2-d^4-r^4+2*r^2*d^2+2*r^2*R^2-R^4)/d^2)^(1/2)-R^2*arcsin(1/2/R*...

>    VA2 := unapply(%, R, r, d):

Finally, we make the visible area a two-piece function of d, where the breakpoint occurs at d = sqrt(R^2-r^2)

>    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.

>    VisibleArea(R, R, d);

Pi*R^2-(4*R^2*d^2-d^4)^(1/2)/d*(R^2-1/4*(4*R^2*d^2-d^4)/d^2)^(1/2)-2*R^2*arcsin(1/2/R*(4*R^2*d^2-d^4)^(1/2)/d)+(4*R^2*d^2-d^4)^(1/2)

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.

>    VA1( R, R, R+R );

Pi*R^2

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) );

Pi*R^2-Pi*r^2

Plot of the area when R = r = 1

>    plot( VisibleArea(1, 1, d), d=0..2, thickness=2, color=blue);

[Maple Plot]

The area when R > r.

>    VisibleArea(R, r, d);

PIECEWISE([Pi*R^2-1/2*(2*R^2*d^2-d^4-r^4+2*r^2*d^2+2*r^2*R^2-R^4)^(1/2)/d*(r^2-1/4*(2*R^2*d^2-d^4-r^4+2*r^2*d^2+2*r^2*R^2-R^4)/d^2)^(1/2)-r^2*arcsin(1/2/r*(2*R^2*d^2-d^4-r^4+2*r^2*d^2+2*r^2*R^2-R^4)^(1...
PIECEWISE([Pi*R^2-1/2*(2*R^2*d^2-d^4-r^4+2*r^2*d^2+2*r^2*R^2-R^4)^(1/2)/d*(r^2-1/4*(2*R^2*d^2-d^4-r^4+2*r^2*d^2+2*r^2*R^2-R^4)/d^2)^(1/2)-r^2*arcsin(1/2/r*(2*R^2*d^2-d^4-r^4+2*r^2*d^2+2*r^2*R^2-R^4)^(1...
PIECEWISE([Pi*R^2-1/2*(2*R^2*d^2-d^4-r^4+2*r^2*d^2+2*r^2*R^2-R^4)^(1/2)/d*(r^2-1/4*(2*R^2*d^2-d^4-r^4+2*r^2*d^2+2*r^2*R^2-R^4)/d^2)^(1/2)-r^2*arcsin(1/2/r*(2*R^2*d^2-d^4-r^4+2*r^2*d^2+2*r^2*R^2-R^4)^(1...
PIECEWISE([Pi*R^2-1/2*(2*R^2*d^2-d^4-r^4+2*r^2*d^2+2*r^2*R^2-R^4)^(1/2)/d*(r^2-1/4*(2*R^2*d^2-d^4-r^4+2*r^2*d^2+2*r^2*R^2-R^4)/d^2)^(1/2)-r^2*arcsin(1/2/r*(2*R^2*d^2-d^4-r^4+2*r^2*d^2+2*r^2*R^2-R^4)^(1...
PIECEWISE([Pi*R^2-1/2*(2*R^2*d^2-d^4-r^4+2*r^2*d^2+2*r^2*R^2-R^4)^(1/2)/d*(r^2-1/4*(2*R^2*d^2-d^4-r^4+2*r^2*d^2+2*r^2*R^2-R^4)/d^2)^(1/2)-r^2*arcsin(1/2/r*(2*R^2*d^2-d^4-r^4+2*r^2*d^2+2*r^2*R^2-R^4)^(1...
PIECEWISE([Pi*R^2-1/2*(2*R^2*d^2-d^4-r^4+2*r^2*d^2+2*r^2*R^2-R^4)^(1/2)/d*(r^2-1/4*(2*R^2*d^2-d^4-r^4+2*r^2*d^2+2*r^2*R^2-R^4)/d^2)^(1/2)-r^2*arcsin(1/2/r*(2*R^2*d^2-d^4-r^4+2*r^2*d^2+2*r^2*R^2-R^4)^(1...
PIECEWISE([Pi*R^2-1/2*(2*R^2*d^2-d^4-r^4+2*r^2*d^2+2*r^2*R^2-R^4)^(1/2)/d*(r^2-1/4*(2*R^2*d^2-d^4-r^4+2*r^2*d^2+2*r^2*R^2-R^4)/d^2)^(1/2)-r^2*arcsin(1/2/r*(2*R^2*d^2-d^4-r^4+2*r^2*d^2+2*r^2*R^2-R^4)^(1...
PIECEWISE([Pi*R^2-1/2*(2*R^2*d^2-d^4-r^4+2*r^2*d^2+2*r^2*R^2-R^4)^(1/2)/d*(r^2-1/4*(2*R^2*d^2-d^4-r^4+2*r^2*d^2+2*r^2*R^2-R^4)/d^2)^(1/2)-r^2*arcsin(1/2/r*(2*R^2*d^2-d^4-r^4+2*r^2*d^2+2*r^2*R^2-R^4)^(1...
PIECEWISE([Pi*R^2-1/2*(2*R^2*d^2-d^4-r^4+2*r^2*d^2+2*r^2*R^2-R^4)^(1/2)/d*(r^2-1/4*(2*R^2*d^2-d^4-r^4+2*r^2*d^2+2*r^2*R^2-R^4)/d^2)^(1/2)-r^2*arcsin(1/2/r*(2*R^2*d^2-d^4-r^4+2*r^2*d^2+2*r^2*R^2-R^4)^(1...

>    plot( [VisibleArea(1, 1, d), VisibleArea(1, .95, d)], d=0..2, thickness=3, color=[blue,orange]);

[Maple Plot]

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);

[Maple Plot]

>   

>   

>