Project Euler #246：切线与切线之间的椭圆角

Question

我使用原点移位，试图解决problem

My Similar question on StackOverflow

现在原点移位后的椭圆有方程(x/a)^2+(y/b)^2=1 所以在 与椭圆相切的条件 (am)^2+b^2=c^2 其中c = y0-mx0 求解m 其中D = (x0, y0) 是变量点仅在+象限中，因为椭圆是关于轴对称的格点数对于切线之间的角度大于给定角度的情况，如果点 D 位于其他轴上，则将乘以 2 且为 4

但现在我无法找到预期的角度 α，该对比程序可能会导致 (180-α) 或 (90-α ) 或其他任何东西，但并不总是提供 α

输入

64817 64819 11420
3
30 1

输出

1.5 1.118033988749895      #values of a, b of ellipse
88.09084756700362          #Given angle
0 2 84.26082952273322      # x, y, and angle
0 3 56.63298703076825
0 4 42.667925494108374
0 5 34.21605113129826
0 6 28.552637360182302
1 1 132.1304147613666
1 2 105.37675990405205
1 3 126.93618989341962
1 4 138.99680040248643
1 5 146.68139315999935
1 6 151.98229919837212
2 0 99.59406822686043
2 1 70.40251628700852
2 2 123.92838964927597      #In figure it is 56.0716 while it's supplementary angle
2 3 135.0795775738836
2 4 143.16521634758482
2 5 149.05461064130475
2 6 153.44570404972666
3 0 133.43253655778977
3 1 45.1685179704227
3 2 138.68062832170193
3 3 143.5441803020086
3 4 148.1966760890815
3 5 152.19282398152126
3 6 155.49812437089727
4 0 146.4426902380793
4 1 33.150772763942314
4 2 148.18918752912822
4 3 150.3525987174534
4 4 152.8971875924072
4 5 155.43753093848753
4 6 157.77424382459193
5 0 153.61567025059207
5 1 26.21616761276448
5 2 154.36403775389698
5 3 155.4135197186744
5 4 156.82314264518547
5 5 158.40584822123196
5 6 160.005114518988
6 0 158.21321070173815
6 1 21.70150969841879
6 2 158.59334015511746
6 3 159.1528282819498
6 4 159.9664352057684
6 5 160.96054819665326
6 6 162.04466838669927

代码

from math import atan
from math import acos
from math import degrees
from math import atan2
from math import sqrt
from operator import gt
from operator import abs
class P_test():
    def __init__(self,x,y,a,b):
        self.x, self.y, self.a, self.b = x, y, a, b
    
    def Region(self):                # Tests if point(x,y) is out of region of ellipse
        if (self.x/self.a)**2 + (self.y/self.b)**2 >1:
            return 1
    
    def T_angle(self):              # Returns angle between tangents but not expected None
        x, y = self.x, self.y
        A, B , C = self.a**2-self.x**2, 2*self.x*self.y, self.b**2-self.y**2
        D = sqrt(B**2-4*A*C)
        m1, m2 = (-B+D)/(2*A), (-B-D)/(2*A)
        alpha1 = degrees(acos(1/(sqrt(1+m1**2))))        
        alpha2 = degrees(acos(1/(sqrt(1+m2**2))))
        if x*y: return 180-degrees(abs(atan2(y,y/m1)-atan2(y,y/m2)))
        return 180-abs(alpha2)-abs(alpha1)
    
class Datasource():
    def __init__(self, cordinate, radius, pq):
        self.xy, self.r, self.pq = cordinate, radius, pq
    
    def a_b_of_ellipse(self):
        x1, x2, y, r = self.xy[0], self.xy[1], self.xy[2], self.r
        H, K = x1 + (x2-x1)/2, y
        mx, my = x1-H, y-K
        gx, gy = x2-H, y-K
        ae = (mx-gx)/2
        a =gx + (r-(x2-x1))/2
        b = sqrt(a**2-ae**2)
        return a, b
    
    def G_angle(self):
        return degrees(atan(self.pq[0]/self.pq[1]))
    
    
    
xy_cord = list(map(int, input().split()))
radius = int(input())
pq = list(map(int, input().split()))
r1 = Datasource(xy_cord, radius, pq)
aa, bb, G = r1.a_b_of_ellipse()[0], r1.a_b_of_ellipse()[1], r1.G_angle()
print(aa, bb); print(G)
for x in range(0,7):
    for y in range(0,7):
        r2 = P_test(x,y,aa,bb)
        if r2.Region():
            print(x,y,r2.T_angle())

我怎样才能找到只有 Alpha 角

Answer 1

但您可以计算切点，因此无需使用斜率 - 只需从三角形 TP1-D-TP2 获取 alpha

Answer 2

from math import atan, acos, degrees, atan2, sqrt
from operator import gt, abs
class P_test():
    def __init__(self, x, y, a,b):
        self.x, self.y, self.a, self.b = x, y, a, b
    
    def Region(self):
        if self.x**2/self.a**2 + self.y**2/self.b**2 >1:
            return 1
        
    def Slope(self):
        x, y = self.x, self.y
        A, B, C = self.a**2-self.x**2, 2*self.x*self.y, self.b**2-self.y**2
        D = sqrt(B**2-4*A*C)
        m1, m2 = (-B+D)/(2*A), (-B-D)/(2*A)
        return m1, m2
        
class Datasource():
    def __init__(self, cordinate, radius, pq):
        self.xy, self.r, self.pq = cordinate, radius, pq

    def a_b_of_ellipse(self):
        x1, x2, y, r = self.xy[0], self.xy[1], self.xy[2], self.r
        H, K = x1+(x2-x1)/2, y
        mx, my = x1-H, K
        gx, gy = x2-H, y-K
        ae = (mx-gx)/2
        a = gx + (r-(x2-x1))/2
        b = sqrt(a**2-ae**2)
        return a, b

    def G_angle(self):
        return degrees(atan(self.pq[0]/self.pq[1]))
    
    
def Touchpoint(x,y,m1,m2,a,b):
    c1,c2 = y-m1*x, y-m2*x
    return -a**2*m1/c1, b**2/c1, -a**2*m2/c2, b**2/c2

xy_cord = list(map(int, input().split()))
radius = int(input())
pq = list(map(int, input().split()))
r1 = Datasource(xy_cord, radius, pq)
aa, bb, G = r1.a_b_of_ellipse()[0], r1.a_b_of_ellipse()[1], r1.G_angle()
#print(aa,bb)
#print(G)
count = 0
for x in range(0, 100000):
    for y in range(0, 100000):
        r2 = P_test(x,y,aa,bb)
        if r2.Region():
            slope1, slope2 = r2.Slope()[0], r2.Slope()[1]
            #if r2.T_angle()>=G:
            Z = Touchpoint(x,y,slope1,slope2, aa, bb)
            px1,py1,px2,py2 = Z[0], Z[1], Z[2], Z[3]
            M = [degrees(atan2((y-py1),(x-px1))), degrees(atan2((y-py2),(x-px2)))]
            if abs(M[1]-M[0])>=G:
            #    print(x,y, abs(M[1]-M[0]))
                if x*y==0:
                    count += 2
                else:
                    count +=4
            else: break
print(count)