[FIXED] Plot illuminated percentage of the moon

Issue

Its easy to get the percentage of the moon illuminated on a given day

import ephem
a=ephem.Moon(datetime.now())
mn=a.moon_phase #illuminated percentage of the moon.

Now I would like to try and plot this to represent the moon phases. But I can't get an idea of how to fill only a percentage of a circle to simulate the earth's shadow on the moon.

My closest try is this one:

First plot the "full" moon

fig, ax = plt.subplots()


theta=np.linspace(0, 2*np.pi, 100)
r = np.sqrt(1)
c1x = r*np.cos(theta)
c1y = r*np.sin(theta)
ax.plot(c1x, c1y,color='k',lw=2,zorder=0*a)

Then add a bigger black circle to simulate the shadow passing over the moon. the offset "xmn" would move the shadow to the side so that only "n" percentage of the full moon is visible

r = np.sqrt(2)
c2x = r*np.cos(theta)
c2y = r*np.sin(theta)
xmn=1.2
ax.fill(c2x+xmn, c2y,color='k',lw=2,zorder=1*a)

Finally hide everything outside the "full moon" circle

#Circle to hide all circles
theta=np.linspace(0, 2*np.pi, 100)
# the radius of the circle
r = np.sqrt(4)
c3x = r*np.cos(theta)
c3y = r*np.sin(theta)
ax.plot(c3x, c3y,color='w')

ax.fill_between(c3x,c3y,c1y,color='w',zorder=2)
ax.fill_betweenx(c1x,c1y,c3y,color='w',zorder=2)

ax.set_xlim(-1.1, 1.1)  
ax.set_ylim(-1.1, 1.1)

plt.show()

I can't find a smart way to do this. How can I find xmn?

Thanks in advance

Solution

Eureka!

I found the answer on this site

The way they solve the problem is by calculating the area of the overlap between the two circles, and then subtract it to the area of the moon circle. For that, the distance between the centre of the two circles must be known, which is exactly what I needed (xmn).

For that they calculate the area of the crescent of the bigger circle, then find the area of the overlap between the two, and finally the area of the crescent of the smaller circle. Then we can find the percentage of the total area.

Somehow, the equation to find the area of the crescent of the bigger circle (lune_1 on the website), gives me the area of the smaller crescent instead. Don't really understand why but it fits my purpose.

If only I still had the math chops to pull it off, solving equation lune_1 in order to find "d" would give me a direct answer, in which lune_1= illuminated area and d=distance between the centre of the two circles.

After a few adjustments, I ended up with this code

        def calc_crescent_pos(r1,r2,mn):
            # mn = Illuminated percentage of the moon
            # r1 = radius of the moon circle
            # r2 = radius of the shadow circle

            lune=0
            d=r2-r1
            area=np.pi*r1**2 # area of the moon circle
            while lune/area < mn: #increase separation between the 2 circles, until the area of the moon crescent is equal to the illuminated percentage of the moon
                d=d+0.01
                lune = 2*(np.sqrt( (r1+r2+d) * (r2+d-r1) * (d+r1-r2) * (r1+r2-d))/ 4) + r1**2 * math.acos( (r2**2-r1**2-d**2) / (2*r1*d) ) - r2**2 * math.acos( (r2**2+d**2-r1**2) / (2*r2*d))                                     
                 
            return d-0.01

        import ephem,math
        import matplotlib.pyplot as plt
        import numpy as np
        from datetime import *

        a=ephem.Moon(datetime.now())
        mn=a.moon_phase #illuminated percentage of the moon.

        fig, ax = plt.subplots()

        #Circle1; full moon
        theta=np.linspace(0, 2*np.pi, 100)
        # the radius of the circle
        r1 = np.sqrt(1)
        c1x = r1*np.cos(theta)
        c1y = r1*np.sin(theta)
        ax.plot(c1x, c1y,color='k',lw=2,zorder=0*a)
    
        #Circle 2; Shadow
        r2 = np.sqrt(2)
        c2x = r2*np.cos(theta)
        c2y = r2*np.sin(theta)

        xmn=calc_crescent_pos(r1,r2,mn)
        
        ax.fill(c2x+xmn, c2y,color='k',lw=2,zorder=1*a)

        #Circle to hide all circles
        theta=np.linspace(0, 2*np.pi, 100)
        # the radius of the circle
        r = np.sqrt(4)
        c3x = r*np.cos(theta)
        c3y = r*np.sin(theta)
        ax.plot(c3x, c3y,color='w')

        ax.fill_between(c3x,c3y,c1y,color='w',zorder=2)
        ax.fill_betweenx(c1x,c1y,c3y,color='w',zorder=2)

        ax.set_xlim(-1.1, 1.1)  
        ax.set_ylim(-1.1, 1.1)

        plt.show()

Now its only a matter of adding a few lines to make sure the shadow goes the right way whenever the moon is waning or crescent.

Answered By - Elcook