API tools faq
paste
Advertisement
Guest User

Untitled

a guest
May 7th, 2025
24
0
Never
Add comment
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
Python 58.40 KB | None | 0 0
raw download clone embed print report
  1. A 2D Barnes-Hut N-Body Gravitational Simulation in PygameI. Introduction to 2D N-Body Gravitational SimulationA. The N-Body Problem: A Brief OverviewThe N-body problem concerns the prediction of the individual motions of a group of objects (or "bodies") that interact with each other gravitationally. This problem is a cornerstone of classical mechanics and astrophysics. N-body simulations are computational models designed to approximate solutions to this problem, tracking the trajectories of particles under the influence of physical forces such as gravity.1 These simulations are indispensable tools in astrophysics, enabling the study of diverse phenomena ranging from the intricate dance of planets within a solar system to the formation and evolution of vast cosmic structures like galaxy filaments and dark matter halos.1The inherent complexity of the N-body problem stems from the fact that every body in the system exerts a gravitational force on every other body. For a system containing N bodies, a direct calculation of all pairwise forces involves computing N(N-1)/2 interactions.4 The computational effort required for this direct summation scales quadratically with the number of bodies, denoted as O(N2).5 While feasible for small N, this quadratic scaling rapidly renders the direct approach computationally intractable for systems involving thousands, millions, or even billions of particles, which are common in astrophysical research.1 This computational bottleneck has been a primary driver for the development of more efficient, approximate algorithms.B. The Barnes-Hut Algorithm: Efficiently Approximating GravityThe Barnes-Hut algorithm, introduced by Josh Barnes and Piet Hut in 1986, is a highly effective approximation method for N-body simulations that significantly reduces computational cost.5 It achieves this by grouping distant particles and treating their collective gravitational influence as if it emanated from a single, more massive pseudo-particle located at their common center of mass. This approximation dramatically cuts down the number of force calculations needed, reducing the overall complexity from O(N2) to O(NlogN) for typical particle distributions.4At the heart of the Barnes-Hut algorithm in two dimensions is a data structure called a quadtree. A quadtree recursively subdivides the 2D simulation space into four equal quadrants. Each node in the tree represents a specific region of space. If a region contains more than one particle (or exceeds a predefined capacity), it is further subdivided into child quadrants. This process continues until each "leaf" node of the tree contains at most one particle, or a minimum subdivision level is reached. This hierarchical structure allows the algorithm to efficiently identify groups of particles that are "sufficiently far away" from a target particle to be approximated. The criterion for "sufficiently far away" is controlled by a parameter commonly denoted as θ (theta). The ratio of a quadrant's width (s) to the distance (d) from the target particle to the quadrant's center of mass is compared to θ. If s/d<θ, the entire group of particles within that quadrant is treated as a single point mass for force calculation purposes.4 A smaller θ value results in higher accuracy but more computations (as fewer approximations are made), while a larger θ speeds up the simulation at the cost of reduced accuracy. This trade-off is a fundamental characteristic of the Barnes-Hut method.C. Simulation Goals: 2D Gravity with Linear Falloff in PygameThis report details the development of a 2D N-body gravitational simulation that employs the Barnes-Hut algorithm. A key feature of this simulation is its use of a non-standard gravitational force law: a linear distance falloff, where the force is proportional to 1/r, rather than the conventional 1/r2 inverse-square law typical of 3D gravity.10 The simulation is implemented in Python using the Pygame library for visualization and user interaction, providing an accessible platform for observing and experimenting with N-body dynamics under this alternative gravitational model. The objective is to produce a single, runnable script that encapsulates the physics engine, rendering logic, and user-tweakable parameters.II. Core Physical and Algorithmic ConceptsA. The Nature of 2D Gravity: Linear Distance Falloff (F∝1/r)In a hypothetical two-dimensional universe, the behavior of gravity differs significantly from its three-dimensional counterpart. Instead of the familiar inverse-square law (F∝1/r2), the gravitational force between two point masses in 2D is inversely proportional to the linear distance (r) separating them.10 The formula for this force can be expressed as:F=rG2D​m1​m2​​where m1​ and m2​ are the masses of the two bodies, r is the distance between them, and G2D​ is the gravitational constant adapted for a 2D system (its value and units would differ from the 3D gravitational constant).This linear falloff can be understood by considering how gravitational field lines would propagate. In 3D, field lines from a point mass spread out over the surface of an expanding sphere. Since the surface area of a sphere is proportional to r2, the density of field lines (representing field strength) decreases as 1/r2. In 2D, however, field lines would spread out along the circumference of an expanding circle. The circumference of a circle is proportional to r, so the field line density, and thus the force, decreases as 1/r.10 This fundamental difference in the force law has profound implications for orbital mechanics within the simulation.One of the most striking consequences pertains to the velocity required for a stable circular orbit. In 3D gravity (1/r2 force), the orbital velocity is v=GM/R​, where M is the central mass and R is the orbital radius.11 However, for a 1/r force law, the situation changes. By equating the gravitational force (Fg​=G2D​Mm/R) with the centripetal force (Fc​=mv2/R) required for circular motion, we get:$RG2D​Mm​=Rmv2​$Multiplying both sides by R and dividing by m (the mass of the orbiting body) yields:$G2D​M=v2Thus,themagnitudeofthecircularorbitalvelocitybecomes:v=G2D​M​$Remarkably, for a 1/r gravitational force, the circular orbital velocity is independent of the orbital radius R and the mass of the orbiting body m. This is a significant deviation from 3D orbital mechanics and will be crucial when initializing the velocities of particles in the simulation.B. The Barnes-Hut Algorithm in Two DimensionsThe implementation of the Barnes-Hut algorithm for this 2D simulation relies on several key components: the quadtree data structure, methods for building this tree by inserting particles, calculation of mass properties for tree nodes, and a recursive force calculation procedure.1. Quadtree Data StructureA quadtree is a tree data structure in which each internal node has exactly four children, corresponding to the four quadrants: Northwest (NW), Northeast (NE), Southwest (SW), and Southeast (SE).8 It is used to recursively partition the 2D simulation space. Each node in the quadtree represents a rectangular region of this space.The attributes of a QuadtreeNode in this simulation will include:
  2. boundary: A pygame.Rect object defining the spatial extent of the node.
  3. particles_in_node: A list to store Body objects. For a leaf node, this list will contain the single particle within its boundary (or be empty). For an internal node that has just been subdivided, this list is temporarily used before particles are passed to children.
  4. children: A list of four elements, each potentially holding a child QuadtreeNode corresponding to the NW, NE, SW, and SE quadrants. Initially None for leaf nodes.
  5. total_mass: The sum of the masses of all particles contained within this node and all its descendant nodes.
  6. center_of_mass: A pygame.Vector2 representing the mass-weighted average position of all particles within this node and its descendants.
  7. is_leaf: A boolean flag indicating whether the node is a leaf (has no children) or an internal node.
  8. depth: The depth of the node in the tree, with the root node at depth 0.
  9. max_particles_per_node: A constant (typically 1 for Barnes-Hut) defining how many particles a leaf node can hold before it must subdivide.
  10. max_depth: A constant to limit the maximum recursion depth of the tree, preventing infinite subdivision in dense regions.
  11. The initial boundary for the root node of the quadtree must be large enough to encompass all particles in the simulation. This can be determined dynamically at the beginning of each tree construction phase by finding the minimum and maximum x and y coordinates among all particles.14 Child node boundaries are then created by halving the parent's boundary in each dimension.2. Building the Quadtree: Particle Insertion and SubdivisionThe quadtree is typically rebuilt from scratch at each physics timestep because particles move.5 The process begins by creating a root node whose boundary encompasses all particles. Then, each particle is inserted into the root node one by one.The insertion logic for a particle into a given QuadtreeNode is recursive:
  12. If the particle does not fall within the node.boundary, the insertion for this branch terminates.
  13. If the node is a leaf (node.is_leaf is true):
  14.  
  15. Add the particle to node.particles_in_node.
  16. If the number of particles in node.particles_in_node now exceeds max_particles_per_node (e.g., > 1) AND the node.depth is less than max_depth:
  17.  
  18. The node must subdivide. Call a subdivide() method.
  19. The subdivide() method creates four new child QuadtreeNode objects, each representing one of the four quadrants of the current node's boundary.
  20. The current node becomes an internal node (is_leaf = False).
  21. All particles previously in node.particles_in_node (including the newly added one) are then re-inserted into the current node, which will now route them to the appropriate child nodes. node.particles_in_node is then cleared.
  22.  
  23.  
  24.  
  25.  
  26. If the node is an internal node (node.is_leaf is false):
  27.  
  28. Determine which of its four child quadrants the particle belongs to using a helper method (e.g., get_quadrant_for_particle(particle)).
  29. Recursively call the insertion method on that child node.
  30.  
  31.  
  32. The get_quadrant_for_particle method compares the particle's position with the center of the current node's boundary to decide if it falls into the NW, NE, SW, or SE child quadrant.3. Calculating Center of Mass (CoM) and Total Mass for NodesAfter all particles have been inserted into the quadtree for a given timestep, or concurrently during the insertion and subdivision process, the total_mass and center_of_mass for each node must be computed.8 This is also a recursive process, typically performed bottom-up or as part of the tree construction:
  33. If a node is a leaf:
  34.  
  35. If it contains a particle, its total_mass is the particle's mass, and its center_of_mass is the particle's position.
  36. If it's an empty leaf, its total_mass is 0, and its center_of_mass can be undefined or set to its geometric center.
  37.  
  38.  
  39. If a node is an internal node:
  40.  
  41. Its total_mass is the sum of the total_mass of its four children (if a child is None or empty, its contribution is zero).
  42. Its center_of_mass is the weighted average of the center_of_mass of its children, where each child's CoM is weighted by its total_mass. The formula for combining two CoMs (CoM1​,CoM2​) with masses (m1​,m2​) to find the new CoM is:
  43. CoMnew​=m1​+m2​CoM1​⋅m1​+CoM2​⋅m2​​
  44. This is applied iteratively for all children.
  45.  
  46.  
  47. 4. Recursive Force Calculation: Applying the θ (theta) CriterionOnce the quadtree is built and mass distributions are computed, the net gravitational force on each particle is calculated by traversing the tree.4 For a given particle P_target, the force calculation starts from the root of the quadtree and proceeds recursively:Let current_node be the node being examined:
  48. If current_node is None or contains no mass, it exerts no force.
  49. If current_node is a leaf node:
  50.  
  51. If it contains a particle P_source that is not P_target itself, calculate the direct gravitational force between P_target and P_source using the 1/r force law (with softening, discussed later). Add this force to the net force on P_target.
  52.  
  53.  
  54. If current_node is an internal node:
  55.  
  56. Calculate the distance d from P_target to current_node.center_of_mass.
  57. Let s be the width of current_node.boundary.
  58. If s/d<θ (the Barnes-Hut criterion):
  59.  
  60. The current_node is considered "sufficiently far away." Approximate the gravitational force it exerts by treating it as a single point mass equal to current_node.total_mass located at current_node.center_of_mass. Calculate this force on P_target (with softening) and add it to the net force.
  61.  
  62.  
  63. Else (the node is too close for approximation):
  64.  
  65. Recursively call the force calculation procedure for each of current_node.children that are not None. Sum the forces returned from these recursive calls.
  66.  
  67.  
  68.  
  69.  
  70. The θ parameter is crucial: θ=0 would mean no approximation, and the algorithm would behave like a direct N² summation (by always recursing to leaves). A typical value for θ is around 0.5 to 0.7, providing a good balance between speed and accuracy.4C. Leapfrog Integration for Particle MotionTo update the positions and velocities of the particles over time, this simulation employs the Leapfrog integration method, specifically the Kick-Drift-Kick (KDK) variant. Leapfrog integration is a numerical technique well-suited for simulating Hamiltonian systems, such as those governed by gravitational forces, because it is symplectic. This means it tends to conserve energy-like quantities over long simulation times, preventing unphysical drift that can occur with simpler methods like Euler integration or even some Runge-Kutta schemes.16 It is also time-reversible and second-order accurate.16The KDK Leapfrog scheme updates positions and velocities in three steps for a time interval Δt:
  71. Kick (Velocity Half-Step): Update velocities by half a timestep using the acceleration at the current time t.
  72. v(t+Δt/2)=v(t)+a(x(t))2Δt​
  73. Drift (Position Full-Step): Update positions by a full timestep using the new half-step velocities.
  74. x(t+Δt)=x(t)+v(t+Δt/2)Δt
  75. Kick (Velocity Half-Step): Update velocities by another half timestep using the acceleration calculated at the new positions x(t+Δt).
  76. v(t+Δt)=v(t+Δt/2)+a(x(t+Δt))2Δt​
  77. Here, x is position, v is velocity, and a is acceleration (calculated as Force/mass). The simulation uses a fixed physics timestep (Δt), which simplifies the KDK implementation while still retaining its benefits of stability and good long-term energy behavior.16
  78. D. Gravitational SofteningIn N-body simulations, when two particles approach each other very closely, the gravitational force between them can become extremely large, especially with a 1/r or 1/r2 force law, leading to a singularity (F→∞ as r→0). This can cause numerical instability, leading to particles being ejected at unrealistically high velocities or other integration errors.1 Gravitational softening is a technique used to prevent these singularities by modifying the force law at very short distances.Instead of treating particles as true point masses, softening effectively gives them a small, finite size or smooths their mass distribution.19 For a standard 1/r2 gravitational force, a common softening method (Plummer softening) involves replacing r2 in the denominator of the force equation with (r2+ϵ2), where ϵ is the "softening length" or "softening factor".1 The potential becomes Φ=−GM/r2+ϵ2​.For the 1/r force law used in this 2D simulation (F=Gm1​m2​/r), a similar approach is adopted. The distance r in the denominator is replaced by an effective distance that does not go to zero. A common way to implement this is to modify the force magnitude calculation to:F=r2+ϵ2​G2D​m1​m2​​This ensures that as r→0, the force approaches G2D​m1​m2​/ϵ, preventing a division by zero and limiting the maximum force. The term ϵ2 is used here (rather than just ϵ) to maintain consistency with the common Plummer-like softening form and its desirable mathematical properties, ensuring the force smoothly transitions from the softened regime to the 1/r regime as r increases beyond ϵ. This can be interpreted as the particles having a smoothed density profile rather than being true point masses, a more physically plausible model at very small scales within the simulation context.19III. Simulation Architecture with PygameThe simulation is structured to run within the Pygame framework, which handles windowing, event processing, and rendering. Key architectural decisions involve managing simulation constants, implementing a game loop that decouples physics from rendering, and correctly visualizing particle motion through state interpolation.A. Managing Simulation ConstantsFor ease of tweaking and experimentation, simulation parameters are centralized. This follows best practices for code maintainability and readability.20 Constants are defined with uppercase names. While a separate config.py file is often recommended for larger projects 20, for a single-script deliverable, these constants will be grouped at the beginning of the script or within a dedicated configuration class.Table 1: Key Simulation ConstantsConstant NameDescriptionDefault Value (Example)Unit (if applicable)SCREEN_WIDTHWidth of the simulation window.1200pixelsSCREEN_HEIGHTHeight of the simulation window.900pixelsCAPTIONTitle of the simulation window."2D Barnes-Hut N-Body"N/AFPS_RENDERTarget frames per second for rendering.60HzPHYSICS_FPSTarget updates per second for physics calculations.100HzPHYSICS_DTFixed timestep for physics integration ($1.0 / \text{PHYSICS_FPS}$).0.01sTIME_SCALEMultiplier for simulation speed (e.g., 2.0 for 2x speed).1.0dimensionlessG_2DGravitational constant for the 2D simulation.500.0(arbitrary units)CENTRAL_BODY_MASSMass of the large central body.10000.0mass unitsNUM_SMALL_BODIESNumber of smaller orbiting bodies to generate.200countSMALL_BODY_MASS_MINMinimum mass for randomly generated small bodies.1.0mass unitsSMALL_BODY_MASS_MAXMaximum mass for randomly generated small bodies.10.0mass unitsCLOUD_RADIUS_MINMinimum initial orbital radius for the cloud of small bodies.150.0length unitsCLOUD_RADIUS_MAXMaximum initial orbital radius for the cloud of small bodies.400.0length unitsBARNES_HUT_THETAApproximation parameter for Barnes-Hut (s/d<θ).0.7dimensionlessSOFTENING_FACTOR_EPSSoftening length ϵ for gravitational force calculation.5.0length unitsGRAVITY_DAMPING_FACTORMultiplier to scale the gravitational force (e.g., < 1 to weaken).1.0dimensionlessMAX_PARTICLES_PER_NODEMax particles in a quadtree leaf before subdivision.1countMAX_TREE_DEPTHMaximum depth for quadtree recursion.20countBODY_COLOR_CENTRALRGB color for the central body.(255, 255, 0)N/ABODY_COLOR_SMALLBase RGB color for small bodies (can be varied).(200, 200, 255)N/ABG_COLORBackground RGB color of the simulation window.(10, 10, 30)N/AQUADTREE_DRAW_COLORRGB color for drawing quadtree boundaries (for debugging).(50, 50, 50)N/ADRAW_QUADTREEBoolean to toggle quadtree visualization.FalseN/AThis structured approach allows users to easily modify the simulation's fundamental parameters, facilitating understanding and exploration of different N-body scenarios.B. The Game Loop: Decoupling Physics and Rendering (GafferOnGames Style)To ensure consistent physics behavior independent of rendering speed, the simulation employs a game loop structure popularized by Glenn Fiedler ("GafferOnGames").22 This loop decouples physics updates, which occur at a fixed timestep, from rendering, which happens as fast as possible up to a target frame rate.The core mechanism involves an accumulator:
  79. At the beginning of each loop iteration, calculate the real_dt – the actual time elapsed since the last frame.
  80. This real_dt is multiplied by TIME_SCALE and added to an accumulator variable. TIME_SCALE allows the simulation to run faster or slower than real-time; a value greater than 1 speeds up the simulation's perceived passage of time by causing more physics steps to be processed per unit of real time.
  81. A while loop then checks if accumulator >= PHYSICS_DT (the fixed timestep for physics).
  82.  
  83. If true, the previous state of all dynamic bodies (e.g., position) is stored for later interpolation.
  84. A single physics update step is performed using the constant PHYSICS_DT. This includes rebuilding the quadtree, calculating forces, and updating positions and velocities via Leapfrog integration.
  85. PHYSICS_DT is subtracted from accumulator.
  86. This inner loop continues until accumulator < PHYSICS_DT, ensuring that physics catches up if rendering lags, or waits if rendering is faster.
  87.  
  88.  
  89. After the physics updates, the rendering step occurs. An interpolation factor, alpha = accumulator / PHYSICS_DT, is calculated. This alpha value (between 0.0 and 1.0) represents how far the current frame is into the next potential physics step.
  90. The game state (specifically, particle positions) is rendered using this alpha value to interpolate between the previous_pos and current_pos of each body. This provides smooth visual motion even if the rendering framerate and physics update rate do not align perfectly.
  91. This architecture is crucial for achieving both physically stable and visually smooth simulations.23
  92. C. Visualizing Motion: State InterpolationWhen the physics update rate (e.g., 100 Hz) differs from the rendering rate (e.g., 60 Hz), simply drawing the particles at their latest physics state can lead to visual stutter or judder. State interpolation mitigates this by displaying particles at positions that are smoothly transitioned between their states at the beginning and end of the most recent (or currently progressing) physics interval.22As described above, the interpolation factor alpha is calculated as accumulator / PHYSICS_DT. The rendered position (render_pos) of a particle is then determined by:p​render​=p​current​⋅α+p​previous​⋅(1.0−α)where p​current​ is the particle's position after the last completed physics step, and p​previous​ is its position before that step. This linear interpolation ensures that particles appear to move smoothly across the screen, effectively bridging the temporal gaps between discrete physics calculations.IV. Implementation DetailsThe simulation is implemented in Python using Pygame. The core components include classes for celestial bodies and quadtree nodes, functions for force calculation and physics updates, scene initialization logic, and the main rendering loop.A. Body ClassThis class encapsulates the properties and behaviors of individual celestial bodies.
  93. Attributes:
  94.  
  95. mass: The mass of the body.
  96. pos: A pygame.math.Vector2 object representing the current position (x,y).
  97. vel: A pygame.math.Vector2 object representing the current velocity (vx​,vy​).
  98. prev_pos: A pygame.math.Vector2 object storing the position from the previous physics step, used for rendering interpolation.
  99. color: An RGB tuple for drawing the body.
  100. radius: The radius of the circle used to represent the body, often scaled with mass.
  101. is_stationary: A boolean flag, true if the body (e.g., the central star) should not move.
  102.  
  103.  
  104. Methods:
  105.  
  106. __init__(mass, pos, vel, color, radius, is_stationary=False): Constructor.
  107. draw_interpolated(surface, alpha): Draws the body on the given Pygame surface. If not stationary, it calculates the interpolated position using alpha, self.pos, and self.prev_pos. Otherwise, it draws at self.pos.
  108. The physics update logic (KDK steps) will be handled externally in the main physics update function but will operate on the pos and vel attributes of Body instances.
  109.  
  110.  
  111. B. QuadtreeNode ClassThis class implements the nodes of the Barnes-Hut quadtree.
  112. Attributes:
  113.  
  114. boundary: A pygame.Rect object defining the spatial region this node covers.
  115. particles_in_node: A list to hold Body objects that fall within this node if it's a leaf and has not yet subdivided, or if it has reached maximum depth or minimum particle capacity for subdivision.
  116. children: A list of four elements, [None, None, None, None], which will hold child QuadtreeNode objects for the NW, NE, SW, SE quadrants respectively, once subdivided.
  117. total_mass: The sum of masses of all particles within this node and its descendants. Initialized to 0.
  118. center_of_mass: A pygame.math.Vector2 for the mass-weighted average position. Initialized to pygame.math.Vector2(0, 0).
  119. is_leaf: Boolean, initially True. Becomes False after subdivision.
  120. depth: Integer representing the node's depth in the tree.
  121. MAX_PARTICLES_PER_NODE: Class or instance constant (e.g., 1).
  122. MAX_DEPTH: Class or instance constant (e.g., 20) to prevent excessive recursion.
  123.  
  124.  
  125. Methods:
  126.  
  127. __init__(self, boundary, depth, MAX_PARTICLES_PER_NODE, MAX_DEPTH): Constructor.
  128. get_quadrant_for_particle(self, particle): Takes a Body object. Returns an integer (0-3) indicating which child quadrant (NW, NE, SW, SE) the particle's position falls into relative to the center of self.boundary. Returns -1 if outside.
  129. subdivide(self): Sets self.is_leaf = False. Creates four new child QuadtreeNode objects. The boundaries for these children are calculated by halving self.boundary.width and self.boundary.height and offsetting them appropriately (e.g., NW child uses top-left quarter of parent boundary). Assigns these to self.children.
  130. insert_particle(self, particle):
  131.  
  132. If particle.pos is not within self.boundary (using self.boundary.collidepoint(particle.pos)), return False.
  133. If self.is_leaf is True:
  134. a.  Add particle to self.particles_in_node.
  135. b.  If len(self.particles_in_node) > self.MAX_PARTICLES_PER_NODE and self.depth < self.MAX_DEPTH:
  136. i.  Call self.subdivide().
  137. ii. For each p_moved in self.particles_in_node:
  138. quadrant_idx = self.get_quadrant_for_particle(p_moved)
  139. if quadrant_idx!= -1:
  140. self.children[quadrant_idx].insert_particle(p_moved)
  141. iii. Clear self.particles_in_node.
  142. c. Return True.
  143. Else (node is internal):
  144. a.  quadrant_idx = self.get_quadrant_for_particle(particle)
  145. b.  If quadrant_idx!= -1:
  146. return self.children[quadrant_idx].insert_particle(particle)
  147. c.  Return False (should not happen if initial check passed and get_quadrant_for_particle is correct).
  148.  
  149.  
  150. calculate_mass_distribution(self):
  151.  
  152. If self.is_leaf:
  153. a.  If self.particles_in_node is not empty:
  154. self.total_mass = sum(p.mass for p in self.particles_in_node)
  155. If self.total_mass > 0:
  156. com_x = sum(p.pos.x * p.mass for p in self.particles_in_node) / self.total_mass
  157. com_y = sum(p.pos.y * p.mass for p in self.particles_in_node) / self.total_mass
  158. self.center_of_mass = pygame.math.Vector2(com_x, com_y)
  159. Else:
  160. self.center_of_mass = pygame.math.Vector2(self.boundary.center)
  161. Else:
  162. self.total_mass = 0
  163. self.center_of_mass = pygame.math.Vector2(self.boundary.center)
  164. Else (internal node):
  165. a.  self.total_mass = 0
  166. b.  weighted_com_x = 0, weighted_com_y = 0
  167. c.  For each child in self.children:
  168. If child is not None:
  169. child.calculate_mass_distribution() (recursive call)
  170. self.total_mass += child.total_mass
  171. weighted_com_x += child.center_of_mass.x * child.total_mass
  172. weighted_com_y += child.center_of_mass.y * child.total_mass
  173. d.  If self.total_mass > 0:
  174. self.center_of_mass = pygame.math.Vector2(weighted_com_x / self.total_mass, weighted_com_y / self.total_mass)
  175. Else:
  176. self.center_of_mass = pygame.math.Vector2(self.boundary.center)
  177.  
  178.  
  179.  
  180.  
  181. The logic for insert_particle and calculate_mass_distribution is based on principles from various sources like.4C. Force Calculation Function (calculate_force_on_body)This function recursively traverses the quadtree to compute the net gravitational force on a specific body.
  182. Parameters: body_to_calc_for (Body object), current_node (QuadtreeNode object), theta (Barnes-Hut threshold), G_2D (gravitational constant), softening_epsilon (softening length), gravity_damping_factor.
  183. Logic:
  184.  
  185. Initialize net_force = pygame.math.Vector2(0, 0).
  186. If current_node is None or current_node.total_mass == 0, return net_force.
  187. If current_node.is_leaf:
  188.  
  189. For each particle_source in current_node.particles_in_node:
  190.  
  191. If particle_source is the same as body_to_calc_for, continue (a body does not exert force on itself).
  192. direction_vec = particle_source.pos - body_to_calc_for.pos
  193. dist_sq = direction_vec.length_squared()
  194. If dist_sq == 0, continue (avoid singularity if somehow at same exact spot, though softening handles this).
  195. effective_dist_denominator_sq = dist_sq + softening_epsilon**2
  196. force_magnitude = (G_2D * body_to_calc_for.mass * particle_source.mass) / math.sqrt(effective_dist_denominator_sq) (for 1/r force with reff​=r2+ϵ2​)
  197. force_vec = direction_vec.normalize() * force_magnitude
  198. net_force += force_vec
  199.  
  200.  
  201.  
  202.  
  203. Else (internal node):
  204.  
  205.  
  206. direction_vec_node = current_node.center_of_mass - body_to_calc_for.pos
  207.  
  208.  
  209. dist_sq_node = direction_vec_node.length_squared()
  210.  
  211.  
  212. If dist_sq_node == 0 and current_node.total_mass > 0: # Body is at the CoM of an internal node
  213. This case is tricky. It might indicate the body is inside the cluster.
  214. Forcing recursion (opening the node) is safer than applying a potentially huge force.
  215. Or, if the node is very small, it might be treated as a leaf.
  216. For simplicity here, we will force recursion if dist_sq_node is zero for an internal node.
  217. This implies s/d would be infinite, so the node must be opened.
  218. pass # Fall through to recursive calls
  219.  
  220.  
  221. s = current_node.boundary.width (assuming square cells for simplicity in s/d)
  222.  
  223.  
  224. If dist_sq_node > 0 and (s * s) / dist_sq_node < theta**2 (using squared values to avoid sqrt):
  225.  
  226. effective_dist_denominator_sq_node = dist_sq_node + softening_epsilon**2
  227. force_magnitude_node = (G_2D * body_to_calc_for.mass * current_node.total_mass) / math.sqrt(effective_dist_denominator_sq_node)
  228. force_vec_node = direction_vec_node.normalize() * force_magnitude_node
  229. net_force += force_vec_node
  230.  
  231.  
  232.  
  233. Else (node is too close, or body is at CoM of node):
  234.  
  235. For each child_node in current_node.children:
  236. net_force += calculate_force_on_body(body_to_calc_for, child_node, theta, G_2D, softening_epsilon, gravity_damping_factor)
  237.  
  238.  
  239.  
  240.  
  241. Return net_force * gravity_damping_factor. The damping factor modifies the strength of gravity globally.
  242.  
  243.  
  244. This logic draws from.4 The application of the gravity_damping_factor addresses the user's request to include such a term.D. Physics Update Step (perform_physics_update)This function orchestrates a single physics step for all bodies.
  245. Parameters: list_of_bodies, root_node_bounds_func (a function that returns the overall bounding box for the quadtree), dt (fixed timestep), G_2D, BARNES_HUT_THETA, SOFTENING_FACTOR_EPS, GRAVITY_DAMPING_FACTOR, MAX_PARTICLES_PER_NODE, MAX_TREE_DEPTH.
  246. Logic:
  247.  
  248. Create a list of dynamic bodies (those not marked is_stationary).
  249. First KDK Kick & Drift:
  250.  
  251. For each body in dynamic_bodies:
  252.  
  253. If body.mass == 0, continue.
  254. force_old = calculate_force_on_body(body, current_root_node, BARNES_HUT_THETA, G_2D, SOFTENING_FACTOR_EPS, GRAVITY_DAMPING_FACTOR)
  255. accel_old = force_old / body.mass
  256. body.vel += accel_old * (dt / 2.0)
  257. body.prev_pos = body.pos.copy() (Store current position as previous for interpolation)
  258. body.pos += body.vel * dt
  259.  
  260.  
  261.  
  262.  
  263. Rebuild Quadtree:
  264.  
  265. Determine overall_bounds by calling root_node_bounds_func(list_of_bodies).
  266. new_root_node = QuadtreeNode(overall_bounds, 0, MAX_PARTICLES_PER_NODE, MAX_TREE_DEPTH)
  267. For body in list_of_bodies:
  268. new_root_node.insert_particle(body)
  269. new_root_node.calculate_mass_distribution()
  270. Update current_root_node to new_root_node.
  271.  
  272.  
  273. Second KDK Kick:
  274.  
  275. For each body in dynamic_bodies:
  276.  
  277. If body.mass == 0, continue.
  278. force_new = calculate_force_on_body(body, current_root_node, BARNES_HUT_THETA, G_2D, SOFTENING_FACTOR_EPS, GRAVITY_DAMPING_FACTOR)
  279. accel_new = force_new / body.mass
  280. body.vel += accel_new * (dt / 2.0)
  281.  
  282.  
  283.  
  284.  
  285. Return the updated current_root_node.
  286.  
  287.  
  288. The quadtree is rebuilt in each physics step as particle positions change significantly.5E. Scene InitializationThis sets up the initial state of the simulation.
  289. Central Body:
  290.  
  291. A Body object is created with CENTRAL_BODY_MASS, positioned at the screen center (SCREEN_WIDTH / 2, SCREEN_HEIGHT / 2).
  292. Initial velocity is (0, 0).
  293. Set is_stationary = True.
  294.  
  295.  
  296. Cloud of Orbiting Bodies:
  297.  
  298. For i from 0 to NUM_SMALL_BODIES - 1:
  299.  
  300. Generate a random mass m_small between SMALL_BODY_MASS_MIN and SMALL_BODY_MASS_MAX.
  301. Generate a random orbital radius r_orbit between CLOUD_RADIUS_MIN and CLOUD_RADIUS_MAX.
  302. Generate a random angle angle_rad from 0 to 2π.
  303. Calculate initial position (x,y) relative to the central body:
  304. pos_x = central_body.pos.x + r_orbit * math.cos(angle_rad)
  305. pos_y = central_body.pos.y + r_orbit * math.sin(angle_rad)
  306. Calculate orbital velocity magnitude for the 1/r force law:
  307. v_mag = math.sqrt(G_2D * CENTRAL_BODY_MASS)
  308. (Note: This velocity is independent of r_orbit and m_small as derived earlier).
  309. Calculate initial velocity components for circular orbit (e.g., counter-clockwise):
  310. vel_x = -v_mag * math.sin(angle_rad)
  311. vel_y = v_mag * math.cos(angle_rad)
  312. To ensure all rotate in the same direction, the sign convention for vel_x and vel_y (or the direction of v_mag) must be consistent.
  313. Create a new Body object with these properties and add it to the list of bodies. This approach is guided by general particle cloud initialization principles 1 and specific orbital velocity calculations 11, adapted for the 1/r force.
  314.  
  315.  
  316.  
  317.  
  318. F. Rendering LogicThis function handles drawing the simulation state to the screen.
  319. Parameters: surface (Pygame display surface), list_of_bodies, alpha (interpolation factor), quadtree_root (optional, for drawing the tree), draw_quadtree_flag.
  320. Logic:
  321.  
  322. Fill the surface with BG_COLOR.
  323. If draw_quadtree_flag is True and quadtree_root is not None:
  324.  
  325. Implement a recursive function draw_quadtree_node(node, surface):
  326.  
  327. If node is None, return.
  328. Draw node.boundary using pygame.draw.rect(surface, QUADTREE_DRAW_COLOR, node.boundary, 1).
  329. If not node.is_leaf, recursively call on node.children.
  330.  
  331.  
  332. Call draw_quadtree_node(quadtree_root, surface).
  333.  
  334.  
  335. For each body in list_of_bodies:
  336.  
  337. Call body.draw_interpolated(surface, alpha). The radius for drawing can be made proportional to body.mass**(1/3) or body.mass**(1/2) for visual representation of mass.
  338.  
  339.  
  340. (Optional) Display simulation information (elapsed time, number of bodies, etc.) using pygame.font.Font, render, and blit.
  341. Call pygame.display.flip() to update the entire screen.
  342.  
  343.  
  344. Basic Pygame drawing functions are used.26V. The Complete Simulation ScriptThe following Python script integrates all the concepts and implementation details discussed. It provides a runnable 2D Barnes-Hut N-body simulation in Pygame.Pythonimport pygame
  345. import math
  346. import random
  347.  
  348. # --- Simulation Constants ---
  349. class Config:
  350.    SCREEN_WIDTH = 1200
  351.    SCREEN_HEIGHT = 900
  352.    CAPTION = "2D Barnes-Hut N-Body Simulation (1/r Gravity)"
  353.    
  354.    FPS_RENDER = 60
  355.    PHYSICS_FPS = 100 # Target updates per second for physics
  356.    PHYSICS_DT = 1.0 / PHYSICS_FPS # Fixed timestep for physics
  357.    
  358.    TIME_SCALE = 1.0 # Multiplier for simulation speed
  359.  
  360.    G_2D = 500.0 # Gravitational constant for 2D (1/r force)
  361.    
  362.    CENTRAL_BODY_MASS = 10000.0
  363.    NUM_SMALL_BODIES = 200
  364.    SMALL_BODY_MASS_MIN = 5.0
  365.    SMALL_BODY_MASS_MAX = 20.0
  366.    
  367.    CLOUD_RADIUS_MIN = 200.0
  368.    CLOUD_RADIUS_MAX = 400.0
  369.    
  370.    BARNES_HUT_THETA = 0.7 # Approximation parameter (s/d < theta)
  371.    SOFTENING_FACTOR_EPS = 10.0 # Softening length epsilon
  372.    GRAVITY_DAMPING_FACTOR = 1.0 # Multiplier for gravitational force
  373.    
  374.    MAX_PARTICLES_PER_NODE = 1 # Max particles in a quadtree leaf
  375.    MAX_TREE_DEPTH = 20 # Max depth for quadtree recursion
  376.    
  377.    BODY_COLOR_CENTRAL = (255, 255, 0) # Yellow
  378.    BODY_COLOR_SMALL_BASE = (200, 200, 255) # Light blue
  379.    BG_COLOR = (10, 10, 30) # Dark blue
  380.    
  381.    QUADTREE_DRAW_COLOR = (50, 50, 50) # Dark grey for quadtree lines
  382.    INFO_TEXT_COLOR = (200, 200, 200)
  383.  
  384.    DRAW_QUADTREE = False # Toggle quadtree visualization
  385.    DRAW_INFO = True
  386.  
  387. # --- Body Class ---
  388. class Body:
  389.    def __init__(self, mass, pos_x, pos_y, vel_x, vel_y, color, is_stationary=False):
  390.        self.mass = mass
  391.        self.pos = pygame.math.Vector2(pos_x, pos_y)
  392.        self.vel = pygame.math.Vector2(vel_x, vel_y)
  393.        self.prev_pos = pygame.math.Vector2(pos_x, pos_y) # For interpolation
  394.        self.color = color
  395.        self.is_stationary = is_stationary
  396.        
  397.        # Visual radius roughly proportional to cube root of mass (for volume)
  398.        if self.mass > 0:
  399.            self.radius = max(2, int(math.pow(self.mass / 10.0, 1.0/2.2)))
  400.        else:
  401.            self.radius = 1
  402.  
  403.  
  404.    def draw_interpolated(self, surface, alpha):
  405.        if self.is_stationary or self.prev_pos == self.pos : # or alpha is problematic for stationary
  406.            render_pos = self.pos
  407.        else:
  408.            render_pos = self.pos * alpha + self.prev_pos * (1.0 - alpha)
  409.        
  410.        pygame.draw.circle(surface, self.color, (int(render_pos.x), int(render_pos.y)), self.radius)
  411.  
  412. # --- QuadtreeNode Class ---
  413. class QuadtreeNode:
  414.    def __init__(self, boundary_rect, depth):
  415.        self.boundary = boundary_rect
  416.        self.depth = depth
  417.        self.particles_in_node =
  418.        self.children = [None, None, None, None] # NW, NE, SW, SE
  419.        self.total_mass = 0.0
  420.        self.center_of_mass = pygame.math.Vector2(boundary_rect.centerx, boundary_rect.centery)
  421.        self.is_leaf = True
  422.  
  423.    def get_quadrant_for_particle(self, particle_pos):
  424.        mid_x = self.boundary.centerx
  425.        mid_y = self.boundary.centery
  426.        
  427.        if particle_pos.x < mid_x: # West
  428.            if particle_pos.y < mid_y: # North-West
  429.                return 0
  430.            else: # South-West
  431.                return 2
  432.        else: # East
  433.            if particle_pos.y < mid_y: # North-East
  434.                return 1
  435.            else: # South-East
  436.                return 3
  437.  
  438.    def subdivide(self):
  439.        self.is_leaf = False
  440.        cx, cy = self.boundary.centerx, self.boundary.centery
  441.        hw, hh = self.boundary.width / 2, self.boundary.height / 2
  442.        
  443.        # NW, NE, SW, SE
  444.        self.children = QuadtreeNode(pygame.Rect(self.boundary.left, self.boundary.top, hw, hh), self.depth + 1)
  445.        self.children = QuadtreeNode(pygame.Rect(cx, self.boundary.top, hw, hh), self.depth + 1)
  446.        self.children = QuadtreeNode(pygame.Rect(self.boundary.left, cy, hw, hh), self.depth + 1)
  447.        self.children = QuadtreeNode(pygame.Rect(cx, cy, hw, hh), self.depth + 1)
  448.  
  449.    def insert_particle(self, particle):
  450.        if not self.boundary.collidepoint(particle.pos):
  451.            return False
  452.  
  453.        if self.is_leaf:
  454.            self.particles_in_node.append(particle)
  455.            if len(self.particles_in_node) > Config.MAX_PARTICLES_PER_NODE and self.depth < Config.MAX_TREE_DEPTH:
  456.                self.subdivide()
  457.                # Move existing particles in this leaf to children
  458.                for p_moved in self.particles_in_node:
  459.                    quadrant_idx = self.get_quadrant_for_particle(p_moved.pos)
  460.                    if self.children[quadrant_idx] is not None:
  461.                         self.children[quadrant_idx].insert_particle(p_moved)
  462.                self.particles_in_node = # Clear particles from this now internal node
  463.            return True
  464.        else: # Internal node
  465.            quadrant_idx = self.get_quadrant_for_particle(particle.pos)
  466.            if self.children[quadrant_idx] is not None:
  467.                return self.children[quadrant_idx].insert_particle(particle)
  468.            return False # Should not happen if children are initialized
  469.  
  470.    def calculate_mass_distribution(self):
  471.        if self.is_leaf:
  472.            if self.particles_in_node:
  473.                self.total_mass = sum(p.mass for p in self.particles_in_node)
  474.                if self.total_mass > 0:
  475.                    com_x = sum(p.pos.x * p.mass for p in self.particles_in_node) / self.total_mass
  476.                    com_y = sum(p.pos.y * p.mass for p in self.particles_in_node) / self.total_mass
  477.                    self.center_of_mass = pygame.math.Vector2(com_x, com_y)
  478.                else: # Should not happen if particles_in_node is not empty and masses > 0
  479.                    self.center_of_mass = pygame.math.Vector2(self.boundary.centerx, self.boundary.centery)
  480.            else:
  481.                self.total_mass = 0.0
  482.                self.center_of_mass = pygame.math.Vector2(self.boundary.centerx, self.boundary.centery)
  483.        else: # Internal node
  484.            self.total_mass = 0.0
  485.            weighted_com_x = 0.0
  486.            weighted_com_y = 0.0
  487.            for child in self.children:
  488.                if child:
  489.                    child.calculate_mass_distribution() # Recursive call
  490.                    self.total_mass += child.total_mass
  491.                    weighted_com_x += child.center_of_mass.x * child.total_mass
  492.                    weighted_com_y += child.center_of_mass.y * child.total_mass
  493.            
  494.            if self.total_mass > 0:
  495.                self.center_of_mass = pygame.math.Vector2(weighted_com_x / self.total_mass, weighted_com_y / self.total_mass)
  496.            else:
  497.                 self.center_of_mass = pygame.math.Vector2(self.boundary.centerx, self.boundary.centery)
  498.    
  499.    def draw(self, surface):
  500.        pygame.draw.rect(surface, Config.QUADTREE_DRAW_COLOR, self.boundary, 1)
  501.        if not self.is_leaf:
  502.            for child in self.children:
  503.                if child:
  504.                    child.draw(surface)
  505.  
  506. # --- Force Calculation ---
  507. def calculate_force_on_body(body_to_calc_for, current_node, theta, G, softening_eps, damping):
  508.    net_force = pygame.math.Vector2(0, 0)
  509.  
  510.    if current_node is None or current_node.total_mass == 0:
  511.        return net_force
  512.  
  513.    if current_node.is_leaf:
  514.        for particle_source in current_node.particles_in_node:
  515.            if particle_source is body_to_calc_for:
  516.                continue
  517.            
  518.            direction_vec = particle_source.pos - body_to_calc_for.pos
  519.            dist_sq = direction_vec.length_squared()
  520.  
  521.            if dist_sq == 0: # Should be rare with float positions
  522.                continue
  523.  
  524.            # Force = G * m1 * m2 / sqrt(r^2 + eps^2) for 1/r gravity
  525.            effective_dist_denominator_sq = dist_sq + softening_eps**2
  526.            if effective_dist_denominator_sq == 0: continue # Avoid sqrt(0)
  527.  
  528.            force_magnitude = (G * body_to_calc_for.mass * particle_source.mass) / math.sqrt(effective_dist_denominator_sq)
  529.            
  530.            # Normalize direction_vec safely
  531.            if direction_vec.length() > 0:
  532.                 force_vec = direction_vec.normalize() * force_magnitude
  533.                 net_force += force_vec
  534.    else: # Internal node
  535.        direction_vec_node = current_node.center_of_mass - body_to_calc_for.pos
  536.        dist_sq_node = direction_vec_node.length_squared()
  537.  
  538.        if dist_sq_node == 0 and current_node.total_mass > 0: # Body at CoM of internal node
  539.            # Force recursion by making s_div_d > theta
  540.            s_div_d_sq = float('inf')
  541.        elif dist_sq_node == 0: # Empty internal node CoM or other issue
  542.             s_div_d_sq = float('inf') # Effectively force recursion or skip
  543.        else:
  544.            s = current_node.boundary.width
  545.            s_div_d_sq = (s * s) / dist_sq_node
  546.  
  547.        if s_div_d_sq < theta**2 : # s/d < theta  => (s/d)^2 < theta^2
  548.            effective_dist_denominator_sq_node = dist_sq_node + softening_eps**2
  549.            if effective_dist_denominator_sq_node == 0: return net_force * damping # Avoid sqrt(0)
  550.  
  551.            force_magnitude_node = (G * body_to_calc_for.mass * current_node.total_mass) / math.sqrt(effective_dist_denominator_sq_node)
  552.            
  553.            if direction_vec_node.length() > 0:
  554.                force_vec_node = direction_vec_node.normalize() * force_magnitude_node
  555.                net_force += force_vec_node
  556.        else: # Node is too close, recurse
  557.            for child_node in current_node.children:
  558.                if child_node:
  559.                    net_force += calculate_force_on_body(body_to_calc_for, child_node, theta, G, softening_eps, 1.0) # Damping applied at the end
  560.  
  561.    return net_force * damping
  562.  
  563.  
  564. # --- Utility: Get overall bounds for Quadtree root ---
  565. def get_simulation_bounds(bodies_list):
  566.    if not bodies_list:
  567.        return pygame.Rect(0, 0, Config.SCREEN_WIDTH, Config.SCREEN_HEIGHT)
  568.  
  569.    min_x = min(b.pos.x for b in bodies_list)
  570.    max_x = max(b.pos.x for b in bodies_list)
  571.    min_y = min(b.pos.y for b in bodies_list)
  572.    max_y = max(b.pos.y for b in bodies_list)
  573.  
  574.    # Add some padding
  575.    padding = 50
  576.    width = max_x - min_x + 2 * padding
  577.    height = max_y - min_y + 2 * padding
  578.    
  579.    # Ensure minimum size
  580.    min_dim = 100
  581.    width = max(width, min_dim)
  582.    height = max(height, min_dim)
  583.  
  584.    # Make it square for simplicity of s in s/d
  585.    side = max(width, height)
  586.    
  587.    center_x = (min_x + max_x) / 2
  588.    center_y = (min_y + max_y) / 2
  589.    
  590.    return pygame.Rect(center_x - side / 2, center_y - side / 2, side, side)
  591.  
  592.  
  593. # --- Main Simulation Logic ---
  594. def main():
  595.    pygame.init()
  596.    screen = pygame.display.set_mode((Config.SCREEN_WIDTH, Config.SCREEN_HEIGHT))
  597.    pygame.display.set_caption(Config.CAPTION)
  598.    clock = pygame.time.Clock()
  599.    font = None
  600.    if Config.DRAW_INFO:
  601.        try:
  602.            font = pygame.font.SysFont("monospace", 15)
  603.        except:
  604.            font = pygame.font.Font(None, 20) # Fallback font
  605.  
  606.  
  607.    bodies =
  608.  
  609.    # Create central body
  610.    central_body = Body(Config.CENTRAL_BODY_MASS,
  611.                        Config.SCREEN_WIDTH / 2, Config.SCREEN_HEIGHT / 2,
  612.                        0, 0, Config.BODY_COLOR_CENTRAL, is_stationary=True)
  613.    bodies.append(central_body)
  614.  
  615.    # Create cloud of smaller bodies
  616.    for _ in range(Config.NUM_SMALL_BODIES):
  617.        mass = random.uniform(Config.SMALL_BODY_MASS_MIN, Config.SMALL_BODY_MASS_MAX)
  618.        r_orbit = random.uniform(Config.CLOUD_RADIUS_MIN, Config.CLOUD_RADIUS_MAX)
  619.        angle = random.uniform(0, 2 * math.pi)
  620.        
  621.        pos_x = central_body.pos.x + r_orbit * math.cos(angle)
  622.        pos_y = central_body.pos.y + r_orbit * math.sin(angle)
  623.        
  624.        # Orbital velocity for 1/r force: v = sqrt(G * M_central)
  625.        # Direction is tangential
  626.        v_mag = math.sqrt(Config.G_2D * Config.CENTRAL_BODY_MASS)
  627.        
  628.        # Counter-clockwise orbit
  629.        vel_x = -v_mag * math.sin(angle)
  630.        vel_y = v_mag * math.cos(angle)
  631.        
  632.        # Vary color slightly
  633.        color_r = max(0, min(255, Config.BODY_COLOR_SMALL_BASE + random.randint(-20, 20)))
  634.        color_g = max(0, min(255, Config.BODY_COLOR_SMALL_BASE + random.randint(-20, 20)))
  635.        color_b = max(0, min(255, Config.BODY_COLOR_SMALL_BASE + random.randint(-20, 20)))
  636.  
  637.        bodies.append(Body(mass, pos_x, pos_y, vel_x, vel_y, (color_r, color_g, color_b)))
  638.  
  639.    # Game loop variables
  640.    running = True
  641.    accumulator = 0.0
  642.    current_real_time = pygame.time.get_ticks() / 1000.0 # Seconds
  643.  
  644.    # Initial Quadtree
  645.    root_bounds = get_simulation_bounds(bodies)
  646.    qtree_root = QuadtreeNode(root_bounds, 0)
  647.    for body in bodies:
  648.        qtree_root.insert_particle(body)
  649.    qtree_root.calculate_mass_distribution()
  650.    
  651.    simulation_time_elapsed = 0.0
  652.  
  653.    while running:
  654.        new_real_time = pygame.time.get_ticks() / 1000.0
  655.        real_dt = new_real_time - current_real_time
  656.        current_real_time = new_real_time
  657.  
  658.        # Cap frame_dt to prevent spiral of death on lag
  659.        if real_dt > 0.25:
  660.            real_dt = 0.25
  661.            
  662.        accumulator += real_dt * Config.TIME_SCALE
  663.  
  664.        for event in pygame.event.get():
  665.            if event.type == pygame.QUIT:
  666.                running = False
  667.            if event.type == pygame.KEYDOWN:
  668.                if event.key == pygame.K_q: # Quit
  669.                    running = False
  670.                if event.key == pygame.K_t: # Toggle quadtree drawing
  671.                    Config.DRAW_QUADTREE = not Config.DRAW_QUADTREE
  672.                if event.key == pygame.K_i: # Toggle info drawing
  673.                    Config.DRAW_INFO = not Config.DRAW_INFO
  674.  
  675.  
  676.        # Fixed physics updates
  677.        while accumulator >= Config.PHYSICS_DT:
  678.            dynamic_bodies = [b for b in bodies if not b.is_stationary]
  679.  
  680.            # KDK Step 1: First Kick (vel) & Drift (pos)
  681.            for body in dynamic_bodies:
  682.                if body.mass == 0: continue
  683.                force_old = calculate_force_on_body(body, qtree_root, Config.BARNES_HUT_THETA, Config.G_2D, Config.SOFTENING_FACTOR_EPS, Config.GRAVITY_DAMPING_FACTOR)
  684.                accel_old = force_old / body.mass
  685.                
  686.                body.vel += accel_old * (Config.PHYSICS_DT / 2.0)
  687.                body.prev_pos = body.pos.copy() # Store before moving
  688.                body.pos += body.vel * Config.PHYSICS_DT
  689.            
  690.            # Rebuild Quadtree for new positions
  691.            root_bounds = get_simulation_bounds(bodies)
  692.            qtree_root = QuadtreeNode(root_bounds, 0)
  693.            for body in bodies:
  694.                qtree_root.insert_particle(body)
  695.            qtree_root.calculate_mass_distribution()
  696.  
  697.            # KDK Step 2: Second Kick (vel)
  698.            for body in dynamic_bodies:
  699.                if body.mass == 0: continue
  700.                force_new = calculate_force_on_body(body, qtree_root, Config.BARNES_HUT_THETA, Config.G_2D, Config.SOFTENING_FACTOR_EPS, Config.GRAVITY_DAMPING_FACTOR)
  701.                accel_new = force_new / body.mass
  702.                body.vel += accel_new * (Config.PHYSICS_DT / 2.0)
  703.  
  704.            accumulator -= Config.PHYSICS_DT
  705.            simulation_time_elapsed += Config.PHYSICS_DT
  706.  
  707.  
  708.        # Rendering with interpolation
  709.        alpha = accumulator / Config.PHYSICS_DT
  710.        
  711.        screen.fill(Config.BG_COLOR)
  712.  
  713.        if Config.DRAW_QUADTREE and qtree_root:
  714.            qtree_root.draw(screen)
  715.            
  716.        for body in bodies:
  717.            body.draw_interpolated(screen, alpha)
  718.        
  719.        if Config.DRAW_INFO and font:
  720.            num_dynamic_bodies = len([b for b in bodies if not b.is_stationary])
  721.            info_text_lines =
  722.            for i, line in enumerate(info_text_lines):
  723.                text_surface = font.render(line, True, Config.INFO_TEXT_COLOR)
  724.                screen.blit(text_surface, (5, 5 + i * 20))
  725.  
  726.  
  727.        pygame.display.flip()
  728.        clock.tick(Config.FPS_RENDER)
  729.  
  730.    pygame.quit()
  731.  
  732. if __name__ == '__main__':
  733.    main()
  734.  
  735. VI. Guide to Tweaking and ExperimentationThe provided simulation script is designed with numerous constants (defined in the Config class) that can be modified to observe different behaviors and explore various aspects of N-body dynamics under 1/r gravity. Understanding how these parameters influence the simulation is key to effective experimentation.Explanation of Key Constants (from Table 1):
  736. Screen Dimensions & FPS:
  737.  
  738. SCREEN_WIDTH, SCREEN_HEIGHT: Affect the visual area. Larger screens may require adjusting particle counts or cloud radii.
  739. FPS_RENDER: Target rendering frame rate. Higher values provide smoother animation if the system can keep up. Does not affect physics accuracy due to decoupling.
  740. PHYSICS_FPS, PHYSICS_DT: PHYSICS_FPS determines how many times per second the physics state is updated. PHYSICS_DT ($1.0 / \text{PHYSICS_FPS}$) is the fixed time interval for each physics step. Higher PHYSICS_FPS (smaller PHYSICS_DT) increases simulation accuracy and stability, especially for fast-moving particles or strong forces, but demands more computational power.
  741.  
  742.  
  743. Simulation Speed & Gravity:
  744.  
  745. TIME_SCALE: A value of 1.0 runs the simulation in "simulation real-time" relative to PHYSICS_DT. Values > 1.0 speed up the simulation's evolution, while < 1.0 slow it down. This directly impacts how quickly the accumulator fills.
  746. G_2D: The gravitational constant for the 2D, 1/r force law. Increasing G_2D makes gravitational forces stronger, leading to faster orbits, tighter clustering, or more chaotic interactions.
  747. GRAVITY_DAMPING_FACTOR: A multiplier applied to the calculated gravitational force. Values < 1.0 weaken gravity, > 1.0 strengthen it. This can be used to simulate scenarios with weaker/stronger fundamental forces or to stabilize a very energetic system.
  748.  
  749.  
  750. Body Properties:
  751.  
  752. CENTRAL_BODY_MASS: The mass of the large, stationary central body. This heavily influences the orbital velocities of smaller bodies (since v=G2D​Mcentral​​ for 1/r gravity).
  753. NUM_SMALL_BODIES: Controls the complexity and visual density of the simulation. Higher numbers significantly increase computational load, especially for Barnes-Hut tree building and force calculations.
  754. SMALL_BODY_MASS_MIN, SMALL_BODY_MASS_MAX: Define the mass range for the orbiting particles. Mass variation can lead to more diverse interactions if collisions or more complex physics were added. In this simulation, mass primarily affects how strongly a particle attracts others.
  755. CLOUD_RADIUS_MIN, CLOUD_RADIUS_MAX: Determine the initial distribution range of the orbiting particles around the central body.
  756.  
  757.  
  758. Barnes-Hut & Softening Parameters:
  759.  
  760. BARNES_HUT_THETA: The opening angle criterion. Values typically range from 0.5 to 1.0.
  761.  
  762. Smaller θ (e.g., 0.3): More nodes are opened, leading to more direct particle-particle force calculations. This increases accuracy but significantly slows down computation, approaching O(N2) as θ→0.
  763. Larger θ (e.g., 0.9): Fewer nodes are opened; more approximations are made using the node's center of mass. This speeds up computation (closer to O(NlogN)) but can reduce accuracy, especially for interactions within dense clusters.
  764.  
  765.  
  766. SOFTENING_FACTOR_EPS: The softening length ϵ. This prevents forces from becoming infinite during very close encounters.
  767.  
  768. Too small ϵ: May not effectively prevent numerical instabilities if particles get extremely close.
  769. Too large ϵ: Can artificially weaken gravitational forces at moderate distances, altering trajectories significantly from the ideal 1/r behavior. A value around the average particle size or slightly smaller is often a good starting point.
  770.  
  771.  
  772. MAX_PARTICLES_PER_NODE, MAX_TREE_DEPTH: Control the structure of the quadtree. MAX_PARTICLES_PER_NODE = 1 is standard for Barnes-Hut. MAX_TREE_DEPTH prevents runaway recursion in extremely dense particle configurations, though well-chosen boundaries and MAX_PARTICLES_PER_NODE usually manage this.
  773.  
  774.  
  775. Visuals:
  776.  
  777. BODY_COLOR_CENTRAL, BODY_COLOR_SMALL_BASE, BG_COLOR: Aesthetic choices.
  778. QUADTREE_DRAW_COLOR, DRAW_QUADTREE: DRAW_QUADTREE = True will render the quadtree boundaries, which is invaluable for debugging the Barnes-Hut implementation and understanding how space is partitioned. However, it adds rendering overhead.
  779. DRAW_INFO: Toggles the display of simulation statistics.
  780.  
  781.  
  782. Suggestions for Further Exploration:
  783. Visualize Quadtree Dynamics: Enable DRAW_QUADTREE and observe how the tree structure adapts as particles move. Pay attention to how nodes are subdivided in denser regions.
  784. Impact of θ: Systematically vary BARNES_HUT_THETA (e.g., 0.2, 0.5, 0.8, 1.0) while keeping other parameters constant. Observe the trade-off between simulation speed (rendering FPS might be an indicator if physics is the bottleneck) and the qualitative behavior of the particle system (e.g., stability of orbits, formation of clumps).
  785. Effect of Softening: Change SOFTENING_FACTOR_EPS. With a very small ϵ, try to induce a scenario where particles have very close, high-energy encounters and observe if instabilities arise. With a very large ϵ, note how orbits might become "mushy" or less defined.
  786. Initial Conditions:
  787.  
  788. Modify CLOUD_RADIUS_MIN and CLOUD_RADIUS_MAX to create different initial cloud shapes (e.g., a thin ring vs. a diffuse cloud).
  789. Experiment with NUM_SMALL_BODIES. How does performance scale?
  790. Implement different initial velocity profiles, e.g., making all particles initially move towards the center or in a non-orbital pattern.
  791.  
  792.  
  793. Force Law Variation: Modify the calculate_force_on_body function to implement a 1/r2 force law (remembering that orbital velocity calculation in main() would also need to change: v=GMcentral​/rorbit​​). Compare the resulting orbital dynamics. Are orbits more stable or less stable? Do they precess differently? (Note: For 1/r2, stable circular and elliptical orbits are expected. For 1/r, the stability of non-circular orbits can be different).
  794. Performance Profiling: Use Python's cProfile module to identify performance bottlenecks in the simulation, especially with high particle counts. Likely candidates are quadtree construction and force calculation.
  795. Adding Collisions: As a more advanced step, one could implement simple collision detection (e.g., based on particle radii) and response (e.g., merging particles and conserving momentum, or elastic bounces).
  796. Energy Conservation: Implement a function to calculate the total kinetic and potential energy of the system. For a 1/r force, the potential energy between two particles is U=G2D​m1​m2​ln(r). Monitor how well the Leapfrog integrator conserves total energy over time with different PHYSICS_DT values.
  797. By systematically adjusting these parameters and observing the outcomes, users can gain a deeper intuition for N-body dynamics, the workings of the Barnes-Hut algorithm, and the nuances of numerical simulation.VII. ConclusionThis report has detailed the design and implementation of a 2D N-body gravitational simulation utilizing the Barnes-Hut algorithm and rendered with Pygame. The simulation uniquely employs a 1/r linear distance falloff for gravity, distinguishing it from typical 1/r2 models and leading to interesting orbital dynamics, such as orbital velocities independent of radius for circular paths.Key features include a robust physics engine built on Kick-Drift-Kick Leapfrog integration for stable and energy-conservative time evolution, and a fixed physics timestep decoupled from rendering to ensure consistent simulation behavior across different hardware. Visual smoothness is maintained through state interpolation. The Barnes-Hut algorithm, with its quadtree-based spatial decomposition and θ-criterion for force approximation, provides an efficient means to handle systems with a moderate number of particles, reducing computational complexity from O(N2) to approximately O(NlogN). Gravitational softening is incorporated to prevent numerical singularities during close particle encounters.The provision of numerous tweakable constants empowers users to experiment with various physical parameters and algorithmic settings, fostering a deeper understanding of the interplay between gravitational forces, particle dynamics, and approximation techniques. The complete, runnable Python script serves as both a demonstration and a foundation for further exploration into computational astrophysics and simulation methods. The balance between physical fidelity and computational efficiency is a central theme in N-body simulations, and this project offers a practical example of navigating that trade-off.
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement
Public Pastes
Advertisement
We use cookies for various purposes including analytics. By continuing to use Pastebin, you agree to our use of cookies as described in the Cookies Policy.  OK, I Understand
Not a member of Pastebin yet?
Sign Up, it unlocks many cool features!