Как найти точку пересечения прямой линии и кривого набора точек, пожалуйста? Я обработал в соответствии с этот вопрос , однако результат неверен. Где ошибка, пожалуйста? Важная часть обозначена #### в коде. Остальное - для генерации списка для конструирования кода и построения результата.
import math
import numpy as np
import matplotlib.pyplot as plt
G = 6.67430e-11
sun_mass = 1.9891e+30
parsec = 3.086e+16
myr_sec = 3.15576e+13
year_sec = 365.25*24*3600
AU = 149597870700
c = 299792458
kmh_ms = 3.6
mBH = 30
f = 1e-2
n = 1e+4 # pc^{-3}
v = 10 # km s^{-1}
a_BBH = 100 # AU
mass_BHs = 10 # sun mass
r_min = n**(1/3) # pc
sigma_foc = np.pi*(r_min*parsec)**2*(1+2*G*2*mass_BHs*sun_mass/(r_min*parsec*(v/3.6)**2))
T_foc = (n/(parsec**3)*sigma_foc*v/3.6)**(-1)
t_mrg_e3_list = []
pos_e = np.linspace(1e-5,0.9999999999,100000)
for e in pos_e:
a = a_BBH
bracket = (1+73/24*e**2+37/96*e**4)**(-1)
tsec = 5/64*(c**5*(a*AU)**4*(1-e**2)**(7/2))/(G**3*(mBH*sun_mass)**2*2*mBH*sun_mass)*bracket
tmrg = tsec/year_sec
t_mrg_e3_list.append(tmrg)
e_1_diff = [1-eccen for eccen in pos_e]
e_1 = np.array(e_1_diff)
t_mrg_e3 = np.array(t_mrg_e3_list)
fig, ax = plt.subplots()
####
gradient = (T_foc - ax.get_xlim()[1]) / (T_foc - ax.get_xlim()[0])
intercept = ax.get_xlim()[1] - gradient * ax.get_xlim()[0]
dis_y_line = (intercept + e_1* gradient) - t_mrg_e3
ix = np.where(dis_y_line[1:] * dis_y_line[:-1] < 0)[0] # index of points where the next point is on the other side of the line
d_ratio = dis_y_line[ix] / (dis_y_line[ix] - dis_y_line[ix + 1]) # similar triangles work out crossing points
cross_points = np.zeros((len(ix), 2)) # empty array for crossing points
cross_points[:,0] = e_1[ix] + d_ratio * (e_1[ix+1] - e_1[ix]) # x crossings
cross_points[:,1] = t_mrg_e3[ix] + d_ratio * (t_mrg_e3[ix+1] - t_mrg_e3[ix]) # y crossings
print('Intersection',cross_points)
####
ax.scatter(cross_points[:,0], cross_points[:,1], marker = 'o', s=500, c='black')
c = ('darkblue', 'green', 'red', 'purple')
ax.semilogx(e_1, t_mrg_e3, c=c[3])
ax.semilogy()
ax.set_ylim(1e+1,1e+21)
ax.set_xlim(5e-7,5e+1)
plt.plot([ax.get_xlim()[0], ax.get_xlim()[1]],[T_foc, T_foc], c=c[0])
plt.tight_layout()
plt.show()