Галина Прилипко

Galina_Prilipko.py

@@ -0,0 +1,86 @@
+# -*- coding: utf-8 -*-
+"""
+Created on Wed Oct 19 16:41:37 2022
+
+@author: galin
+"""
+
+
+#from scipy.stats import poisson
+from decimal import Decimal
+
+import numpy as np
+import matplotlib.pyplot as plt
+
+
+def poisson_d(N,  mu):
+    #n = np.arange(0,N)
+   # y = poisson.pmf(n, mu) встроенная функция 
+    n1 = np.arange(1,N)
+    n_l=np.log(n1)
+    Y = -mu+n1*np.log(mu)-np.cumsum(n_l) #берем логарифм от пуассона, чтобы было проще считать
+    y = np.exp(Y)
+    y = np.insert(y,0,np.exp(-mu)) #добавляет значение в нуле
+    #y = np.power(mu,n)*np.exp(-mu)/sp.factorial(n) работает на маленьких N
+    return y
+
+
+
+def moments(n,y):
+    try:
+        N=y.size
+
+        x = np.arange(0, N)
+        x=np.power(x,n)
+        return np.dot(x, y)
+    except: print("error: y should be an array")
+
+try:
+    N = float(input('enter N '))
+    mu = float(input('enter average '))
+    assert(mu>0)
+    assert(N>0)
+    assert(N>mu)
+
+except:
+    print("N and mu should be int>0 and N>mu")
+
+x = np.arange(0, N)
+plt.plot(x,poisson_d(N, mu),'bo')
+plt.xlim(0,3*mu)
+plt.suptitle('Float')
+plt.show()
+
+N_dec=Decimal(N)
+mu_dec=Decimal(mu)
+y = poisson_d(N,mu)
+x = np.asarray([Decimal(i) for i in x])
+y = np.asarray([Decimal(i) for i in y])
+
+plt.plot(x,y,'bo')
+plt.xlim(0,3*mu)
+plt.suptitle('Decimal')
+plt.show()
+try:
+   x1 = np.arange(0, N)
+   n=int(input('enter n '))
+   assert(n>0)
+   Y=poisson_d(N, mu)
+   print('moments for floats=',moments(n,Y))
+   a=np.dot(x1, poisson_d(N,mu))
+   print('theoretical average for floats =', a)
+   b = np.power((x1 - a),2)
+   k = np.dot(b,poisson_d(N,mu))
+   print('dispersion for floats =', k)
+
+except: print("n must be an integer>0")
+n_dec=Decimal(n)
+m = moments(n,Y)
+m = Decimal(m)
+print('moments for dec=',m)
+a=np.dot(x, y)
+print('theoretical average for dec =', a)
+b = np.power((x - a),2)
+k = np.dot(b,y)
+
+print('dispersion for dec =', k)

prog-coarse-work.py

