绘制月球的照明百分比

Question

很容易得到某一天被照亮的月亮百分比

import ephem
a=ephem.Moon(datetime.now())
mn=a.moon_phase #illuminated percentage of the 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)

然后添加一个更大的黑色圆圈来模拟经过月球的阴影。偏移量“xmn”会将阴影移到一边，这样只有“n”百分比的满月可见

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)

终于把“满月”圈外的一切都隐藏了

#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()

我找不到一个聪明的方法来做到这一点。如何找到 xmn？

提前致谢

Answer 1

尤里卡！

我在this site找到了答案

他们解决问题的方法是计算两个圆重叠的面积，然后减去月球圆的面积。为此，必须知道两个圆的中心之间的距离，这正是我所需要的（xmn）。

为此，他们计算大圆的月牙面积，然后找到两者重叠的面积，最后计算出小圆的月牙面积。然后我们可以找到总面积的百分比。

不知何故，找到较大圆圈（网站上的 lune_1）的新月面积的等式给了我较小的新月面积。不太明白为什么，但它符合我的目的。

如果我仍然有数学能力来完成它，求解方程 lune_1 以找到“d”会给我一个直接的答案，其中 lune_1 = 照明区域，d = 两个圆中心之间的距离.

经过一些调整，我最终得到了这段代码

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

现在只需添加几行即可确保在月亏或新月时阴影移动正确。