Как найти точку пересечения прямой линии и подобного кривой набора точек? - PullRequest
0 голосов
/ 28 апреля 2020

Как найти точку пересечения прямой линии и кривого набора точек, пожалуйста? Я обработал в соответствии с этот вопрос , однако результат неверен. Где ошибка, пожалуйста? Важная часть обозначена #### в коде. Остальное - для генерации списка для конструирования кода и построения результата.

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()
Добро пожаловать на сайт PullRequest, где вы можете задавать вопросы и получать ответы от других членов сообщества.
...