@@ -0,0 +1,141 @@
+#!/usr/bin/env python
+# coding: utf-8
+
+# In[19]:
+
+
+import random
+import numpy as np
+import matplotlib.pyplot as plt
+import math
+import itertools
+
+
+class Cell:           
+#класс клетка, в котором хранится количество и характеристики дроплетов
+    bodies = []
+    Dt = 0.1    
+    T=180 
+    def J(t):         
+#метод, определяющий поток, который попал в клетку в момент времени t
+        J0=5 
+        alpha=0.1  
+        return J0*np.sin(t*alpha)
+    def formula(t):     
+#метод, считающий производную обьема каждого дроплета в момент времени t
+        dV_dt = np.zeros(len(Cell.bodies),dtype=float)
+        v=[] 
+        for i in Cell.bodies:
+            v.append(i.volume)
+        v=np.array(v)
+        for j in range(len(Cell.bodies)):
+            dV_dt[j]=Cell.J(t)*np.power(Cell.bodies[j].volume,0.7)/np.sum(np.power(v,0.7))
+        return dV_dt
+    def add_body(body):
+#метод добавления дроплета в клетку. он вызывается, когда мы инициализируем новый дроплет
+        Cell.bodies.append(body)
+    def merge():
+#метод, определяющий слились ли какие-то два дроплета за промежуток времени dt
+        p=np.random.uniform(0,0.0225)   #random probability 
+        m=np.random.randint(len(Cell.bodies), size=2)  #which of them merged
+        if(m[0]!=m[1] and len(Cell.bodies)>1):
+            if random.choices([True,False],weights = [p,1-p])[0]:
+                if Cell.bodies[m[0]].volume>Cell.bodies[m[1]].volume:
+                    Cell.bodies[m[0]].volume=Cell.bodies[m[0]].volume+Cell.bodies[m[1]].volume
+                    Cell.bodies.pop(m[1])
+                else:
+                    Cell.bodies[m[1]].volume=Cell.bodies[m[0]].volume+Cell.bodies[m[1]].volume
+                    Cell.bodies.pop(m[0])
+    def add():
+#метод, определяющий появились ли новые дроплета за промежуток времени dt
+        p=np.random.uniform(0,0.005)   #random probability 
+        if random.choices([True,False],weights = [p,1-p])[0]:
+            Droplets(5)
+    def check():
+#метод, который вызывается каждую итерацию chаnge() чтобы проверить, 
+#что все обьемы положительны и при необходимости удалить отрицательные
+        ind=[]
+        for i, v in enumerate(Cell.bodies):
+            if v.volume<=0: Cell.bodies.pop(i)
+
+    def change():
+#метод, считающий, как изменились обьемы дроплетов за dt
+        n_steps = int(round((Cell.T)/Cell.Dt))
+        t_arr = np.zeros(n_steps + 1)  
+        for i in range (1, n_steps + 1):
+            t = t_arr[i-1]
+            dVdt = Cell.formula(t) 
+            V0=np.zeros(len(Cell.bodies)+1)
+            for idx,val in np.ndenumerate(Cell.bodies):
+                V0[idx]=val.volume
+            for idx,val in np.ndenumerate(Cell.bodies):
+                val.volume=V0[idx] +Cell.Dt*dVdt[idx] 
+            t_arr[i] = t + Cell.Dt
+            Cell.check()
+            Cell.merge()
+            Cell.add()
+            T=np.zeros(len(Cell.bodies))+t  
+            v=[] 
+            col=[]
+            for i in Cell.bodies:
+                v.append(i.volume)
+                col.append(i.colour)
+            v=np.array(v)
+            plt.xlabel('t (in miliseconds)', fontsize = 12)
+            plt.ylabel('V(t)', fontsize = 12)
+            plt.scatter(T, v, s=5 ,c=col)
+
+
+class Droplets:
+
+    colours = itertools.cycle(['#3366CC','#DD4477','#009292','#FF6DB6','#DC3912','#FF9900','#109618','#990099',
+                               '#0099C6','#004949','#7E0021','#314004','#7375B5','#4B1F6F','#6B452B','#A89985','#370335',
+                               '#A50E82'])
+    def __init__(
+        self,
+        volume
+    ):
+        self.volume = volume
+        self.colour = next(Droplets.colours)
+        Cell.add_body(self)
+
+#инициализируем 7 дроплетов. их число ограниченно размерами клетки 
+#(то есть суммарный обьем не должен быть больше клетки, подробнее - в тз)
+Droplets(20)
+Droplets(20)
+Droplets(20)
+Droplets(20)
+Droplets(20)
+Droplets(20)
+Droplets(20)
+Droplets(20)
+Droplets(20) 
+Cell.change()
+
+
+# In[20]:
+
+
+import unittest
+
+class TestCell(unittest.TestCase):
+  def setUp(self):
+    self.cell = Cell()
+
+  def test_J(self):   #тестируем метод, выдающий обьем
+    self.assertEqual(Cell.J(7), 3.2210884361884555)
+
+unittest.main(argv=[''], verbosity=2, exit=False)
+
+
+# In[ ]:
+
+
+
+
+
+# In[ ]:
+
+
+
+

solar.py

