欧拉计划 #246：切线与椭圆的切线之间的夹角
Project Euler #246: Tangents to an ellipse Angle between Tangents
我使用原点偏移并试图解决problem
现在原点移动后的椭圆有方程 (x/a)^2+(y/b)^2=1
所以在椭圆相切条件 (am)^2+b^2=c^2 where c = y0-mx0 solved for m where D = (x0, y0) is variable point in only +象限，因为椭圆关于轴对称，如果点 D 在任一轴上，则切线之间的角度大于给定角度的格点数将乘以 2，否则为 4
但现在我找不到预期的角度 α，这个对比程序可能会导致 (180-α)或 (90-α) 或其他但不给出 α always
输入
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 角度
但是您可以计算切点，因此无需使用斜率 - 只需从三角形 TP1-D-TP2 获取 alpha
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)
我使用原点偏移并试图解决problem
现在原点移动后的椭圆有方程 (x/a)^2+(y/b)^2=1
所以在椭圆相切条件 (am)^2+b^2=c^2 where c = y0-mx0 solved for m where D = (x0, y0) is variable point in only +象限，因为椭圆关于轴对称，如果点 D 在任一轴上，则切线之间的角度大于给定角度的格点数将乘以 2，否则为 4
但现在我找不到预期的角度 α，这个对比程序可能会导致 (180-α)或 (90-α) 或其他但不给出 α always
输入
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 角度
但是您可以计算切点，因此无需使用斜率 - 只需从三角形 TP1-D-TP2 获取 alpha
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)