Maxim Blokh

Celestial_Mechanics

@@ -0,0 +1,104 @@
+import numpy as np
+import matplotlib.pyplot as plt
+import math
+import time
+from tqdm import tqdm
+
+G = 6.7e-2
+zerovec = np.array([0, 0, 0], dtype='float')
+
+
+class Star:
+    def __init__(self, mass, vec_p):
+        self.mass = mass
+        self.destroyed = False              
+        self.vec_p = vec_p
+
+
+class CosmicBody:
+    x_cor = []
+    y_cor = []
+    z_cor = []
+
+    def __init__(self, mass, vec_v, vec_p):
+        self.mass = mass
+        self.vec_v = vec_v
+        self.vec_p = vec_p
+        self.destroyed = False
+        self.x_cor.append(self.vec_p[0])
+        self.y_cor.append(self.vec_p[1])
+        self.z_cor.append(self.vec_p[2])
+
+    def gravitate(self, bodys):
+        dv = 0
+        for body in bodys:
+            if body.destroyed == False and body != self:
+                r = - self.vec_p + body.vec_p
+                dv += G * body.mass * dt * r / math.pow(np.linalg.norm(r), 3)
+                self.vec_v += dv
+
+    def destroy(self):
+        self.destroyed = True
+        raise Exception('game over')
+
+    def step(self, dt):
+        self.vec_p += self.vec_v * dt
+        self.x_cor.append(self.vec_p[0])
+        self.y_cor.append(self.vec_p[1])
+        self.z_cor.append(self.vec_p[2])
+
+
+star = Star(10000, zerovec)
+body1 = CosmicBody(1000, np.array([-5., 0., 0.]), np.array([-4., 4., 4.]))
+body2 = CosmicBody(1000, np.array([-4., 0., 0.]), np.array([6., 10., 0.]))
+body3 = CosmicBody(10, np.array([5., 10., 5.]), np.array([5., 0., 0.]))
+bodys = [body1, body2, body3]
+stars = [star]
+
+%matplotlib notebook
+fig = plt.figure(figsize=(5, 5))
+ax = fig.add_subplot(111, projection='3d')
+ax.axes.set_xlim3d(-11, 11)
+ax.axes.set_ylim3d(-5, 10)
+ax.axes.set_zlim3d(-8, 8)
+fig.show()
+fig.canvas.draw()
+
+dt = 0.05
+t0 = 6
+
+for t in tqdm(np.arange(0., t0, dt)):
+    for body in bodys:
+        if body.destroyed == False:
+            body.gravitate(stars)
+            body.gravitate(bodys)
+
+    bodys_destroying = []
+
+    for main_body in bodys:
+        if main_body.destroyed == False:
+            for body in bodys:
+                if body != main_body:
+                    if body.destroyed == False:
+                        if np.linalg.norm(main_body.vec_p - body.vec_p) < 0.5:
+                            bodys_destroying.append(body)
+                            break
+
+    for body in bodys_destroying:
+        body.destroy()
+
+    for body in bodys:
+        if body.destroyed == False:
+            body.step(dt)
+
+    ax.clear()
+    ax.axes.set_xlim3d(-11, 11)
+    ax.axes.set_ylim3d(-5, 10)
+    ax.axes.set_zlim3d(-8, 8)
+    ax.scatter(0, 0, 0, s=100)
+
+    for body in bodys:
+        if body.destroyed == False:
+            ax.scatter(body.x_cor, body.y_cor, body.z_cor,
+                       color='green', marker='o')
+            fig.canvas.draw()

h/h/hw_numpy

@@ -0,0 +1,57 @@
+import matplotlib.pyplot as plt
+import numpy as np
+import scipy.special as sci
+
+
+def poisson_desribution(lmbd, N):
+    if lmbd < 0:
+        raise ValueError("lambda not negative")
+    numpy_arr = np.arange(N)
+    arr = lmbd ** numpy_arr * np.exp(-lmbd) / sci.factorial(numpy_arr)
+    return arr
+
+
+def initial_Moment(arr, k):
+    if not isinstance(k, int):
+        raise ValueError("k must be integer")
+    elif not isinstance(arr, np.ndarray):
+        raise ValueError("Wrong type of array")
+    return np.sum(arr * np.arange(len(arr)) ** k)
+
+
+def Average(arr):
+    return initial_Moment(arr, 1)
+
+
+def Dispersion(arr):
+    return initial_Moment(arr, 2) - initial_Moment(arr, 1) ** 2
+
+
+def Plot(arr):
+    plt.plot(arr)
+    plt.show()
+
+
+def Compare(x, value, error=0.001):
+    if abs(x-value) <= error:
+        return "matches within the error " + str(error)
+    return "not matches within the error " + str(error)
+
+
+if __name__ == '__main__':
+    print("Enter lambda: ")
+    lmbd = 10
+    # also we can use: lmbd = int(input())
+    print("Enter N: ")
+    N = 100
+    # also we can use: N = int(input())
+    print("Enter k: ")
+    k = 2
+    # also we can use: N = int(input())
+    arr = poisson_desribution(lmbd, N)
+    Plot(arr)
+    print("Initial Moment for k =", k, "is: ", initial_Moment(arr, k))
+    print("Average is: ", Average(arr), "and it is",
+          Compare(Average(arr), lmbd), "with real value", lmbd)
+    print("Disrepsion is: ", Dispersion(arr), "and it is",
+          Compare(Dispersion(arr), lmbd), "with real value", lmbd)

h/h/numpy_decimal

@@ -0,0 +1,65 @@
+import matplotlib.pyplot as plt
+import numpy as np
+from decimal import Decimal, getcontext
+import scipy.special as sci
+
+getcontext().prec = 8
+
+
+@np.vectorize
+def factorial_D(x):
+    return Decimal(sci.factorial(int(x)))
+
+
+def poisson_desribution(lmbd, N):
+    if lmbd < 0:
+        raise ValueError("lambda not negative")
+    numpy_arr_D = np.asarray([Decimal(i) for i in range(N)])
+    P = Decimal(lmbd) ** numpy_arr_D * Decimal(-lmbd).exp()
+    return P / factorial_D(numpy_arr_D)
+
+
+def initial_Moment(arr, k):
+    if not isinstance(k, int):
+        raise ValueError("k must be integer")
+    elif not isinstance(arr, np.ndarray):
+        raise ValueError("Wrong type of array")
+    return np.sum(arr * np.arange(len(arr)) ** k)
+
+
+def Average(arr):
+    return initial_Moment(arr, 1)
+
+
+def Dispersion(arr):
+    return initial_Moment(arr, 2) - initial_Moment(arr, 1) ** 2
+
+
+def Plot(arr):
+    plt.plot(arr)
+    plt.show()
+
+
+def Compare(x, value, error=0.001):
+    if abs(x-value) <= error:
+        return "matches within the error " + str(error)
+    return "not matches within the error " + str(error)
+
+
+if __name__ == '__main__':
+    print("Enter lambda: ")
+    lmbd = 50
+    # also we can use: lmbd = int(input())
+    print("Enter N: ")
+    N = 100
+    # also we can use: N = int(input())
+    print("Enter k: ")
+    k = 2
+    # also we can use: k = int(input())
+    arr = poisson_desribution(lmbd, N)
+    Plot(arr)
+    print("Initial Moment for k =", k, "is: ", initial_Moment(arr, k))
+    print("Average is: ", Average(arr), "and it is",
+          Compare(Average(arr), lmbd), "with real value", lmbd)
+    print("Disrepsion is: ", Dispersion(arr), "and it is",
+          Compare(Dispersion(arr), lmbd), "with real value", lmbd)