From 95c359a47ef580a201724c559b8e042506337cd1 Mon Sep 17 00:00:00 2001
From: jakut77 <jakut77@zedat.fu-berlin.de>
Date: Tue, 6 Feb 2024 12:54:11 +0100
Subject: [PATCH] updated to Octree
---
tasks/bht/2D/bht.py | 316 +++++++++++++++
tasks/bht/2D/galaxysimulation.py | 147 +++++++
tasks/bht/2D/simulate.py | 5 +
tasks/bht/__pycache__/bht_oct.cpython-311.pyc | Bin 0 -> 19590 bytes
.../galaxysimulation.cpython-311.pyc | Bin 10722 -> 12599 bytes
tasks/bht/bht_oct.py | 362 ++++++++++++++++++
tasks/bht/galaxysimulation.py | 103 +++--
tasks/bht/print_tree.py | 2 +-
tasks/bht/simulate.py | 6 +-
tasks/bht/system.py | 57 +++
10 files changed, 957 insertions(+), 41 deletions(-)
create mode 100644 tasks/bht/2D/bht.py
create mode 100644 tasks/bht/2D/galaxysimulation.py
create mode 100644 tasks/bht/2D/simulate.py
create mode 100644 tasks/bht/__pycache__/bht_oct.cpython-311.pyc
create mode 100644 tasks/bht/bht_oct.py
create mode 100644 tasks/bht/system.py
diff --git a/tasks/bht/2D/bht.py b/tasks/bht/2D/bht.py
new file mode 100644
index 0000000..30f345a
--- /dev/null
+++ b/tasks/bht/2D/bht.py
@@ -0,0 +1,316 @@
+import numpy as np
+
+
+class MainApp:
+
+ def __init__(self):
+ """Initialize the MainApp with a root node."""
+ self.rootNode = TreeNode(x=0, y=0, width=100, height=100)
+
+ def BuildTree(self, particles):
+ """Build the Quadtree by inserting particles.
+
+ Args:
+ particles (list): A list of Particle objects to be inserted into the Quadtree.
+ """
+ self.ResetTree() # Empty the tree
+ for particle in particles:
+ self.rootNode.insert(particle)
+
+ def ResetTree(self):
+ """Reset the Quadtree by reinitializing the root node."""
+ self.rootNode = TreeNode(x=0, y=0, width=100, height=100)
+
+
+class Particle:
+
+ def __init__(self, x, y, mass):
+ """Initialize a Particle with x and y coordinates and mass."""
+ self.x = x
+ self.y = y
+ self.mass = mass
+ self.vx = 0.0 # Velocity component in x direction
+ self.vy = 0.0 # Velocity component in y direction
+ self.fx = 0.0 # Force component in x direction
+ self.fy = 0.0 # Force component in y direction
+
+
+class TreeNode:
+
+ def __init__(self, x, y, width, height):
+ """Initialize a TreeNode representing a quadrant in the Quadtree.
+
+ Args:
+ x (float): x-coordinate of the node.
+ y (float): y-coordinate of the node.
+ width (float): Width of the node.
+ height (float): Height of the node.
+ """
+ self.x = x
+ self.y = y
+ self.width = width
+ self.height = height
+ self.particle = None # Particle contained in this node
+ self.center_of_mass = np.array([(x + width) / 2, (y + width) / 2])
+ self.total_mass = 0
+ self.children = np.empty(4, dtype=object) # Children nodes (SW, SE, NW, NE)
+
+ def contains(self, particle):
+ """Check if the particle is within the bounds of this node.
+
+ Args:
+ particle (Particle): The particle to be checked.
+
+ Returns:
+ bool: True if the particle is within the node's bounds, False otherwise.
+ """
+ return (self.x <= particle.x < self.x + self.width and self.y <= particle.y < self.y + self.height)
+
+ def insert(self, particle):
+ """Insert a particle into the Quadtree.
+
+ Args:
+ particle (Particle): The particle to be inserted.
+
+ Returns:
+ bool: True if the particle is inserted, False otherwise.
+ """
+ if not self.contains(particle):
+ return False # Particle doesn't belong in this node
+
+ if self.particle is None and all(child is None for child in self.children):
+ # If the node is empty and has no children, insert particle here.
+ self.particle = particle
+ #print(f'particle inserted: x={round(self.particle.x,2)}, y={round(self.particle.y,2)}')
+ return True # Particle inserted in an empty node
+
+ if all(child is None for child in self.children):
+ # If no children exist, create and insert both particles
+ self.subdivide()
+ self.insert(self.particle) # Reinsert existing particle
+ self.insert(particle) # Insert new particle
+ self.particle = None # Clear particle from this node
+ else:
+ # If the node has children, insert particle in the child node.
+ quad_index = self.get_quadrant(particle)
+ if self.children[quad_index] is None:
+ # Create a child node if it doesn't exist
+ self.children[quad_index] = TreeNode(self.x + (quad_index % 2) * (self.width / 2),
+ self.y + (quad_index // 2) * (self.height / 2), self.width / 2,
+ self.height / 2)
+ self.children[quad_index].insert(particle)
+
+ def subdivide(self):
+ """Subdivide the node into four quadrants."""
+ sub_width = self.width / 2
+ sub_height = self.height / 2
+ self.children[0] = TreeNode(self.x, self.y, sub_width, sub_height) # SW
+ self.children[1] = TreeNode(self.x + sub_width, self.y, sub_width, sub_height) # SE
+ self.children[2] = TreeNode(self.x, self.y + sub_height, sub_width, sub_height) # NW
+ self.children[3] = TreeNode(self.x + sub_width, self.y + sub_height, sub_width, sub_height) # NE
+
+ def get_quadrant(self, particle):
+ """Determine the quadrant index for a particle based on its position.
+
+ Args:
+ particle (Particle): The particle to determine the quadrant index for.
+
+ Returns:
+ int: Quadrant index (0 for SW, 1 for SE, 2 for NW, 3 for NE).
+ """
+ mid_x = self.x + self.width / 2
+ mid_y = self.y + self.height / 2
+ quad_index = (particle.x >= mid_x) + 2 * (particle.y >= mid_y)
+ return quad_index
+
+ def print_tree(self, depth=0):
+ """Print the structure of the Quadtree.
+
+ Args:
+ depth (int): Current depth in the tree (for indentation in print).
+ """
+ if self.particle:
+ print(
+ f"{' ' * depth}Particle at ({round(self.particle.x,2)}, {round(self.particle.y,2)}) in Node ({self.x}, {self.y}), size={self.width}"
+ )
+ else:
+ print(f"{' ' * depth}Node ({self.x}, {self.y}) - Width: {self.width}, Height: {self.height}")
+ for child in self.children:
+ if child:
+ child.print_tree(depth + 2)
+
+ def ComputeMassDistribution(self):
+ """Compute the mass distribution for the tree nodes.
+
+ This function calculates the total mass and the center of mass
+ for each node in the Quadtree. It's a recursive function that
+ computes the mass distribution starting from the current node.
+
+ Note:
+ This method modifies the 'mass' and 'center_of_mass' attributes
+ for each node in the Quadtree.
+
+ Returns:
+ None
+ """
+ if self.particle is not None:
+ # Node contains only one particle
+ self.center_of_mass = np.array([self.particle.x, self.particle.y])
+ self.total_mass = self.particle.mass
+ else:
+ # Multiple particles in node
+ total_mass = 0
+ center_of_mass_accumulator = np.array([0.0, 0.0])
+
+ for child in self.children:
+ if child is not None:
+ # Recursively compute mass distribution for child nodes
+ child.ComputeMassDistribution()
+ total_mass += child.total_mass
+ center_of_mass_accumulator += child.total_mass * child.center_of_mass
+
+ if total_mass > 0:
+ self.center_of_mass = center_of_mass_accumulator / total_mass
+ self.total_mass = total_mass
+ else:
+ # If total mass is 0 or no child nodes have mass, leave values as default
+ pass
+ #self.center_of_mass = np.array([(x+width)/2, (y+width)/2])
+ #self.total_mass = 0
+
+ def CalculateForceFromTree(self, target_particle, theta=1.0):
+ """Calculate the force on a target particle using the Barnes-Hut algorithm.
+
+ Args:
+ target_particle (Particle): The particle for which the force is calculated.
+ theta (float): The Barnes-Hut criterion for force approximation.
+
+ Returns:
+ np.ndarray: The total force acting on the target particle.
+ """
+
+ total_force = np.array([0.0, 0.0])
+
+ if self.particle is not None:
+ # Node contains only one particle
+ if self.particle != target_particle:
+ # Calculate gravitational force between target_particle and node's particle
+ force = self.GravitationalForce(target_particle, self.particle)
+ total_force += force
+ else:
+ if self.total_mass == 0:
+ return total_force
+
+ r = np.linalg.norm(np.array([target_particle.x, target_particle.y]) - self.center_of_mass)
+ d = max(self.width, self.height)
+
+ if d / r < theta:
+ # Calculate gravitational force between target_particle and "node particle" representing cluster
+ node_particle = Particle(self.center_of_mass[0], self.center_of_mass[1], self.total_mass)
+ force = self.GravitationalForce(target_particle, node_particle)
+ total_force += force
+ else:
+ for child in self.children:
+ if child is not None:
+ # Recursively calculate force from child nodes
+ if target_particle is not None: # Check if the target_particle is not None
+ force = child.CalculateForceFromTree(target_particle)
+ total_force += force
+
+ return total_force
+
+ def CalculateForce(self, target_particle, particle, theta=1.0):
+ """Calculate the gravitational force between two particles.
+
+ Args:
+ target_particle (Particle): The particle for which the force is calculated.
+ particle (Particle): The particle exerting the force.
+
+ Returns:
+ np.ndarray: The force vector acting on the target particle due to 'particle'.
+ """
+ force = np.array([0.0, 0.0])
+ print('function CalculateForce is called')
+ if self.particle is not None:
+ # Node contains only one particle
+ if self.particle != target_particle:
+ # Calculate gravitational force between target_particle and node's particle
+ force = self.GravitationalForce(target_particle, self.particle)
+ else:
+ if target_particle is not None and particle is not None: # Check if both particles are not None
+ r = np.linalg.norm(
+ np.array([target_particle.x, target_particle.y]) - np.array([particle.x, particle.y]))
+ d = max(self.width, self.height)
+
+ if d / r < theta:
+ # Calculate gravitational force between target_particle and particle
+ force = self.GravitationalForce(target_particle, particle)
+ else:
+ for child in self.children:
+ if child is not None:
+ # Recursively calculate force from child nodes
+ force += child.CalculateForce(target_particle, particle)
+ return force
+
+ def GravitationalForce(self, particle1, particle2):
+ """Calculate the gravitational force between two particles.
+
+ Args:
+ particle1 (Particle): First particle.
+ particle2 (Particle): Second particle.
+
+ Returns:
+ np.ndarray: The gravitational force vector between particle1 and particle2.
+ """
+ #G = 6.674 * (10 ** -11) # Gravitational constant
+ G = 1
+
+ dx = particle2.x - particle1.x
+ dy = particle2.y - particle1.y
+ cutoff_radius = 5
+ r = max(np.sqrt(dx**2 + dy**2), cutoff_radius)
+
+ force_magnitude = G * particle1.mass * particle2.mass / (r**2)
+ force_x = force_magnitude * (dx / r)
+ force_y = force_magnitude * (dy / r)
+
+ return np.array([force_x, force_y])
+
+ # Helper method to retrieve all particles in the subtree
+ def particles_in_subtree(self):
+ """Retrieve all particles in the subtree rooted at this node.
+
+ Returns:
+ list: A list of particles in the subtree rooted at this node.
+ """
+ particles = []
+ if self.particle is not None:
+ particles.append(self.particle)
+ else:
+ for child in self.children:
+ if child is not None:
+ particles.extend(child.particles_in_subtree())
+ return particles
+
+ def compute_center_of_mass(self):
+ """Compute the center of mass for the node."""
+ print('Function compute_center_of_mass is called')
+ if self.particle is not None:
+ self.center_of_mass = np.array([self.particle.x, self.particle.y])
+ self.mass = self.particle.mass
+ else:
+ total_mass = 0
+ center_of_mass_accumulator = np.array([0.0, 0.0])
+
+ for child in self.children:
+ if child is not None:
+ child.compute_center_of_mass()
+ total_mass += child.mass
+ center_of_mass_accumulator += child.mass * child.center_of_mass
+
+ if total_mass > 0:
+ self.center_of_mass = center_of_mass_accumulator / total_mass
+ self.mass = total_mass
+ else:
+ self.center_of_mass = np.array([0.0, 0.0])
+ self.mass = 0
diff --git a/tasks/bht/2D/galaxysimulation.py b/tasks/bht/2D/galaxysimulation.py
new file mode 100644
index 0000000..f8aae98
--- /dev/null
+++ b/tasks/bht/2D/galaxysimulation.py
@@ -0,0 +1,147 @@
+import matplotlib.pyplot as plt
+import numpy as np
+from bht import MainApp, Particle
+
+
+# Simulation class for galaxy collision using Barnes-Hut and Euler method
+class GalaxyCollisionSimulation:
+
+ def __init__(self, num_particles):
+ # Initialize the simulation with a given number of particles
+ self.num_particles = num_particles
+ self.particles = self.initialize_particles(initial='random') # Generate particles
+ self.barnes_hut = MainApp() # Initialize Barnes-Hut algorithm instance
+ self.barnes_hut.BuildTree(self.particles) # Build the Barnes-Hut tree with particles
+ self.barnes_hut.rootNode.ComputeMassDistribution() #Compute the center of mass of the tree nodes
+ self.particle_positions = {
+ i: [(particle.x, particle.y)]
+ for i, particle in enumerate(self.particles)
+ } # Store particle positions for each time step
+ self.particle_velocities = {
+ i: [(particle.vx, particle.vy)]
+ for i, particle in enumerate(self.particles)
+ } # Store particle velocities for each time step
+ self.particle_forces = {
+ i: [(particle.fx, particle.fy)]
+ for i, particle in enumerate(self.particles)
+ } # Store particle velocities for each time step
+ self.system_center_of_mass = {
+ 0: (self.barnes_hut.rootNode.center_of_mass[0], self.barnes_hut.rootNode.center_of_mass[1])
+ }
+
+ def initialize_particles(self, initial='random'):
+
+ if initial == 'two':
+ particles = [Particle(80, 50, 1000), Particle(20, 50, 1000)]
+ #Particle(0.5, 0, 10000),
+ #Particle(-0.5, 0, 10000)]
+ particles[0].vx, particles[0].vy = 0, 3
+ particles[1].vx, particles[1].vy = 0, -3
+ #particles=particles[:2]
+ return particles
+
+ elif initial == 'random':
+ # Generate random particles within a specified range
+ particles = []
+ for _ in range(self.num_particles):
+ x = np.random.uniform(0, 100)
+ y = np.random.uniform(0, 100)
+ #vx = 0
+ #vy = 0
+ #mass = np.random.uniform(10, 100)
+ mass = 10000
+ particles.append(Particle(x, y, mass))
+ return particles
+
+ def simulate(self, num_steps, time_step):
+ # Simulate the galaxy collision for a certain number of steps
+ for step in range(num_steps):
+ # For each step, build the tree
+ self.barnes_hut.BuildTree(self.particles)
+ self.barnes_hut.rootNode.ComputeMassDistribution()
+ if step in np.arange(0, num_steps, 1000):
+ print(f'===Quadtree at step {step}===')
+ self.barnes_hut.rootNode.print_tree()
+ self.evaluate_forces() # Evaluate forces acting on each particle
+ self.integrate(time_step, 'velocity_verlet') # Integrate to update particle positions
+
+ self.system_center_of_mass[step] = (self.barnes_hut.rootNode.center_of_mass[0],
+ self.barnes_hut.rootNode.center_of_mass[1])
+ # Store particle positions, velocities and forces for each time step
+ for i, particle in enumerate(self.particles):
+ self.particle_positions[i].append((particle.x, particle.y))
+ self.particle_velocities[i].append((particle.vx, particle.vy))
+ self.particle_forces[i].append((particle.fx, particle.fy))
+
+ def evaluate_forces(self):
+ # Evaluate forces acting on each particle using Barnes-Hut algorithm
+ for particle in self.particles:
+ force = self.barnes_hut.rootNode.CalculateForceFromTree(particle)
+ particle.fx, particle.fy = force
+
+ def calculate_force(self, particle):
+ # Calculate force acting on a single particle
+ force = self.barnes_hut.rootNode.CalculateForceFromTree(particle)
+ particle.fx, particle.fy = force
+
+ def integrate(self, time_step, method):
+ # Integrate to update particle positions based on calculated forces
+ if method == 'explicit_euler':
+ for particle in self.particles:
+ particle.vx += particle.fx * time_step / particle.mass # Update velocity components
+ particle.vy += particle.fy * time_step / particle.mass
+ particle.x += particle.vx * time_step # Update particle positions using velocity
+ particle.y += particle.vy * time_step
+
+ elif method == 'velocity_verlet':
+ for particle in self.particles:
+ particle.vx += particle.fx * time_step / (2 * particle.mass) # First half-step
+ particle.vy += particle.fy * time_step / (2 * particle.mass)
+ particle.x += particle.vx * time_step # Full step for position
+ particle.y += particle.vy * time_step
+
+ self.calculate_force(particle) # Recalculate force for the particle
+
+ particle.vx += particle.fx * time_step / (2 * particle.mass) # Second half-step
+ particle.vy += particle.fy * time_step / (2 * particle.mass)
+
+ def print_particle_positions(self):
+ # Print the positions of all particles for each time step
+ for i, positions in self.particle_positions.items():
+ print(f"Particle {i + 1} Positions:")
+ for step, pos in enumerate(positions):
+ print(f" Time Step {step}: Position ({pos[0]}, {pos[1]})")
+
+ def display_snapshots(self, num_shots, fix_axes=True):
+ # Display pictures for each time step
+ num_time_steps = len(next(iter(self.particle_positions.values())))
+ print('===Quadtree at the end===')
+ self.barnes_hut.rootNode.print_tree() # print the tree
+ axes_type = "fixed" if fix_axes else "unfixed"
+ #Fixed axis limits
+ print(f"For {axes_type} axis limits, num_particles={self.num_particles}")
+ for step in np.arange(0, num_time_steps, int(num_time_steps / num_shots)):
+ fig, ax = plt.subplots()
+ if fix_axes:
+ ax.set_xlim(0, 100)
+ ax.set_ylim(0, 100)
+ ax.set_xlabel('X-axis')
+ ax.set_ylabel('Y-axis')
+ ax.set_title(f'Time Step {step}')
+ ax.grid()
+
+ print(f'Step {step}')
+ for particle_index, positions in self.particle_positions.items():
+ x, y = positions[step]
+ vx, vy = self.particle_velocities[particle_index][step]
+ fx, fy = self.particle_forces[particle_index][step]
+ #mass = self.particles[particle_index].mass
+ ax.scatter(x, y, c='blue', s=20, alpha=0.5) # Plot particle position for the current time step
+ c_x, c_y = self.system_center_of_mass[step]
+ ax.scatter(c_x, c_y, c='orange', s=20)
+ print(
+ f'i={particle_index}: x={round(x,2)}, y={round(y,2)}, vx={round(vx,2)}, vy={round(vy,2)}, fx={round(fx,2)}, fy={round(fy,2)}'
+ )
+ #print(y)
+
+ plt.show()
diff --git a/tasks/bht/2D/simulate.py b/tasks/bht/2D/simulate.py
new file mode 100644
index 0000000..c703e42
--- /dev/null
+++ b/tasks/bht/2D/simulate.py
@@ -0,0 +1,5 @@
+from galaxysimulation import GalaxyCollisionSimulation
+
+sim = GalaxyCollisionSimulation(num_particles=6)
+sim.simulate(num_steps=7000, time_step=0.001)
+sim.display_snapshots(200, fix_axes=True)
diff --git a/tasks/bht/__pycache__/bht_oct.cpython-311.pyc b/tasks/bht/__pycache__/bht_oct.cpython-311.pyc
new file mode 100644
index 0000000000000000000000000000000000000000..040eedbf568386874364ac28f4d803d8393a67f5
GIT binary patch
literal 19590
zcmZ3^%ge>Uz`)=xdnh%>mVx0hhy%kcP{!vN1_p-d3@HpLj5!QZj7$vf3@J=43@OaZ
zm>3vVGeK1`L@|M7SXvlTSkYvd!7^+u3@PkrGAt>K!3>%lFF}U*C4)#9W`r_6D}ZfJ
zWr$)-VTfW%VQgWDVg}hC#gf9@!Vtxp!V=7&$$E?3H!(BMv7jIs?jQyR1`Y-W24)5Z
zhR;@DLrP$VGB7ZtFoGdd3R4Qx8s=qe3=FH;pjI-}FvP=ELX5A{@XX6B$xO`2tV&fV
z$w*az+NMyRS(2fUs8E!jU!stgpOUIq1hW6-2?hp+m29^dD{e7X-eRn}#af=3Qj&3t
zEh9BEJ)`6nYf5TCNrs;$(=CpWqSREs{FKyN93V45EKR0cEXAogX}35F5{pVQlXFsw
zi$R7eDEzY04=qkDD%Q^~O-j@+&o9a@E=WvH)lVx*%}mu#E-om_&q_@$(Jx6X&MwwZ
z$|wQh`26G&y@JYH9P#m)d6^~g@l^s)$LqmFKprnvVqjosV0bE_u)uHy<BGry$p@5A
z1YeK{yeJWPMIx|+rH8kJw}baHD2!0T7?e!FVVnmJ;~ItqASD>MhH)7q1H)=~xH2--
zFuE|r%G5H|Fl50+7l1s9sh0^|Z;vcE;V@Orbt=uwNdZTRe{xAtYN|p~r9x(2acWUX
zW?s4iJOuT)6rcdka4bqMwu1BEDit(xGK)(ztrQ#;Kn#WaG=%_YJSgNRfx@>~p(I}+
zDHUo>YKlT;UP(TZ6A*U!X)@p93`#9dEdfPy5kCV1!!0(5>ol2*1Q-|?iUdJo9B|t~
zA*rCCP$UEr<^+2Vq_s*A7H4oFP{u9x07u#f9tKX{e!ecg36&Q)6t8e7Uf@u?A))wz
zfss=iOmukM;1THYz0RX>iAQ09#6=#hD?C~kz~}=rOi_^t0|Ns})WWk*3nQomK+Qfa
z;0RpCz`(E?9$~1_RLg{-gNX<o%qTjTiO|7P0tzyy3sP8!(7}qLgOvy!Y$!U|h|s}a
z!&bvw#l^r-%TdEx!vtn=*09&GfLUB6tWXCtFfi0G)iBqv)Ns~t)Nrj~$0&&4zDwZ<
zX3*rUG6p9naB2%IO-z9$`Jz-%<^~r~pmYx^qCn*dq)bYN=Pzh^C=4nSWx)A0ouP&y
zma~?zhA|87rin~FT)_-08CEhDfs(*Vrdv#U21ODK3=GhM2%J2BaoK<}NOFE|fnAj<
zEPaBV2THPfHaYppi8;k~dI)9WAbX|2Y4(DM6$CBtm=G}~3c|i1Vg*TsWI4D9<fE00
zMVz2~C5PtPB2Ya-q%%Q2Em8!T0?rZ?IR=zXizso7B*-1qa7>Xb$RtXf0*X1RgiDb!
zrB0Cs89<d#Q3aVqNvMEQ9U@eU*cliY{4`l^G3RFH-D1v7tSI6G*#Ij0HMx<>eAe8|
zy!eV+V5$;ARo!CEO{@S5CRW6Q1QRRbt8Q_Fa&kPh>A(f%K$;T}Hn=TO%+0{S0IlIb
z6gY<#i7_xRAgT^{)&!;MVo;l_f#E9~gP;f`^KIa~AYye<#OjKO)kOj8D+1OXtak;3
z2aVhfG3gnp6LUMbZ-`3G2%Q+)!3AgBl~7n<y25os=mD(@p$Qi<iY`bLUz8}mB2f%$
zwt!0&KTS4JD=059H#I)~7FT?HZhlH>4v5VYA75CSm;;e10+o(M;A{s9KX9o6)(;6&
zP<y8cM1XLyC%7<cV1U36EKICHj2{?~NDdBG?GFraLWGC4gz*Cdb`oL>$Xx7AVNkqH
zXFzWX2Qz3g-{OF^QIe6W76wpj3RHo7o&+xIYZz-7;^C!y4HI%Zw}!cfAs(b2EMLP?
z!w?VB3Ff7M^nsj(;K6MH%Yhn@2suy=K=42{2!aPHn;|?Vkh`laklNLWh=vxpSzV!{
zP^qI(rI4JTUzC!WmspZotdN+OqL7<dT&(vTYzjE9ae-?oNU;wpZf~)GRNP`JtGLBf
zR(Xr5tm+n1TE#7<w8~pdX;qr6kW2`27eY^wHK@J-`3jVKAnh-3LZ}jig)F?etOio*
z1CICyf<hfEJ#05nDG`YdmLARrLZTfkJ?uA7D6E1V>>cdjxcAd!DKY_-7N#J=3`Bqv
zBe;kG#XcknT7V=#gbK(dpt2N_{6ORf7G74d4-A+D*n%REE|kOti*jf$4?S_UFd+I_
zEeug?DXc9FQS2!KDeNtbQ5-28EeuhdDS|0nEsRlIDS|27EsRmzDLgF<Q9LQUEeuh-
zDSRyqQG6-<!3>&0;DRt2o;ny97(nR;l*&E}fg70Qrq3GI8isgy`mABAVTcFCBshhj
zrq2>kF$@-GV8{Zc05H3Tv4*7z)F))B;s!U!Sdn<lHEal8Ix}h$w}v4ek)lC42$cT7
zW}?^xDxAThDJ&(R#s!$oz_13bqg=xf&kL1e0QYa$5+fNI7*MhZsA&c2B`Or978Iox
zr{;nCIf)AS$t8(-B?_5&;HDv@FQmr>@7|zxb}AG!(sJ?>LA{&`T}1W*^>sikM{q|L
zuC)@OxRRjaDum)Hf{MYl7(zohsE>@@4UozYq0|H1kH)SRTpuEox`4Xc*j4(y1f_mR
zj~$d4G#PKPrj%3`q-ye^6)8oa7%bue<vUOef!o5jnDPp4u_hK3C05?zOHR!zNiB-c
zPm2eaJzOREC5br@Hfw5bK}qE;Hc0337DsYMW=={`YMv%L)`AGtreaWb0T)ToP9d~V
zf+YiZp#;j<#a-Y+>4O-9kT8~FNmQzXrHAW=m`n#t52zLe^CUn#-n()Nb23&4t_WQz
zxgl{S{|?0+X*;z}FkZCuJmGjDAoQf?g~+Ir_7~KnFUrMSk&EeIxxyiRLrS`XrHB8n
zgw!041%hkUE=U+&lrX#^VTe-cqI4)hc^!m5zk$?842V7la{nWlse&mR-0IF3szFwR
zG*SS{lMpp|Ovy})4Eeq_3|XM)1ZSWcMg$wpG=v(4ERb8k>M=}Ht6|7Olqie{HkxS&
zH4ItsG6&tXU<OUbsxQtNsma+2nQ7n*3hVGIWEO+Rkf8ZCDZex?r5KW0Gm9a;H+V*b
zGf*=#+$arLO{Hn25R#D!AN+t0S0sb%PfbB86oOJqN{jN420@bY^K+~eLW)XL3Ahtf
zLZ}xjK%A?i;Fg$EoT`vtl95_eo>`oV7ywzx2FXa^L=7q_A$bT?N)~}Kk|q-*4}c0q
z<i6!Cj^zBjlElos;wnLSQiloofoi?EkfbieASgP8ZHC4itqTIy7X+*?fyostYuGl3
zY>~RC0TqBqeqdk}bYMKda)j-K#u=>(4mlSbaxQ_%3lT9_BH}J2rd~<RxflTx1WQ3Q
z-4KzO;x;2=PSynx`wJrWm%!u-mo;u1LbgO*)PV{>Btb(ej*Q^?7A5tA{0~ZTpZ_s|
zhEs?gDJy0I)dk3-WDV-b0|=mD&>Ez)?81O=P^O3zT+4!U7bxk2#y3$+0}Y4}Z<<dH
z!vc670GS8FsBU0ns9{(D&mS-~AR5_UsC6qIJs4w&DJX5(6sBO1eW2C_$O;sUJ`M>U
ze+0z}SXB)ZsND}{GBk7}*WoE_$YG2+jLpbU0vg7HXk(~h$l`=D8ETkP*sItWKw_XY
z4;D^g2emS)SQ!{V<wOlb7ASLoWx%SLz^Xx24Ok3AH7MJFWsp>ZdYE7_4Ap41)-YtD
znOwt=g{B*Q2ntb3)-c(C2CR@$2Xc7fh+iym1sd|NI`0V{rBz5olvA0w#!)FO`5+^%
zloWk16YvyuZhjcUY{sDBF3`||Hv=QXbkIN>2WX%zg`t!ojA1$>BSQ~o3PUi%N=83$
zCvzoJkvwSFmlfPXKrRk{aoVJ(=A~8?6xmfN!OKI4qmahgU~1Ap1th4K@`B+7Lj%JF
zji3t}K|4%Wh;C3>Aq(Llli-pngk0Ab1%p~e;2uf{$X-zW3Tvx`fyN(*aBvPtBa(yp
zEH3a_ED-G|ox##kdjZTtCczG_;`PleE(Q&Df+{1cA{Pb*22K8=0FY}y)yXa9#GITW
zUr>N@7MCWaWR_*7q!xijAa8M{r<Q=GMH2H$Zm|{=W#*NDyCwvxD(HYc*lf=H<dXQz
zyp+_6Vo=!u8?gr)U&W6WLs1~#O=1L1uY3?<5D@9F?W&y-agkr`3cng?cwX{`u*3w@
z2i*KUwimb+mIg!O=77r)-wPT+7d3*eXarqU4Z5Nlbb(vpfs9IrOONjZVettjPr3OZ
zp|V0{jqU|Li;H{~SNJS0a$8*Cw)nsxz-jh?N1)%Y%kLtO%oQG)3p_G6cm(>ryS!(x
z%;CMrqjrTy?E;S)OpN0qkKz>`#S1)&H)Iq*R&@Aw_<mqxkd&RzKZ}0_2efJDd4pfL
zr*=l<6@K*#{Q7qVL?`4gC|Oc_QNZ+yfawJR)2G_T8(g-y?Z`Z#xHI=cK**JVs0*=i
z7X#uh+QwfriN9fD0hTn}nRg*D^h#jNh4_Svfe9CF6EB)1qDw|!h>O1%7=O_=;i5^x
zg`~7ANjVqt3oa%VfE0QkP~7N$Amm8Yg`nVraUa>4bwU}xFo4JjEK_)AWL^}|_`t!y
zE7;>ZL-Puc%5@&?OFY^uJTCH>Ug0sl07KxEgp%?>{QyvuezszURMLnhY6>H$;I3sv
zuCkzFph+xH!iN;X5OvHYptd_$1YD%UtH~OqreY047QC9OVFb+uU=+V_Q&PaiFHTb#
zq16eznkC6>W>C?N(`;lHGePYk-rbCN&Bk;$={B(7wE@}n$l=Wl_6ea-XT@tart3+!
zfeo(>$gal>8RGoFj@JfEKagt=2gsfj&R_;jt||rRqSV9^$aFclMk`9qgY;}tAsriV
zkqzm6gPNl+LD}vlD3`nhCE=H#^z{;y`293F(EHz@CJUsv1~s*cK&>2bQLo7Y?oZqT
zmH6@C;e}gV#idEG!E>-MczmE3WDlrc4H;331XWZ*@NyeoU(|t$?v>!8`+<nW6t4w>
z6a6|kZ-`1xiC$nhF}{Nf)GvUF-j!Bf;IhOG+$jJRfu0AdE(FJ1l#ab39SiCbPe{DT
zApuhqc~Qk|2g@EdaA$<`LU8Pb_>7CvnOCGUF_eNj8K5r40g)qOCnPU~#9c_pyeOS@
zMLG*ZDcn@4ix%D|Le7L;2#vdtkaa^v_aif-jQ$q}Mj3s$rl)A`jl5{$b^^as;O2w6
zP?8rNf-Z!_UWm^mY(Cumnim~|iE$L%Cb&<+;xGEfUr5ZlkW+@;z6%@@7dRv!B|5w;
zW?*0d#WW~seLjLY=zuh!04d2qB?7oh05<`-*^i`#8D)R~DUTwlVL=|$sA0+mC2NKh
z&{QkJbY_SstUHjxRD(3wP{V>MTgzI*3K}*~%wS|-a7ismEy~TzgXA#y_+LtDg+f|>
z5pugaDX};;MIk>=A+w}dp&-9Fvm`S=579}dq;;Kw*Ae)}8Zz@rtQ7oVqn037!`<Pe
zpsV1NUs95vtD_L2psNs)UjUcWFaZ0`NjErLN5K%rcGXcZg0lU>brg(YY*!t)Mia1D
z2w0yfjP0tUU<PG_^qIrhu9}G8ffQ2-pmwD{sEx@6E@B`<fULQhDd73K%oOl^U1myr
z)h({V(!>;K>kw3kfr=Seix|`vtP+M7G>C@i1W-wn4;j>uf(&Yi%#m7<I8*MTfO-e(
z4FQoUY%@gWNG~XyDSuHw10jN>Moex>)&kc{Vwx)gH)w3pI>36#$n~PG+eJ~g6M`o)
z&SYJPNxbBjc+oTIqG(bF*A0Hj3mlT*`~oh7QN}7jtz2-XNy2Cb6}?7T4#8N%2ui%O
z8B&;##}San?GQak^hIN{8B$nq>H=k4aPx^IUEoXx(uPGBE1F#(AA*!(EM~*vE>c{B
z&kby7ZUChaTy9`TQwypMajE4%QwwUy;!rzR7=5HAg*TW1N33$9nL7|+hPq$|sey$h
zj=9m?i!aQ0(A44!GhXyCLmuaWu4H2>lBi)?0BVhZLmru^VN3xx)zD4pkp+dlCSO%w
z0H`Ah9vm(%DJo4aDJ@EcHJ4!1IgpY>K>=Q<p^gNDTQ>0u8ku>ZnH}fSqN3Ei5(Tg*
zY+?w!UPuE}KY)6#c_oRUssJPj?gk^3dcPPIs(4`o$_j}k3L4;Yu!>1Xp^94*qysz)
zrBTJHsiOdHY1vkB>L@@aYHX`GbRhQHR<S|Vace5*DuAbWtQ4xabQB;{I93WEUEnzy
zD+NDIZg8<$6b<UNvlf8{!_ixLpusZ8ATnqm4BY6u#RYbJJjid*Wp%gMz`l+z3I!Fb
zpd<_LqJzp<q!JjRbumb(1T(ms&cnbX&|~|7fs51dE|1_1VX+xv6Ky^)FbKx7K*+VJ
z8(g;dZt*|hcu~*oik=%p6hfX=1$TU7uLQ(jNJzZsm2|}`2_gd_ucl>R$SJszQ+Xw)
z`a(_Z#k9IBX>|}O2>F$P!G|f3=>r2p5K|b_7Z4Qz;zcq=F?|HdMS}!=n4&>~a3_6b
zV2Ec*0STotWr7r>GG&8!IZU}AUGYqLAVClXu?0eYWnd^`DhH{nV5$Krs9>rC@#>iz
zm_LH#9GSmBD2RFp334HGAT!8?%wZq}LCi5AUMzDQNS6<DJV+2kLDWOYuM7-X%=sX7
z1<WNN1qIC2AYKh~El5`ua~((!L_yR;$gd0x9xT3~aQ9;gV)??r;KvdH;zhDVv3vx{
zMS}!ASfW9K@X+|k%*4z0g@K8e?*l6Fm4ks-2q{KuLB18NWm;OgB4kb22FHc<OX@)u
z!^HQhg4*Z@3U})7(MJ|P8+su){7P`-g{bJ0@n_<ZrLJb?UC1xKl3#hDs`_GP&6Ug=
zWZAC_44zDWAY1*Jf<Ynf&lCmXMKi^K;>nXK79<Fw&~#-m<$~1ZF%^Ln<S|u%c$G|5
zAYB<u)gVC-g{%ub6xbM~<mSkNLc|bEOklksECHI~FuWlmd0j;1l8DNJuoa05qb`aV
zToEz2AY$<211l3RL@~IjfHGJDs^>tN>vIQW{vK(ruF)5{K0(aJVQ&|ph{9`H6j6{}
z-~tuP96Y=ZMbQB-`D+;B5hXigFhdGcEhBc<r8A^3FG3r!26qQQ-Au6UE)20Iwahil
zjRrN$*h_khnGePqX2e>9DqaSLTJVqqLac@vF$e@|OM%re^_bQ$rZ8oLL!(HzhN%WT
ztW?7gkMJ#MT@I?xK+EY_64x;@FgWMu7L-D}-=H~Zg_O+VlA_F{QgDe6ssTV#{h%@v
z+#4@OYO9B2WELx=mF6Xbl_n?VB$wuZSKokjfM&MAYjD7eaX><l$!}0)1+oZUIf3j*
zO-#-Jn-8inKx&~?nSy7DdNF9>d}?xOQE_HjD#E^!jKmVS8KBjrrJ#jBAiW^dQN2-I
z0$Q4#m#&ajl%ET(JCmWc9%xV)$xD9uC8>zT-e3>rrj}&nrzqs+r(~vOrWS*Bt0U&^
z)xq9VhxiRPv9GR>SON*e)MB_B@%a<60s~{R)h|CU6>f$fxCy?Jxd^mV0YShjM0Zf{
z1GLy1y*dP~KLf8_D)Is;0S|YJL!;g|vAEa;DbiO$21mduZ%LtgEIu(gxil9P&G|);
zN)wdsK<y3407)@ufJ6b_Z^2>1R!|js8PbW70au|HxK*AisVvc0;j+eejsK32J+XV@
zPo$j5J(GVyJMN-V{1v764woxD@;Bs^=Hy<HGw5*X@$c~O@c+QZASK=5byrh+jl>3h
z5WJ}AbVbwY0*~4QNm)=4#>yMT_<;dLf(x+=!YU6yMVR3w5yOij##cm)FNhf55RqG;
zdPT|TqKMH1uNxwI8&t0tIUNYS;u?9;DDt96<OLDE9}oB?XE=7%b=2KZb2wmiQO$P(
z%M}6T8%jD8Sf=n#;Ge+%fsH{>bOH;w^n<QC`(U38D)zty0s{jBsBi%BKlgzrc9GYx
zBepiAGo&zL6fK}p1o&7$qYmm4hHQ|T3`Ghx3=6<4m~}83G?9xaSW%}WYgs_6%E1K*
zbKY;{v5FMtY_N?*Of}5G3?(qb7#NVoFc!dz7$j*1?Eb1@%t92jj7WtVBVrm9v&hbZ
zH<Cds;!zvLwX7weX&R{2$WxXz3|XL7A6N!?;&K@?1H)=~7}c`Xu+^}F#_DUBvcYPL
zf^dgni9A>dn5bbuxGRMjF%^s=&Voan6-6ABeZXceV_{%ejo_i|nL>#rjT(k5(3TLe
zE*FMag<AF+_C}c+#zv_c_60~2%wPpz0%I~8Db{-wLGi1}mKX%83}9tAxRguFFG@}Y
zjRz$vlq43Vr<Ne5U`xR>YM{c_DX}OowOH4qv_v5>Cq2I?vm_%IsX#^Tc|eSbC%hB|
zR92Q}WF}`IIV7`K0ba+XAjVxvGEz$t5$m`>D_W4;lw6cql3J9Rp9dbYf%qt~pr9zf
zA~P4<sX=lpM$wyBpqG~dUbzKXDFZI`K?WtKDkLU@SJdT07ScgH1`i*E>xzm%Yfr(Y
zc2OLt=mpO|3KbzrYS0kMEg|=!#Ij6CH!Cs64HWQ@@|P_qGcPeG{T54Jeo-#gQe6vV
zDq433w8{(8=Mi&;MJd=&x1#)9Q0c760Un*Y#Sd{fJVtJ@f`juGH`uH3;PAS|3ZWQ_
zZZW1n3Vl#x0#xXOFr?tG1{M7BAfrM3I6dsP-T{^S(coe7uVM`1l2bBf)PvwfG5srI
z`W@WRxlNS^pvl7p{FnIkFYxO>kdV19p?^t2e*^ykkt4DfCH${Q_)lQHAs{`Y^rC>u
z1nwUo<@E~`FYzm1;8(sOu5>{_>8XIw6wc{<llblm2u<Obp?^h6dqw;eUAqJPR~!N_
zNCsUL2)-f^d?76AN?7Jafy|H0jDp#WUl>4ShvQQL;R!BNycV!5;a(B4ChDSq#T5aI
z4#x-l0uu~7ODB|e)ppc^r#vC<x*?)`S43t;$Q4<G4T4u>jW>i`G4(hhc*WHFM93B2
zgbRsDSA3H%q-I`8Ex9OCdPSu4BP)|+IpY@w5IKPrRBy;Yy)uFOfvn=3*bfY>f`(vX
zg6j=&nG50?52Up&NSobI(Y&tWa7o4Cfc!<3z$+?&GhA=TsIM@+C}S|g`^N)*sS6xZ
z;QA*SUXUQ~ZVV!|+Cdp}OTpM(hdkunsE*95VT4vn(zw?BfGZ_<QH)+G<z22}Y`j>*
z)OdbiDjpnF7!y(z#)MRbG2yPluvS~5;8Gb(psedbT6Ke5h2gBWn4#5{IG#{u!U$Vc
z<SNmHAyyV#fk@V{)UYo=TK5cg9+&`4IHQ^ZZZNR)$kl*0)iNh~VN@OI$i*47#7s&p
zDNjw!Qz$9VhZmH{-3i)P6@=@Z)C%a1P*BwdF7`;NAHW(E%2Ja{KocwYDhP#?(p1oT
zOLdr&)X}PpDn)q53biZ+tprU@%*jbj(G)<cML<JgkU_BuP<>DdB0!T==v4-`+5+4g
z;6t?rT0s=mfD8gvQAM>NRy~Ml01=HK0z7PnSs65gDg$vGl>rAM14FR}q!JKe5EP!m
zF+&~%FA8X05zy>lz01LSfzR+FhtU-dqYE5HHz4a>h^iBQJm8nRz^}EU;)1^OT@kS<
zp40s%`GFW7GxD#<>2Hv~V(4_h|B6fK1=+BRBH>p=!Xc6wbFx-wtkK$#u_fz*tm8!y
zrz;{(ADJ0Nof$teGl{x@ifU08#tDpH*%(A6X0XiRTHvz88@zB-#O8{K%>>2=!Xgt2
zCz{SMon$+~_JNGTbs6JJGR7PH4@6v)@w+19H^Fs+>jOpAC9)qFm_<VwA><6!8`26F
zq;($1C@u)Su4;Kn)pCda35$!WAy-sGF3N<=aQ$(EU!;QvQli7lN(Kf7P(vET{~Q7y
zPf1~{Wv*dDEI>h?P6U-dwJbGEhz>>#BZ3W<WvyXCbPteMM1y77QkZJlkr#ZFfVwf@
zasjz3P{WeKT*bk_P{RtLL20js4Z%ijy4Rp>@UP{lVXk3B9o?v5LgRBlP2*-@sO7BT
ztYN9)05#r0BO$e1HJml9NPO-Z&Kfo(J`akkLE#7XK@C?8cMT7Ew<DNAlO<7@24xy7
zvl}WPufcZ9EGkActr6t_Otle;>fqGm{Ja!+0ZK~ogso797FMvb8D^&;cvJ(%GeS5D
z(uhq2Wqxo=mIG@618UvgVks^xDuE0-fQBSAc_B?s_;Ce>h?5J9ZZW#wVoIsF#gtNc
ziz%h*7H@KCNq$;dd{JUbW@&K|sKmL&4=!Hga}(3^GD}KRQg5+C1S(*RN*JT67}Wg(
zHBRBU2((ka2vjy!Ny7^!%tp5ssBk(A8JQ4c;N|OJe84N#6FtLnj`c-ewGO5m_;?bM
z9lUpW_@^>WVV}V?!*M3-0;UCy3t4CKu3%c>xRUjvgziNFy(<EG8+b4B*j?eV>tK8!
zB|Rf^PToZ+ofVa9>Mlyz9uPd>cu?e`lv@Y?4I#w^K`Rs&hFuiW>0rmld&)1^Q!zt)
zf%*#b3lhc``Ax3yn_S>9fs`AdB{Sf(&cMI`8XE<d8<)UK9MVA}jL2iL(lv|=kk!NV
zfz~C%*{C_0ks;3*c{sv_Ay&JVsfMXhy@nCTa0C-EE5Xri>Cq_x^>LuCK=07jFs6a_
z2{KmA4oWR4%1kXwRY=UqfoGCp*r;N0X%cuO6SO%tHANw@M4<#UhKV)^h|v!O?O#Ve
z=mC~5h%*<K1~u8h74$8(#Dap<yc9_Pj4idIBsDMPmI&O;_{_X`s5@^l=cMLAr>ejw
zKtOt5peTgTQ%wQYwX*Oejm;QWP@*n{Bx)W8QSlB=SV!vykLYzCnM*t}GeYLXUgXif
z!lQiwjG%*OOuWvZZK}M^jIhD88zQpM@&C~4s`i&u?Juf2UQvbCnHNNye%$33nqb*k
z*Hbs6@&doc1r7~7$?z4p`a*7Or=Sh0GlG&HJV`eCB2UC3mJ1-Sv_x#dK@kNdeQ@GI
z5rvOmp;rMl4DpCc0%`RRQllSvvVbX=A%(e?5xeWs8B$m<M)W`{K{>$=VqjpP(PRPK
z2cRkv++>GL7BFL(EP$&8b<B{aL70LWi1!t8sHK3mIf9yoU{8X!T4=IXX+Q^r!J|5m
zk!9F)0H~7y76xa=Dor=|yaHr`A|9d~Ho}f*21AC#!Q(8N+(n3XFsKeH0#!WdH56hj
z9Nh2~gU%hGSx^M(d7>@<n-0orh(QkQ8bd+3?Hf3^Jq0yoEiZCdUE#30z+nYz%HAbv
zq+3E}4(}BStqzwS&koNH&j*rHc*nU}ctaRJFn~xzR+YRiVsJ^s;G&4p6%nHgB1YhG
z?j>3mMRcIU+*`DElwPq9zGxJDQ6w18D7TE-1ePhB6F4VuK2_6LlCh#51TU)DUs1D%
z41>$a!W4kByPqZ_c*skWv8W8R7K$J1c4hDye9#bX5oldy5opa*5op@I2sDdO1R7B)
z0`(-owRRDxoGSuF0C*1s*cM2g4qoF4B0&3~i|>I09<;+5gg>x|uo`_}fD<}gtd5Ky
z7?4N_C04T!3~<7OQGiwd0|T5;5oDG6zyK%I1XvA0vdDyrGOOqZ1~}ozsKqM%fdQ3>
zU{nApLM2=nC0I>AFu(~jK32^S3~)k0lvVBn1116SI>^JCthZS6N^=V;A^kwmm=Ad4
zb`eMtEtr0B*g%dkuq&F$z`y`1Qj0+|kRO;C85wUd2wi}o4{S1wprZ${lOGvC;$Ohz
q2MFoF$j>PDfdP~F2on1OA|SFVB`l0`9~iKcAHm{Zz$B&`u%7`>cLziO
literal 0
HcmV?d00001
diff --git a/tasks/bht/__pycache__/galaxysimulation.cpython-311.pyc b/tasks/bht/__pycache__/galaxysimulation.cpython-311.pyc
index 217cfe8eeb559192f49de5628658625e9ff730ff..1a9755185c2906f2b2ed6f525202579068b0e3dc 100644
GIT binary patch
delta 5767
zcmaD9ygi9;IWI340|NuYDd|J06O=acF*9<pGB7YOGcYiGj+-pa7&<wMNosQ)<0i(*
zEX<0E9Fic#V4TLB$ymb>&kN=-FfgPr1T$!IPIhBf6lG&zU=RkG;K0DZFr6Wjp@tz=
zdh!7lCFL5XEGC$SH0E^1$xJ;m!3--IG?{NPR@`E&yv0~Gc@y(WM!U%|EE-%ZB-Th>
z5HPvGZ#ub)1&^_|*o`gXW?%rB%v4r!i>a*g7E@VO5kCV1LovvA3JMBEAU76?Oy*~G
zV6>YYz^dV}!eWi}1p)I5{1zXWnN_$xFu;f+Q3eKvWJZ|N85kIF`jL=h(kgB-rB&Wy
zO2g$CYc>gxBkd-;uxV(axJUuyA}9fNQ4uIGUxJ+Vl4)`;+e=;rkWP>w$aj;Q+07Vr
zChua`;98)vMDL=M;T0*v$<Nqt)ZY*gnZkcXK;w#l`3}Y_0#+ZG83k-UFf*%hb+BCF
zP`Us{5F0o^9tW8S^E)Rjc?hwBBDaPyg)ti}UsO^8atTyj3R4Pm3QG#>G8P7g)htj!
z28aq32sec#g*k-@T}e7qEmI0x4D)1fPNT`Gocz3?*aB-v;fP_L+{kHSlfsEehVgLe
z8m3@|8paeZh-)S=7D=Zw)H0QT+zZyuz)-`G#X5PTtQe0A!vw}y##-hQ9tMWVKV>zT
zQkW+zatSj-xe_VNDcsAL7#LQAJqFTU%L-GAO_B{E$&<p1rjs3?BnM0~D+N<0X9{aA
z7bM6)ZeU=jVX0xQVXNUn4b)oh8ul9QczEE~aMW<e!(Cp(S;LJSNDDY8KV%fHXQ<&`
z#>l|18m^0xVFF_hGc*wSAc0u)1le|cQ3`P%_83Jm0~F)@DXbQ#UO-dAk|KbSELdup
zQUtLk3n56dP{5unQiNle85nAqQbZ;Pa*8vGOrFmvrjR0nEd`~BLY0bvm5PCs)~ASJ
zQwj<vaET}GcZ<2CJpZK)0|UcLBM_khB3>|qik_FCvh*b=g%+7HFff3#t0r@iDM-cN
z|Ns9(MCCwIpnUdn8c5<dNCKqf7E4-wX?@X4H;^buktQ?P1c)w>fiMAckO@!$FdG#6
zMPQaZxD147L5PCYautE{Wsx{YGpNWZ0(tZnYf)ledTJ4<B)P?uS1{R}TRf7dNEf6C
zl$DB<7#J9CF~%41fJ8v0e-R&ug=Pv%Zenq9kr1dj0x48L0+3=?k%57sN^i0Pi`3-B
z+#>ZLZgDUZ149GD18$)Ssh7B=FK|oWl~i1yv><7P(u$OelEzmgjkiYZNIIZ&AnAn4
zMN|JPrv4Ws0y;czIJjTrvF-4jkl5qX<3A&ELE)U3Iq@qpcQCHWS(AUm!2^%%17Vp9
zJTedJc|~S$T;f%@$gA)`LaM{_uBg-ulPl7?7e)21i0XY{;N{i+2qGph-W3s>;(1*}
z^^%C{f}|D77e(~0i0F|lEGj-F^17(nB~i5nMJo(1iq`915!Ic*bQfJ5ukZtL=?N@%
zrIi;rE+|}~xFY1Dw9yr5qX`^0<P;`wOyNeO1SGEuC|nXySRipxK<kQt)&&8rivn6m
z{%5<yD|dlc?x}#p451l83xpPgKos63UEl*Jt0Z5C=La4JM@A1uJi!l(oO)bw=ShTS
z?0O+S`yk1{Cx3xg4qTZfgK{FIq+(!Tht*?k%nS^ZFYqX{tgmIRnOq~HG`U8Wd-5$l
zzR5>;B{zHXCNoZM;#1^h2jvck<}~KX8~CglXHR~^=ciP|ngwb{fMwE{(aY2thAffE
z0sOk=pauk3O$}?77?{Dpz)-`Q1<E&IE>^t~AieHzn`>CJq`-<HdZl4(OuapGCNTDF
z1eZ`8evnebWU_#Ofi|cdDFPP?GDWrw3=BoeAOch)gL8I~Du}BFBJ3uo30MbfgP1xX
zLLEea3b!IqvkpD~Yl37!o+*MCB}G~wF_0sRWI!xgP*KEIq&xYkfP$DHNCaFyfLc<J
z;-N@*vb3O~q7SIlSkDM53!aKd&d|IlqH;w<<pYBtX9VL%5IH$l@Tm0-1CuTMR}4I^
z7=&Mlh`eGD^?{YiAcnES<qD75f|83oI$t>$6f{?etdYAYZ*fK5;sXOSZv-QloE#&h
zV0psgjMW9Bu!{oWR|LW@2!!7d5V<0tu|n{QfcA9(qe}ut8$7mnACNdAeInyb_66gp
zivrPC1fnkpL{DBTBpG%><BawNeT41{0^uLo!4}AZoC$F(B+<cZSOx|LP<;SOgrA$h
zE#n%7ELapWxG=<;O_mULW@MfmBP^v{!x#?=7_fAYIVd7mGWuz<Lm~+rXR4EXg=I`p
zgLfq}q@)G~>@C)`{G#Mk)Yh&CsMys66}y_LWr;bZi6yD=5M_y=Q1@g2h57>tnK^tP
z7??P%z{F%35ot>k%ZogQ9UeX2GbH9nFUVMuy&(GoD@+B%qoC>q>|O>222jlc_GsSZ
z91%BA8?J^C)(Qdn&jRG5B7aDVfOrR#AOb)Hs+WoaK@J0%1#YrW61A^S2Z<O!oC<Bp
zUKg;xBw&Al<D!7)6#>r<)*A+<7ddo0IC{7zL`;cZps_@If%XDzh&w@M!Q2T7V-Wu{
z6AQRC%LH%DZb6HjRkh4D%oU7Tj0}@^iFPrvPd+NDIyp*=kC9`trI;c=sQv~U3Rc89
zIa5qA9#mf<#K8S55VwXo3uG;b$xy>s#mc}?#mm5u%oGn-(X$#9JT*)yOptci1jZt3
zROd`!j5V3OQB-O2Tk)L9#$tku?34M$weh>hmTdRbFf8B(wHhby5}%K-jTi2V35-1^
zlcObM>Rmwb2QFQBASo4GB7h^*2_ywB8^DDDs7xyYl?dnsNeC!Sut1WrH>j<_mYZ6V
zk)Hx-On}lAsIY)EB0$NtO4&UzC$XZ^IX@?7vKEVgZf0IdYI;#(Noo!#ndm^<6oMB-
zRWEU?Uf@=}A)_$I?*juHuh|C%w#h1zE|a|^JtoUYiEB*Y`oIej<~92uz`!SRL0sb!
zuf|1Q4X7R-USGzK3_O!}N%~LDlaimzDy27BMyhafmXxy5MD7pL41Agwcr`vSfz*9q
z0FjUshg<}JnoQsn_igfiDGSE9$v>om8I33VNZZ@8*RZ4@^{{6%q;SsVM6XeLOnTx$
zsY8>?uZlAuv8W_7IVV-2NN(~rY11HZ$^MH`Q<J3#R8rq!%`8dHEiMAJAd5idw<aqj
zU_qttEzW}c;>?oF{Ji2KP?39!xgfu|7*s-m0vOUez9lg^M@C|DfQ$)S8z?llOzx7A
zW)A}uqnu%r*T|UEgV>C`VT=<wC$LONzAGSlK`dlJ$OSPF+EBd1a7!74eL(m?;1Mwh
z`-Jg{!ZT(N_64z!ivpon1VS%>(FbObVINqTIY9=4NN@;-pj6$Ud;$)_9+u6)vbl_m
zeUlH$sj=0v)vzt#2c=VKh7|S`jy0U<J+T_5BB#mn^3sf4lda`7IoVRUASxzj%PS-C
zrphZba!=kbuOXJggC6`f><dIeR){mC@M4h@1IfuUr0~_S&Ssd)ghhooNX29YIibli
z3OtMwP#!;v3qx${WPb&5P+3{aRm)w=Q^SxYIk`|l5iG)6%U8=^!;mEf7ZIoxtQD$Z
zUmy)K!j++htA@LVCxt&nAVqME&@wg#hSl)=$q4HDg0dZ?%wYv(J}{F3)ZwWWt`(^j
zg_)$wP$OI;QX`rooFcGBWEm?1!)j2P04rx?=xHmF1~b7#mh9w<a$-E-uHyu7q7<E+
zs}$<-G9Q$2nT%638E<hF6y;~7CWF%1%M+m5hb=-ku_Cj$iY*dKM?q+=kj&gvh2WCZ
z0)<;FNjasdFGVJQQc?q#_C*<>QbFz(b53gBEtb60ijrF_nI)-3kW|hFs>xEJ$sAFc
zxPwX@P(P-~4#aW=xrsS5ujCeUK~BjnwzSOj(xTK`+=(eE@x`S{1v&X8w>XMZOX4eX
zGIJrU%ACwxb_la-@-Af=HZG8S)#MAx;%ZSKGr_GMP&4NidvS7N3CKmDw$d$@;*9)q
zO|c@~$s#Hyj6Rb+Rb*qMK^j1vp(0QbfUVb51PWS2Tj>^aa(o4ds00yJ#h_**IBQ}M
zMZ%L;t2ndug0iUX<S#1AjifonKQc&jig&o&6%w7IHBo1ZPKWaiL6Pb5ljIi&ZxERv
ze^JovilAAC(+z%cP@_m?@V$rn{cCSO#wp8QMIto{QVgPze9U6<>+UYB&eF6#PR
z(e?Sjz|Ly}CN2tSPT*U>c|%Hhfy*5IiR=>?Cm2phydkT2Mb>mf-W6Hb11VQzy(e-_
zV4C1Gg9$8dxIyoVtm6TdE3)n|@dpAz)43*b%`lvixIl10+6JZ_$`iOQ3fNr{unV{#
zV0S}7<%)vU4wEYi4ikANFi!}cp?m|X!*~Jng5VV%JD5P~t*;1JUjU=8JPf=7JszE2
zJzgDNon9Zzr8&hvfC8dObn-QI(|SEnB!D`bMZO>ws5w*ws&R@uLEKmn(GDU&Lx<p!
z2U5|1a#9h90QDM*7l4}{4Ga+Yfklj!>jMLvFwkU`{lI`qgfR-Uf*POo5DH9cv9MZw
zV1N@m%&cY~7~lk#1gkIO2L>e4gq0OkCL)o1%B<oa7;q2}*Mr=w$$3jKH?ahiz$bHP
z3QH7#G6HK}X>LI!IA`2qPs%8X&rhDLuPG)9>ZgGH{ENdTH$SB`C)KW~e{z(jVYuoA
b1}x|UM<x@a_y-2;<VUdh7chyb2JA!tr))4b
delta 4599
zcmdm<^eC8bIWI340|NuYq1tt+!ZI8Am>IcP7#J9s85kHot4)?>44r&Sm~V3z<0eM2
zG^R|(8isgYka7@CVF+f><d|H?tSHLLz`!63(qh8Ez%ZR5lc9zoR&+8qtCC0>Q##{h
zrXI0ihLsGOOt%;-ZZTF){>i+O(Q5Ky77exw0vZ?iH78$X!DEIl!i*wr1_qEZOl1|f
zn93@P_!t-%ia~BwP*5lmVqjn>5}usJ>cD6<c{Qts%>@Ci3;fz2n3+|$J}|(DA`u1#
zhGa&N2cQ_I<FQ$sR&k3d4VT5%Y!V;`Sxw%|rlAP8T>)e}lmOda#KgeB@RD&dEBi}c
z8IT4LALN0_zuC<gMJ5|@Xt2!GoTE88h2uv31^)C8%uH%r9V}NklrDe~#1IaUKR`O6
z-e3R;fWqSQ6-H1b)G($nW`pI6qDqkbo5GaBoWin<nSo(7NHthDiaaZ_JSc=27#Px-
zYMD~lVwfi1=QLtT;fP_HtjlFGxsr=Vtf(=ap_ZwH9c&JmsA0%rnY@8Xf+IGgmbrwB
zfno9uMh%G+_7v`AObiUG!A=4xtYv{oN@0;?P2s6!t6`dakxNVh*?BcAHLNvk5a&<c
zz$8}B$k3Bn!xYQ_PL8~Ox0p-H^IrylLf8mIykKTvcnOMgO(saVgFFpp!36&P|Np;8
z8kBfIA^0+dfq@}YK%Id>ld}kv6^g__QsN8@47XT|67$kii)27-ro4h8kow6r{6g03
zMVcT*oFGCDB+D3I!~^1hGRrNN+{EJIVo(SvKma7E$}%u8ROwC5XO)`FC?Ha=07?lT
z85kHE7#?s7O>ph9>#(~ksklICLDCA9i<0_RB=xsS>`>Z~bU@{zq5Bm>_X`pp9iAW9
z7z|7<^5}GUPDt$W>G7Wtxu9@P%$)ce`q)Gt2+LgHk$J!?G9lv<uiOP*xd#$b9iDeZ
zrDmAaUy;_mD5`fwRPO@=JFhvIn80{fL~M%ZbrID|BB~1#R|sDe(Y+#~i&OXm2ZMy{
z1eO~DlGg<kE(s_sl~|#;B5(uCMMcvqil!F@%&rKST>zsSV26AVXW)~+z$*tz&=L&E
zpqK$A3=jr|7Ks1ZY;v!V8e{wBEke5(nbMdh&k?a^?3{c<#7``ZX&EB}!)kcaW@M;g
z$P%8sMpRcD6yG3I7;0FvM8OOO1_rF^#X#z<L9q)~R>PVF%EDkaruv?)35-3R;DpBR
z2Tnsp`jZvJ4AelWM^my0RQeUkGcYg|fg-aAloE;*L0lyeVLG{5%-TgA#MA&0${<1o
zM3{pJczP*P1##i2u1F2Ul>`w|pwf@6NOLl$xPq7fNJJ2nXxJDS7+?u&vc0&Wq9sUl
z5+f)nJr$9hp?Oh6<%)>P2L=I7Z^n-xa`Hs+qgo%B84Y3>J6x{ts4XbD$fNU>jX^<k
zg~%GYi}Ds%<SjlhF!Opdg2~AfB@}or7<pY3@V+A8eL=u`@<9p7`U}Rs7X|#T2>4wP
z@Vg-(az#L8LEseu_3Hw9mjv`SaBSh;5wR!wg1+ZP0k10pUKhaVBRkkoX^=}G?g7O;
zYIzLGbl|8j16L+B3|Zjl01++>u?Dq_H4Isd3=ESKC9N5mCeM|W67Mkt#mP!WKTS4B
zpo8N;aq=xm83TC0tz?E|Qc#fHVol30N=`+sI_yFDQ4^dWU8VS?d_jTe3dxTWGIRJo
zFfeiIf{Dp_QqrdCIv06VIy`#3XGqMEo+15#6($Su5WJLwR;FMNB~3mg<>tr4z)-^&
z4-x=*3B(5Z7|aHFqsR>u%dC)e1WrBfpqPgHqsRl~Fpzd|r8hZU+MX>KB&0KWjkKJw
zvc^RYg$|A$?g<f7q9;U8h%OSD{7_m_=>ikDQec8t3N>f}TQpflX7c1pIWc}vdICGA
zhPi?<i+%DuIYn(yegup46oVXD!<51Vtz2%SnR9irjEvIcbMms2<K+`@na2e;uV(TE
zF7e3%^7EK_u1x+RFOy~sDz!nur^yM4OmI#Bdmogpi$Hm`2$T~LC2Wy5D9BhK@#_RC
zf7x<VOEU6Pia?QF3<_RQ{(uy_4xngO21T=OW?o5ZdQoCYY6K|!W`Q%z18%_!qN<m;
zRWER>-jGq4<M)Aqh1cu@1Iy$K3NDjXl(aM^aDCu}2=SVI;9=ksxgf4_iC5zSuLe{-
zC$BH#2L{f`vlaa(^D1dg_EU;7n!x=*lz~t40<XpgCa65GFC!$az)LU&1_o|udinfh
z@?Rwj#?Z-T%E63Ule?7dwK?XpqZfER+C5>Q=+NZ!tKtkuEGo%N&Pi1$lAg?_V(R;g
zQB#wp2voq_V$Cc`%`Gmn1tla<v7^Zf2`f+;a*MMdzc{lbGe57m$QcxN%mw+y#h^wJ
zD9j+us#_A152{E^UaDfkmJ1516_f9)NVEHZN(WA#$&9Kd^&mDQuMguy&Iv3NlJ5$L
zUJ&zI5OP7x>w=iqhT<KDTgo8p1HuOakBC9o7sL=H$^`-M56mC~KCm)#f(!$Z;NS^C
zsTV+{IVkOa7GnjaXIRT%^8wXFM#jsNmDJS)Vwe~hYT0VpYnY03YS<U>PoBslB9_9I
z!j9VfU}UIapL{@6j0Ge$d4{?sI|n3qC!bPRX6J;mzo{#;bJeiUW|+$~SxQ4=5|e;Y
z2`G-hd4Yi;g&QIVQ4JMj0M|qySAqo=h=N=&iHT3mg(3Dk*mWGWoHd*^te`{<R>A-l
z<*MbbVaO7jtf#5S4e|p|EpH7&mN=ZxSIb|+zCdE~L>^Ijt{UzXo)q2`zBT;ISQ!{r
zgW?Np6xeRb$p<t=!Da~53c~c%3)Bdv2&C|?5nRT?z_1#WE1)`1d?ZxEF`FSpxWp7J
z2_|Z|QbbS#Z#F}UD2fVhEGooMRPfYeQ6i3_gcpkn2^1B4SX4-&sNly`F_#Y<J<{-e
z#>mj~qXblEgMG)qkR?6&gN7InBf|v7o~z*eBsJMZE7YLK7*tlUMd&70WENMkMM7w<
zkj&gvh2WCZ0);A8Fm;P1DW^2G$YAnAEw%b6P$r88X=KW@tzxlKsIV>a194e&6e?}2
zSalT2Dj*bumj>abRoZIG++xm2&AY{tms(MBizTxpwFr`Z*~$`gN>icv7gD!@EU^O-
zpvGtsD6WcZz?oM$Gq2<pb3snYEso;Sq=KCMlHyw&#i=Fn6*-x?5LV^nUTvAlTeSss
zgFq&Ln;W1;#x3^Z<irw?t)TYFE!Lv^(!7*gEX5i5<(fi8nv>sXn=m?0R@9Nz1JyT3
z#V@#o2A9jA6b&w%Z!sswSAdAhVxh?mI-*QPLX&6fIJ1?2ijTXKpXe-Cl;9Nq$RNQf
z-r;grOma%j#Jnkao$ejZ6NGN?i+^BX;#HY^Q8!?+v7So(Qy#$ydNWie=uObO$fJ0L
zNAUuW;tgHHE4nV%b-ga>dR^4@xuWayfq{!x2uxfQ(44@xfb)j9!UB;g^%GepFitSM
zAt`%B(r|;`6-mbfDpw@kC$ddon&1Q$)m@QyMbdUh$`wgxi0A_Wq3K+cxMmp6NL*mH
zBB*{t@dU1m0@haqtS<;y-;kESB5k@M<%+cRM2-o}6M}C*RU0p0UQn^Y2&5F+(7h`l
z`hkf-G>B=2+y@2_4I$TtZZO;sxH4)@6od~U_l6!YJP^1uY7dI&+0+Xm^<fvnBTi<Y
z$%ZI|kXI8@FQjE$$jrK!kbNZqRjR1`LPgbu>Y9s1wO5L2A?hLIR|W=GCNHKB3=G~(
zz97nv$)D*XNHBou3xt9wgpePZnFNKtXfiMf3VnbRU%41~1$sO>y?VSlygI!;7)o%8
ze*i^fk;r5<ebag^P@<9r5iTGC)T{zm&_$q{w<r`ORt_RS#VELhht$iUR!0$t05uJY
z=Ytzw4Ga+Yfkl*+>jMLv&=X^o{lI`q7zwe0nyW(f5DG+Uu(0ZWV1N@m%&eLp7~lk#
zFe|8@Mk2LXSwU3~63M5)D*k~12LW+B$jO?Vw*+$&OF+3TXL6K*kbN>JFS6#9<`z_f
zGv_Vlq>PdxP?Z4bQ-U;s2QwgMfQBx9aoFVMr<CTT+7;DLUT<I+u6Tg~3;Mv}%fu-D
SfdM=D5iI@%Ok%15TLu6W3my~z
diff --git a/tasks/bht/bht_oct.py b/tasks/bht/bht_oct.py
new file mode 100644
index 0000000..c55b43f
--- /dev/null
+++ b/tasks/bht/bht_oct.py
@@ -0,0 +1,362 @@
+import numpy as np
+
+
+class MainApp:
+
+ def __init__(self, particles):
+ """Initialize the MainApp with a root node."""
+ self.rootNode = TreeNode(x=0, y=0, z=0, width=200, height=200, depth=200)
+
+ def BuildTree(self, particles):
+ """Build the Octree by inserting particles.
+
+ Args:
+ particles (list): A list of Particle objects to be inserted into the Octree.
+ """
+ self.ResetTree(particles) # Empty the tree
+ for particle in particles:
+ self.rootNode.insert(particle)
+
+ def ResetTree(self, particles):
+ """Reset the Quadtree by reinitializing the root node."""
+ # Define the size of the rootNode based on the positions of the particles
+ #min_x = min([particle.x for particle in particles])
+ #min_y = min([particle.y for particle in particles])
+ #min_z = min([particle.z for particle in particles])
+ #max_x = max([particle.x for particle in particles])
+ #max_y = max([particle.y for particle in particles])
+ #max_z = max([particle.z for particle in particles])
+ #root_height = max_y - min_y
+ #root_width = max_x - min_x
+ #root_depth = max_z - min_z
+ #self.rootNode = TreeNode(x=min_x, y=min_y, z=min_z, width=root_width, height=root_height, depth=root_depth)
+ self.rootNode = TreeNode(x=0, y=0, z=0, width=200, height=200, depth=200)
+
+
+class Particle:
+
+ def __init__(self, x, y, z, mass):
+ """Initialize a Particle with x, y, z coordinates and mass."""
+ self.x = x
+ self.y = y
+ self.z = z
+ self.mass = mass
+
+ self.vx = 0.0 # Velocity component in x direction
+ self.vy = 0.0 # Velocity component in y direction
+ self.vz = 0.0 # Velocity component in z direction
+
+ self.fx = 0.0 # Force component in x direction
+ self.fy = 0.0 # Force component in y direction
+ self.fz = 0.0 # Force component in z direction
+
+
+class TreeNode:
+
+ def __init__(self, x, y, z, width, height, depth):
+ """Initialize a TreeNode representing a octant in the Octree.
+
+ Args:
+ x (float): x-coordinate of the node.
+ y (float): y-coordinate of the node.
+ z (float): z-coordinate of the node.
+ width (float): Width of the node.
+ height (float): Height of the node.
+ depth (float): Depth of the node.
+ """
+ self.x = x
+ self.y = y
+ self.z = z
+ self.width = width
+ self.height = height
+ self.depth = depth
+ self.particle = None # Particle contained in this node
+ self.center_of_mass = np.array([(x + width) / 2, (y + height) / 2, (z + depth) / 2])
+ self.total_mass = 0 # Total mass contained in the node
+ self.children = np.empty(8, dtype=object) # Children nodes (F-SW, F-SE, F-NW, F-NE, B-SW, B-SE, B-NW, B-NE,)
+
+ def contains(self, particle):
+ """Check if the particle is within the bounds of this node.
+
+ Args:
+ particle (Particle): The particle to be checked.
+
+ Returns:
+ bool: True if the particle is within the node's bounds, False otherwise.
+ """
+ return (self.x <= particle.x <= self.x + self.width and self.y <= particle.y <= self.y + self.height
+ and self.z <= particle.z <= self.z + self.depth)
+
+ def insert(self, particle):
+ """Insert a particle into the Octree.
+
+ Args:
+ particle (Particle): The particle to be inserted.
+
+ Returns:
+ bool: True if the particle is inserted, False otherwise.
+ """
+ if not self.contains(particle):
+ return False # Particle doesn't belong in this node
+
+ if self.particle is None and all(child is None for child in self.children):
+ # If the node is empty and has no children, insert particle here.
+ self.particle = particle
+ #print(f'particle inserted: x={round(self.particle.x,2)}, y={round(self.particle.y,2)}, , z={round(self.particle.z,2)}')
+ return True # Particle inserted in an empty node
+
+ if all(child is None for child in self.children):
+ # If no children exist, create and insert both particles
+ self.subdivide()
+ self.insert(self.particle) # Reinsert existing particle into subnode
+ self.insert(particle) # Insert new particle
+ self.particle = None # Clear particle from this node
+ else:
+ # If the node has children, insert particle in the child node.
+ oct_index = self.get_octant(particle)
+ if self.children[oct_index] is None:
+ print('Missing node:', self.children)
+ # Create a child node if it doesn't exist
+ self.children[oct_index] = TreeNode(self.x + (oct_index % 2) * (self.width / 2),
+ self.y + (oct_index // 2) * (self.height / 2),
+ self.z + (oct_index // 2) * (self.depth / 2), self.width / 2,
+ self.height / 2, self.depth / 2)
+ self.children[oct_index].insert(particle)
+
+ def subdivide(self):
+ """Create the children of the node."""
+ sub_width = self.width / 2
+ sub_height = self.height / 2
+ sub_depth = self.depth / 2
+ self.children[0] = TreeNode(self.x, self.y, self.z, sub_width, sub_height, sub_depth) # B-SW
+ self.children[1] = TreeNode(self.x + sub_width, self.y, self.z, sub_width, sub_height, sub_depth) # B-SE
+ self.children[2] = TreeNode(self.x, self.y + sub_height, self.z, sub_width, sub_height, sub_depth) # B-NW
+ self.children[3] = TreeNode(self.x + sub_width, self.y + sub_height, self.z, sub_width, sub_height,
+ sub_depth) # B-NE
+ self.children[4] = TreeNode(self.x, self.y, self.z + sub_depth, sub_width, sub_height, sub_depth) # T-SW
+ self.children[5] = TreeNode(self.x + sub_width, self.y, self.z + sub_depth, sub_width, sub_height,
+ sub_depth) # T-SE
+ self.children[6] = TreeNode(self.x, self.y + sub_height, self.z + sub_depth, sub_width, sub_height,
+ sub_depth) # T-NW
+ self.children[7] = TreeNode(self.x + sub_width, self.y + sub_height, self.z + sub_depth, sub_width, sub_height,
+ sub_depth) # T-NE
+
+ def get_octant(self, particle):
+ """Determine the octant index for a particle based on its position.
+
+ Args:
+ particle (Particle): The particle to determine the octant index for.
+
+ Returns:
+ int: Octant index
+ B - Bottom, T - Top
+ (0 for B-SW, 1 for B-SE, 2 for B-NW, 3 for B-NE,
+ 4 for T-SW, 5 for T-SE, 6 for T-NW, 7 for T-NE).
+ """
+ # Determine separating planes
+ mid_x = self.x + self.width / 2
+ mid_y = self.y + self.height / 2
+ mid_z = self.z + self.depth / 2
+ quad_index = (particle.x >= mid_x) + 2 * (particle.y >= mid_y) + 4 * (particle.z >= mid_z)
+ return quad_index
+
+ def print_tree(self, depth_=0):
+ """Print the structure of the Octree.
+
+ Args:
+ depth_ (int): Current depth in the tree (for indentation in print).
+ """
+ if self.particle:
+ print(
+ f"{' ' * depth_}Particle at ({round(self.particle.x,2)}, {round(self.particle.y,2)}, {round(self.particle.z,2)}) in Node ({round(self.x,2)}, {round(self.y,2)}, {round(self.z,2)}), width={round(self.width,2)}, height={round(self.height,2)}, depth={round(self.depth,2)}"
+ )
+ else:
+ print(
+ f"{' ' * depth_}Node ({round(self.x,2)}, {round(self.y,2)}, {round(self.z,2)}) - Width: {round(self.width,2)}, Height: {round(self.height,2)}, Depth: {round(self.depth,2)}"
+ )
+ for child in self.children:
+ if child:
+ child.print_tree(depth_ + 2)
+
+ def ComputeMassDistribution(self):
+ """Compute the mass distribution for the tree nodes.
+
+ This function calculates the total mass and the center of mass
+ for each node in the Octree. It's a recursive function that
+ computes the mass distribution starting from the current node.
+
+ Note:
+ This method modifies the 'total_mass' and 'center_of_mass' attributes
+ for each node in the Octree.
+
+ Returns:
+ None
+ """
+ if self.particle is not None:
+ # Node contains only one particle
+ self.center_of_mass = np.array([self.particle.x, self.particle.y, self.particle.z])
+ self.total_mass = self.particle.mass
+ else:
+ # Multiple particles in node
+ total_mass = 0
+ center_of_mass_accumulator = np.array([0.0, 0.0, 0.0])
+
+ for child in self.children:
+ if child is not None:
+ # Recursively compute mass distribution for child nodes
+ child.ComputeMassDistribution()
+ total_mass += child.total_mass
+ center_of_mass_accumulator += child.total_mass * child.center_of_mass
+
+ if total_mass > 0:
+ self.center_of_mass = center_of_mass_accumulator / total_mass
+ self.total_mass = total_mass
+ else:
+ # If total mass is 0 or no child nodes have mass, leave values as default
+ pass
+ #self.center_of_mass = np.array([(x+width)/2, (y+height)/2, (z+depth)/2])
+ #self.total_mass = 0
+
+ def CalculateForceFromTree(self, target_particle, theta=1.0):
+ """Calculate the force on a target particle using the Barnes-Hut algorithm.
+
+ Args:
+ target_particle (Particle): The particle for which the force is calculated.
+ theta (float): The Barnes-Hut criterion for force approximation.
+
+ Returns:
+ np.ndarray: The total force acting on the target particle.
+ """
+
+ total_force = np.array([0.0, 0.0, 0.0])
+
+ if self.particle is not None:
+ # Node contains only one particle
+ if self.particle != target_particle:
+ # Calculate gravitational force between target_particle and node's particle
+ force = self.GravitationalForce(target_particle, self.particle)
+ total_force += force
+ else:
+ if self.total_mass == 0:
+ return total_force
+
+ r = np.linalg.norm(
+ np.array([target_particle.x, target_particle.y, target_particle.z]) - self.center_of_mass)
+ d = max(self.width, self.height, self.depth)
+
+ if d / r < theta:
+ # Calculate gravitational force between target_particle and "node particle" representing cluster
+
+ force = self.GravitationalForce(
+ target_particle,
+ Particle(self.center_of_mass[0], self.center_of_mass[1], self.center_of_mass[2], self.total_mass))
+ total_force += force
+ else:
+ for child in self.children:
+ if child is not None:
+ # Recursively calculate force from child nodes
+ if target_particle is not None: # Check if the target_particle is not None
+ force = child.CalculateForceFromTree(target_particle)
+ total_force += force
+
+ return total_force
+
+ def CalculateForce(self, target_particle, particle, theta=1.0):
+ """Calculate the gravitational force between two particles.
+
+ Args:
+ target_particle (Particle): The particle for which the force is calculated.
+ particle (Particle): The particle exerting the force.
+
+ Returns:
+ np.ndarray: The force vector acting on the target particle due to 'particle'.
+ """
+ force = np.array([0.0, 0.0, 0.0])
+ print('function CalculateForce is called')
+ if self.particle is not None:
+ # Node contains only one particle
+ if self.particle != target_particle:
+ # Calculate gravitational force between target_particle and node's particle
+ force = self.GravitationalForce(target_particle, self.particle)
+ else:
+ if target_particle is not None and particle is not None: # Check if both particles are not None
+ r = np.linalg.norm(
+ np.array([target_particle.x, target_particle.y, target_particle.z]) -
+ np.array([particle.x, particle.y, particle.z]))
+ d = max(self.width, self.height, self.depth)
+
+ if d / r < theta:
+ # Calculate gravitational force between target_particle and particle
+ force = self.GravitationalForce(target_particle, particle)
+ else:
+ for child in self.children:
+ if child is not None:
+ # Recursively calculate force from child nodes
+ force += child.CalculateForce(target_particle, particle)
+ return force
+
+ def GravitationalForce(self, particle1, particle2):
+ """Calculate the gravitational force between two particles.
+
+ Args:
+ particle1 (Particle): First particle.
+ particle2 (Particle): Second particle.
+
+ Returns:
+ np.ndarray: The gravitational force vector between particle1 and particle2.
+ """
+ #G = 6.674 * (10 ** -11) # Gravitational constant
+ G = 1
+
+ dx = particle2.x - particle1.x
+ dy = particle2.y - particle1.y
+ dz = particle2.z - particle1.z
+ cutoff_radius = 5
+ r = max(np.sqrt(dx**2 + dy**2 + dz**2), cutoff_radius)
+
+ force_magnitude = G * particle1.mass * particle2.mass / (r**2)
+ force_x = force_magnitude * (dx / r)
+ force_y = force_magnitude * (dy / r)
+ force_z = force_magnitude * (dz / r)
+
+ return np.array([force_x, force_y, force_z])
+
+ # Helper method to retrieve all particles in the subtree
+ def particles_in_subtree(self):
+ """Retrieve all particles in the subtree rooted at this node.
+
+ Returns:
+ list: A list of particles in the subtree rooted at this node.
+ """
+ particles = []
+
+ if self.particle is not None:
+ particles.append(self.particle)
+ else:
+ for child in self.children:
+ if child is not None:
+ particles.extend(child.particles_in_subtree())
+ return len(particles), particles
+
+ def compute_center_of_mass(self):
+ """Compute the center of mass for the node."""
+ print('Function compute_center_of_mass is called')
+ if self.particle is not None:
+ self.center_of_mass = np.array([self.particle.x, self.particle.y, self.particle.z])
+ self.mass = self.particle.mass
+ else:
+ total_mass = 0
+ center_of_mass_accumulator = np.array([0.0, 0.0, 0.0])
+
+ for child in self.children:
+ if child is not None:
+ child.compute_center_of_mass()
+ total_mass += child.mass
+ center_of_mass_accumulator += child.mass * child.center_of_mass
+
+ if total_mass > 0:
+ self.center_of_mass = center_of_mass_accumulator / total_mass
+ self.mass = total_mass
+ else:
+ self.center_of_mass = np.array([0.0, 0.0, 0.0])
+ self.mass = 0
diff --git a/tasks/bht/galaxysimulation.py b/tasks/bht/galaxysimulation.py
index f8aae98..03d532f 100644
--- a/tasks/bht/galaxysimulation.py
+++ b/tasks/bht/galaxysimulation.py
@@ -1,6 +1,6 @@
import matplotlib.pyplot as plt
import numpy as np
-from bht import MainApp, Particle
+from bht_oct import MainApp, Particle
# Simulation class for galaxy collision using Barnes-Hut and Euler method
@@ -10,33 +10,32 @@ class GalaxyCollisionSimulation:
# Initialize the simulation with a given number of particles
self.num_particles = num_particles
self.particles = self.initialize_particles(initial='random') # Generate particles
- self.barnes_hut = MainApp() # Initialize Barnes-Hut algorithm instance
+ self.barnes_hut = MainApp(self.particles) # Initialize Barnes-Hut algorithm instance
self.barnes_hut.BuildTree(self.particles) # Build the Barnes-Hut tree with particles
self.barnes_hut.rootNode.ComputeMassDistribution() #Compute the center of mass of the tree nodes
self.particle_positions = {
- i: [(particle.x, particle.y)]
+ i: [(particle.x, particle.y, particle.z)]
for i, particle in enumerate(self.particles)
} # Store particle positions for each time step
self.particle_velocities = {
- i: [(particle.vx, particle.vy)]
+ i: [(particle.vx, particle.vy, particle.vz)]
for i, particle in enumerate(self.particles)
} # Store particle velocities for each time step
self.particle_forces = {
- i: [(particle.fx, particle.fy)]
+ i: [(particle.fx, particle.fy, particle.fz)]
for i, particle in enumerate(self.particles)
} # Store particle velocities for each time step
self.system_center_of_mass = {
- 0: (self.barnes_hut.rootNode.center_of_mass[0], self.barnes_hut.rootNode.center_of_mass[1])
+ 0: (self.barnes_hut.rootNode.center_of_mass[0], self.barnes_hut.rootNode.center_of_mass[1],
+ self.barnes_hut.rootNode.center_of_mass[2])
}
def initialize_particles(self, initial='random'):
if initial == 'two':
- particles = [Particle(80, 50, 1000), Particle(20, 50, 1000)]
- #Particle(0.5, 0, 10000),
- #Particle(-0.5, 0, 10000)]
- particles[0].vx, particles[0].vy = 0, 3
- particles[1].vx, particles[1].vy = 0, -3
+ particles = [Particle(60, 50, 40, 1000), Particle(40, 50, 60, 1000)]
+ particles[0].vx, particles[0].vy, particles[0].vz = 5, 3, 3
+ particles[1].vx, particles[1].vy, particles[1].vz = 5, -3, 3
#particles=particles[:2]
return particles
@@ -44,13 +43,34 @@ class GalaxyCollisionSimulation:
# Generate random particles within a specified range
particles = []
for _ in range(self.num_particles):
- x = np.random.uniform(0, 100)
- y = np.random.uniform(0, 100)
- #vx = 0
- #vy = 0
+ x = np.random.uniform(50, 150)
+ y = np.random.uniform(50, 150)
+ z = np.random.uniform(50, 150)
+ vx = np.random.uniform(-5, 5)
+ vy = np.random.uniform(-5, 5)
+ vz = np.random.uniform(-5, 5)
#mass = np.random.uniform(10, 100)
- mass = 10000
- particles.append(Particle(x, y, mass))
+ mass = 1000
+ particle = Particle(x, y, z, mass)
+ particle.vx = vx
+ particle.vy = vy
+ particle.vz = vz
+ particles.append(particle)
+
+ return particles
+
+ elif initial == 'four':
+ particles = [
+ Particle(60, 50, 40, 1000),
+ Particle(40, 50, 60, 1000),
+ Particle(50, 60, 70, 1000),
+ Particle(50, 40, 30, 1000)
+ ]
+ particles[0].vx, particles[0].vy, particles[0].vz = 0, 3, 3
+ particles[1].vx, particles[1].vy, particles[1].vz = 0, -3, 3
+ particles[2].vx, particles[2].vy, particles[2].vz = -3, 0, -3
+ particles[3].vx, particles[3].vy, particles[3].vz = 3, 0, -3
+ #particles=particles[:2]
return particles
def simulate(self, num_steps, time_step):
@@ -66,23 +86,24 @@ class GalaxyCollisionSimulation:
self.integrate(time_step, 'velocity_verlet') # Integrate to update particle positions
self.system_center_of_mass[step] = (self.barnes_hut.rootNode.center_of_mass[0],
- self.barnes_hut.rootNode.center_of_mass[1])
+ self.barnes_hut.rootNode.center_of_mass[1],
+ self.barnes_hut.rootNode.center_of_mass[2])
# Store particle positions, velocities and forces for each time step
for i, particle in enumerate(self.particles):
- self.particle_positions[i].append((particle.x, particle.y))
- self.particle_velocities[i].append((particle.vx, particle.vy))
- self.particle_forces[i].append((particle.fx, particle.fy))
+ self.particle_positions[i].append((particle.x, particle.y, particle.z))
+ self.particle_velocities[i].append((particle.vx, particle.vy, particle.vz))
+ self.particle_forces[i].append((particle.fx, particle.fy, particle.fz))
def evaluate_forces(self):
# Evaluate forces acting on each particle using Barnes-Hut algorithm
for particle in self.particles:
force = self.barnes_hut.rootNode.CalculateForceFromTree(particle)
- particle.fx, particle.fy = force
+ particle.fx, particle.fy, particle.fz = force
def calculate_force(self, particle):
# Calculate force acting on a single particle
force = self.barnes_hut.rootNode.CalculateForceFromTree(particle)
- particle.fx, particle.fy = force
+ particle.fx, particle.fy, particle.fz = force
def integrate(self, time_step, method):
# Integrate to update particle positions based on calculated forces
@@ -90,27 +111,32 @@ class GalaxyCollisionSimulation:
for particle in self.particles:
particle.vx += particle.fx * time_step / particle.mass # Update velocity components
particle.vy += particle.fy * time_step / particle.mass
+ particle.vz += particle.fz * time_step / particle.mass
particle.x += particle.vx * time_step # Update particle positions using velocity
particle.y += particle.vy * time_step
+ particle.z += particle.vz * time_step
elif method == 'velocity_verlet':
for particle in self.particles:
particle.vx += particle.fx * time_step / (2 * particle.mass) # First half-step
particle.vy += particle.fy * time_step / (2 * particle.mass)
+ particle.vz += particle.fz * time_step / (2 * particle.mass)
particle.x += particle.vx * time_step # Full step for position
particle.y += particle.vy * time_step
+ particle.z += particle.vz * time_step
self.calculate_force(particle) # Recalculate force for the particle
particle.vx += particle.fx * time_step / (2 * particle.mass) # Second half-step
particle.vy += particle.fy * time_step / (2 * particle.mass)
+ particle.vz += particle.fz * time_step / (2 * particle.mass)
def print_particle_positions(self):
# Print the positions of all particles for each time step
for i, positions in self.particle_positions.items():
print(f"Particle {i + 1} Positions:")
for step, pos in enumerate(positions):
- print(f" Time Step {step}: Position ({pos[0]}, {pos[1]})")
+ print(f" Time Step {step}: Position ({pos[0]}, {pos[1]}, {pos[2]})")
def display_snapshots(self, num_shots, fix_axes=True):
# Display pictures for each time step
@@ -121,27 +147,30 @@ class GalaxyCollisionSimulation:
#Fixed axis limits
print(f"For {axes_type} axis limits, num_particles={self.num_particles}")
for step in np.arange(0, num_time_steps, int(num_time_steps / num_shots)):
- fig, ax = plt.subplots()
+ fig = plt.figure()
+ ax = fig.add_subplot(111, projection='3d')
if fix_axes:
- ax.set_xlim(0, 100)
- ax.set_ylim(0, 100)
+ ax.set_xlim(0, 200)
+ ax.set_ylim(0, 200)
+ ax.set_zlim(0, 200)
ax.set_xlabel('X-axis')
ax.set_ylabel('Y-axis')
+ ax.set_zlabel('Z-axis')
ax.set_title(f'Time Step {step}')
ax.grid()
- print(f'Step {step}')
+ #print(f'Step {step}')
for particle_index, positions in self.particle_positions.items():
- x, y = positions[step]
- vx, vy = self.particle_velocities[particle_index][step]
- fx, fy = self.particle_forces[particle_index][step]
+ x, y, z = positions[step]
+ vx, vy, vz = self.particle_velocities[particle_index][step]
+ fx, fy, fz = self.particle_forces[particle_index][step]
#mass = self.particles[particle_index].mass
- ax.scatter(x, y, c='blue', s=20, alpha=0.5) # Plot particle position for the current time step
- c_x, c_y = self.system_center_of_mass[step]
- ax.scatter(c_x, c_y, c='orange', s=20)
- print(
- f'i={particle_index}: x={round(x,2)}, y={round(y,2)}, vx={round(vx,2)}, vy={round(vy,2)}, fx={round(fx,2)}, fy={round(fy,2)}'
- )
+ ax.scatter(x, y, z, c='blue', s=20, alpha=0.5) # Plot particle position for the current time step
+ c_x, c_y, c_z = self.system_center_of_mass[step]
+ ax.scatter(c_x, c_y, c_z, c='orange', s=20)
+ #print(
+ #f'i={particle_index}: x={round(x,2)}, y={round(y,2)}, z={round(z,2)}, vx={round(vx,2)}, vy={round(vy,2)}, vz={round(vz,2)}, fx={round(fx,2)}, fy={round(fy,2)}, fz={round(fz,2)}'
+ #)
#print(y)
plt.show()
diff --git a/tasks/bht/print_tree.py b/tasks/bht/print_tree.py
index 9390f83..5ed8745 100644
--- a/tasks/bht/print_tree.py
+++ b/tasks/bht/print_tree.py
@@ -1,5 +1,5 @@
import numpy as np
-from bht import MainApp, Particle
+from tasks.bht.bht_oct import MainApp, Particle
# Function to generate random particles
diff --git a/tasks/bht/simulate.py b/tasks/bht/simulate.py
index c703e42..1de867a 100644
--- a/tasks/bht/simulate.py
+++ b/tasks/bht/simulate.py
@@ -1,5 +1,5 @@
from galaxysimulation import GalaxyCollisionSimulation
-sim = GalaxyCollisionSimulation(num_particles=6)
-sim.simulate(num_steps=7000, time_step=0.001)
-sim.display_snapshots(200, fix_axes=True)
+sim = GalaxyCollisionSimulation(num_particles=4)
+sim.simulate(num_steps=10000, time_step=0.001)
+sim.display_snapshots(200, fix_axes=True)
\ No newline at end of file
diff --git a/tasks/bht/system.py b/tasks/bht/system.py
new file mode 100644
index 0000000..264dcfd
--- /dev/null
+++ b/tasks/bht/system.py
@@ -0,0 +1,57 @@
+"""
+This module contains base class for a n-body gravitational system
+with 2 methods of simulations: direct and approximated force field.
+"""
+
+from collections.abc import Callable
+from dataclasses import dataclass
+
+import numpy as np
+
+
+@dataclass
+class GravitationalSystem:
+ r0: np.ndarray[float] # shape = (N,d)
+ v0: np.ndarray[float] # shape = (N,d)
+ m: np.ndarray[float] # shape = N
+ t: np.ndarray[float] # shape = n_time_steps
+ force: Callable
+ solver: Callable
+
+ def __post_init__(self):
+ "Checking dimensions of inputs"
+ assert self.r0.shape == self.v0.shape
+ assert self.m.shape[0] == self.r0.shape[0]
+
+ solvers = ['verlet']
+ assert self.solver.__name__ in solvers
+
+ def direct_simulation(self):
+ """
+ Using integrator to compute trajector in phase space
+ """
+ p = np.zeros((len(self.t), *self.v0.shape))
+ q = np.zeros((len(self.t), *self.r0.shape))
+ q[0] = self.r0
+ p[0] = self.m * self.v0
+
+ for i in range(1, len(self.t)):
+ dt = self.t[i] - self.t[i - 1]
+ p[i], q[i] = self.solver(self.force, q[i - 1], p[i - 1], self.m, dt)
+
+ return self.t, p, q
+
+ def force_field(self):
+ """
+ Using force field method to compute trajector in phase space
+ """
+ p = np.zeros((len(self.t), *self.v0.shape))
+ q = np.zeros((len(self.t), *self.r0.shape))
+ q[0] = self.r0
+ p[0] = self.m * self.v0
+ for i in range(1, len(self.t)):
+ # something with dt = self.t[i] - self.t[i-1]
+ # p = m*v
+ pass
+
+ return self.t, p, q
--
GitLab