diff --git a/Celestial_Mechanics b/Celestial_Mechanics new file mode 100644 index 0000000..1afe8b4 --- /dev/null +++ b/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() diff --git a/h/h/hw_numpy b/h/h/hw_numpy new file mode 100644 index 0000000..5ad4aa8 --- /dev/null +++ b/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) diff --git a/h/h/numpy_decimal b/h/h/numpy_decimal new file mode 100644 index 0000000..1817db8 --- /dev/null +++ b/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)