Dateien nach "/" hochladen

Ray.py

@@ -0,0 +1,105 @@
+from scipy.integrate import solve_ivp
+import math
+import numpy as np
+from scipy.constants import pi
+from scipy.interpolate import interp1d
+from functools import cache
+import copy
+from PIL import Image
+
+
+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
+
+        n = None
+        m = None
+
+        angle = math.degrees(math.atan2(math.sqrt(n**2 + m**2)))
+        return self.interpolation(angle)