115 lines
3.4 KiB
Python
115 lines
3.4 KiB
Python
#py -m pip install scipy numpy
|
|
|
|
from scipy.integrate import solve_ivp
|
|
import math
|
|
import numpy as np
|
|
from scipy.constants import pi
|
|
from scipy.interpolate import interp1d
|
|
|
|
|
|
class Ray():
|
|
interpolation = None
|
|
|
|
def __init__(self, x, y, D, Rs, fov, res, accuracy):
|
|
self.x = x
|
|
self.y = y
|
|
self.D = D
|
|
self.Rs = Rs
|
|
self.fovX = fov[0]
|
|
self.fovY = fov[1]
|
|
self.res = res
|
|
self.accuracy = accuracy
|
|
|
|
if Ray.interpolation is None:
|
|
self.interpolate()
|
|
|
|
self.interpolation = Ray.interpolation
|
|
|
|
def interpolate(self):
|
|
maxA = self.fovX
|
|
minA = self.apperentSize(0, 180)
|
|
|
|
linspace = np.linspace(0, 1, self.accuracy)
|
|
angles = minA + (maxA - minA) * linspace**6
|
|
|
|
array = []
|
|
for i in range(self.accuracy):
|
|
sol = self.solver(angles[i])
|
|
deviatedAngle = sol[1][len(sol[1]) - 1]
|
|
array.append([angles[i], np.degrees(deviatedAngle)])
|
|
|
|
#interpolierung der punkte
|
|
points = np.array(array)
|
|
Ray.interpolation = interp1d(points[:, 0], points[:, 1], kind='cubic')
|
|
|
|
def solver(self, alpha):
|
|
if alpha == 0: # skip divided by 0 error
|
|
return [0], [0] # r and phi=0
|
|
|
|
if alpha == 180:
|
|
return [self.D], [0] # if angle= pi then, tan(pi)=0 so 1/tan=1/0
|
|
|
|
# initial value for position and angular speed
|
|
y0 = [1/self.D, 1/(self.D*math.tan(math.radians(alpha)))]
|
|
sol = solve_ivp(fun=self._diff_eq, t_span=[0, 10*pi], y0=y0, method='DOP853', events=[self._eventRs]) #, self._eventR])#,t_eval=np.linspace(0, t_max, 10000)) #dense_output=False
|
|
|
|
if sol.t[-1] == 10*pi:
|
|
raise StopIteration("solver error, alpha reached computation limit (loop number)")
|
|
|
|
phi = np.array(sol.t)
|
|
r = np.abs(1/sol.y[0, :]) # must use this because solver can't be stop before infinity because negative
|
|
|
|
return r, phi
|
|
|
|
def _diff_eq(self, phi, u):
|
|
"""Represent the differential equation : d²u(ɸ)/dɸ²=3/2*Rs*u²(ɸ)-u(ɸ)
|
|
"""
|
|
v0 = u[1] #correspond to u'
|
|
v1 = 3/2*self.Rs*u[0]**2 - u[0] #correspond to u"
|
|
return v0, v1
|
|
|
|
def _eventRs(self, phi, u):
|
|
"""stop solver if radius < black hole radius"""
|
|
with np.errstate(all='ignore'):
|
|
return 1/u[0] - self.Rs
|
|
_eventRs.terminal = True
|
|
|
|
def apperentSize(self, a1, a2):
|
|
|
|
if a2-a1 <= .01: #!!!!!!!!!!!!MaxDeviation später einfügen maybe
|
|
return a2
|
|
|
|
points = np.linspace(a1, a2, 3)
|
|
point = self.solver(points[1])
|
|
|
|
test = point[0][len(point[0]) - 1 ]
|
|
|
|
if test > self.Rs:
|
|
return self.apperentSize(a1, points[1])
|
|
|
|
|
|
elif test <= self.Rs:
|
|
return self.apperentSize(points[1], a2)
|
|
else:
|
|
print("Fehler")
|
|
|
|
def angle(self):
|
|
x = self.x
|
|
y = self.y
|
|
|
|
fovX = math.radians(self.fovX)
|
|
fovY = math.radians(self.fovY)
|
|
|
|
maxX = math.tan(fovX / 2)
|
|
maxY = math.tan(fovY / 2)
|
|
|
|
resX = self.res
|
|
resY = self.res * (maxY / maxX)
|
|
|
|
n = (x / resX) * maxX
|
|
m = (y / resY) * maxY
|
|
|
|
angle = math.degrees(math.atan(math.sqrt(n**2 + m**2)))
|
|
print(angle)
|
|
return self.interpolation(angle)
|