@@ -0,0 +1,148 @@
+# -*- coding: utf-8 -*-
+"""
+Created on Fri Dec 23 13:52:40 2022
+
+@author: galin
+"""
+
+import numpy as np
+import matplotlib.pyplot as plt
+import math
+import itertools
+from matplotlib import animation
+from tqdm import tqdm
+%matplotlib notebook
+
+class SolarSystem:
+    def __init__(self, size):
+        self.size = size
+        self.bodies = []
+
+        self.fig, self.ax = plt.subplots(
+            1,
+            1,
+            subplot_kw={"projection": "3d"},
+            figsize=(self.size / 50, self.size / 50),
+        )
+        self.fig.tight_layout()
+        self.ax.view_init(1, 1)
+    def add_body(self, body):
+        self.bodies.append(body)
+
+    def update_all(self):
+        for body in self.bodies:
+            body.move()
+            body.draw()
+    def draw_all(self):
+        self.ax.set_xlim3d(-self.size / 2, self.size / 2)
+        self.ax.set_ylim3d(-self.size / 2, self.size / 2)
+        self.ax.set_zlim3d(-self.size / 2, self.size / 2)
+        plt.pause(0.01)
+        self.ax.clear()
+
+
+    def calculate_all_body_interactions(self):
+        bodies_copy = self.bodies.copy()
+        for idx, first in enumerate(bodies_copy):
+            for second in bodies_copy[idx + 1:]:
+                first.accelerate_due_to_gravity(second)
+
+
+
+class SpaceBody:
+    min_display_size = 10
+    display_log_base = 1.3
+    def __init__(
+        self,
+        solar_system,
+        mass,
+        position=np.array([0, 0, 0]),
+        velocity=np.array([0, 0, 0]),
+    ):
+        self.solar_system = solar_system
+        self.mass = mass
+        self.position = position
+        self.velocity = velocity
+        self.display_size = max(
+            math.log(self.mass, self.display_log_base),
+            self.min_display_size,
+        )
+        self.colour = "black"
+        self.solar_system.add_body(self)
+    def move(self):
+        self.position = (
+            self.position[0] + self.velocity[0],
+            self.position[1] + self.velocity[1],
+            self.position[2] + self.velocity[2],
+        )
+    def draw(self):
+        self.solar_system.ax.plot(
+            *self.position,
+            marker="o",
+            markersize=self.display_size,
+            color=self.colour
+        )
+
+    def accelerate_due_to_gravity(self, other):
+        distance = np.array(other.position) - np.array(self.position)
+        distance_mag = np.linalg.norm(distance)
+        force_mag = self.mass * other.mass / (distance_mag ** 2)
+        force = distance * force_mag/distance_mag
+        reverse = 1.
+        for body in self, other:
+            acceleration = force / body.mass
+            body.velocity =body.velocity + acceleration * reverse
+            reverse = -1.
+class Sun(SpaceBody):
+    def __init__(
+        self,
+        solar_system,
+        mass=10000,
+        position=np.array([0, 0, 0]),
+        velocity=np.array([0, 0, 0]),
+    ):
+        super(Sun, self).__init__(solar_system, mass, position, velocity)
+        self.colour = "yellow"
+class Planet(SpaceBody):
+    colours = itertools.cycle([(1, 0, 0), (0, 1, 0), (0, 0, 1)])
+    def __init__(
+        self,
+        solar_system,
+        mass=10,
+        position=np.array([0, 0, 0]),
+        velocity=np.array([0, 0, 0]),
+    ):
+        super(Planet, self).__init__(solar_system, mass, position, velocity)
+        self.colour = next(Planet.colours)
+
+
+
+solar_system = SolarSystem(400)
+sun = Sun(solar_system)
+planets = (
+    Planet(
+        solar_system,
+        position=(150, 50, 0),
+        velocity=(0, 5, 5),
+    ),
+    Planet(
+        solar_system,
+        mass=20,
+        position=(100, -50, 150),
+        velocity=(5, 0, 0)
+    )
+)
+t0=4
+dt=0.01
+
+for t in tqdm(np.arange(0.,t0,dt)):
+    angle = 60 + 60 * t / t0
+    solar_system.draw_all()
+    solar_system.calculate_all_body_interactions()
+    solar_system.update_all()
+    solar_system.ax.set_xlim3d(-200,200)
+    solar_system.ax.set_ylim3d(-200,200)
+    solar_system.ax.set_zlim3d(-200,200)
+    solar_system.fig.canvas.draw()
+    solar_system.ax.view_init(30 - angle * 0.2, angle)
+