Eureka ！
我在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()

现在只需要增加几条线，以确保阴影走正确的道路时，月亮是减弱或新月。