diff --git a/Poisson.py b/Poisson.py new file mode 100644 index 0000000..79f1148 --- /dev/null +++ b/Poisson.py @@ -0,0 +1,55 @@ +import numpy as np +import scipy.special +import matplotlib.pyplot as plt + +def poisson(aver, N): + n = np.arange(0, N) + p = (np.power(aver, n)*np.exp(-aver))/scipy.special.factorial(n) + if aver >= 0: + return p + else: + return 0 + +def poisson_1(aver, N): + p = (np.power(aver, N) * np.exp(-aver)) / scipy.special.factorial(N) + if aver >= 0: + return p + else: + return 0 + +def initial_moment(n, k): + N = np.size(n) + arr = np.arange(N) + x = np.power(arr, k) + y = n + m = np.dot(y, x) + if isinstance(k, int) == False or not(np.size(np.array(n)) > 1): + return 0 + else: + return m + +def average(n): + m = initial_moment(n, 1) + return m + +def dispersion(n): + m = initial_moment(n, 2) - initial_moment(n, 1) ** 2 + return m + +def test(aver, N): + a = poisson(aver, N) + b = average(a) + c = dispersion(a) + plt.plot(a) + plt.scatter(b, poisson_1(b, b)) + plt.scatter(c, poisson_1(c, c)) + plt.show() + +#test +a = poisson(3, 100) +b = average(a) +c = dispersion(a) +print('a = ', 3) +print('b = ', b) +print('c = ', c) +test(3, 100) \ No newline at end of file diff --git a/Poisson_decimal.py b/Poisson_decimal.py new file mode 100644 index 0000000..756ad85 --- /dev/null +++ b/Poisson_decimal.py @@ -0,0 +1,57 @@ +import numpy as np +import matplotlib.pyplot as plt +from decimal import Decimal + + +def exp(x): + return x.exp() + +@np.vectorize +def fact(n): + f = Decimal('1') + while n > Decimal('1'): + f = f * n + n = n - Decimal('1') + return f + +def poisson(aver, N): + n = np.arange(N+1) + p = ((aver**n)*exp(-aver))/fact(n) + if aver >= 0: + return p + else: + return 0 + +def initial_moment(n, k): + arr = np.arange(np.size(n)) + m = n*arr*k + if isinstance(k, int) == False or not(np.size(np.array(n)) > 1): + return 0 + else: + return m.sum() + +def average(n): + m = initial_moment(n, 1) + return m + +def dispersion(n): + m = initial_moment(n, 2) - initial_moment(n, 1)**2 + return m + +def test(aver, N): + a = poisson(aver, N) + b = average(a) + c = dispersion(a) + plt.plot(a) + b1 = poisson(b, b) + plt.scatter(b, b1[np.size(b1)-1]) + plt.show() + +#test +a = poisson(Decimal('3'), Decimal('100')) +b = average(a) +c = dispersion(a) +print('a = ', Decimal('3')) +print('b = ', b) +print('c = ', c) +test(Decimal('3'), Decimal('100')) \ No newline at end of file diff --git a/Solar_system3d.py b/Solar_system3d.py new file mode 100644 index 0000000..2db5621 --- /dev/null +++ b/Solar_system3d.py @@ -0,0 +1,128 @@ +import numpy as np +import matplotlib.pyplot as plt +from matplotlib import animation + + +class Star: + def __init__(self, mass: float): + self.mass = mass + self.vec_P = np.array([[0.0, 0.0, 0.0],]) + + +class CosmicBody: + G = 6.67430e-11 + sec_in_day = 86400 + + def __init__(self, mass: float, vec_v: np.ndarray, vec_P: np.ndarray): + self.mass = mass + self.vec_v = vec_v + self.vec_P = vec_P + + def gravitate(self, bodys: list, t0: float): + v = self.vec_v[-1] + v_delt = np.array([[0.0, 0.0, 0.0]]) + + for affect_bod in bodys: + r_vec = self.vec_P[-1] - affect_bod.vec_P[-1] + r = sum(r_vec**2)**0.5 + + v_delt += (-self.G * self.mass * + affect_bod.mass * + r_vec / r**3) * self.sec_in_day * t0 / self.mass + + r_new = self.vec_P[-1] + (v + v_delt) * self.sec_in_day * t0 + self.vec_P = np.append(self.vec_P, r_new, axis=0) + self.vec_v += v_delt + + def destroy(self, bodys): + for affect_bod in bodys: + if sum(((self.vec_P[-1] - affect_bod.vec_P[-1])**2))**0.5 < 1.75e9: + return True, affect_bod + + +def random_vec(): + neg_or_pos = [1 if np.random.random() > 0.5 + else -1 for _ in range(3)] + + return np.random.rand(1, 3) * neg_or_pos + +Sun = Star(2.5e31) + +first_body = CosmicBody(np.random.random() * 1.0e29, + random_vec()*1.0e4, + random_vec()*5.0e11) + +all_bodyes = [Sun, first_body] +bodyes_for_plot = [{first_body: 1}] + +t0 = 0.0 +delta_t = 1.0 +t = 365 + +while t0 < t and len(all_bodyes) > 1: + for i, body in enumerate(all_bodyes[1:], start=1): + affect_bodyes = all_bodyes[:i] + all_bodyes[i+1:] + body.gravitate(affect_bodyes, delta_t) + destroy = body.destroy(affect_bodyes) + if destroy: + all_bodyes.remove(body) + if destroy[1] is not Sun: + all_bodyes.remove(destroy[1]) + bodyes_for_plot += [{body: body.vec_P.shape[0] for body in all_bodyes[1:]}] + if np.random.random() > 0.9: + new_body = CosmicBody(np.random.random() * 1.0e29, + random_vec()*5.0e4, + random_vec()*5.0e11) + all_bodyes.append(new_body) + + + t0 += delta_t + + + +fig = plt.figure(figsize=(8, 8)) +ax = fig.add_subplot(projection="3d") + + +ploted_bodyes = {} +def animate(i): + global ploted_bodyes + for body in bodyes_for_plot[i]: + index = bodyes_for_plot[i][body] + coord = body.vec_P + if body not in ploted_bodyes: + ploted_bodyes[body] = (ax.plot(coord[index-1][0], coord[index-1][1], coord[index-1][2], + marker='o', + markersize=6, + markeredgecolor="black", + markerfacecolor="black")[0], + ax.plot([], [], [], '-g', + lw=0.2, + c="blue")[0]) + else: + ploted_bodyes[body][0].set_data_3d(coord[index-1][0], coord[index-1][1], coord[index-1][2]) + ploted_bodyes[body][1].set_data_3d(coord[:index, 0], coord[:index, 1], coord[:index, 2]) + + for remove_el in set(ploted_bodyes) - set(bodyes_for_plot[i]): + remove_plot = ploted_bodyes.pop(remove_el) + remove_plot[0].remove() + remove_plot[1].remove() + + sun_point = ax.plot([0], [0], [0], + marker='o', + markersize=10, + markeredgecolor="black", + markerfacecolor="yellow")[0] + + ax.set(xlim3d=(-1.2*5.0e11, 1.2*5.0e11), xlabel='X') + ax.set(ylim3d=(-1.2*5.0e11, 1.2*5.0e11), ylabel='Y') + ax.set(zlim3d=(-1.2*5.0e11, 1.2*5.0e11), zlabel='Z') + return (sun_point,) + tuple(ploted_bodyes.values()) + + +solar_system_anim = animation.FuncAnimation(fig, + func=animate, + frames = len(bodyes_for_plot), + interval=1) + +plt.show() \ No newline at end of file