Skip to article frontmatterSkip to article content
Site not loading correctly?

This may be due to an incorrect BASE_URL configuration. See the MyST Documentation for reference.

Deeltjes model

Zoom je heel ver in, dan zie je deeltjes rond vliegen. Elk met een eigen massa, een eigen snelheid en richting. De deeltjes botsen onderling, wisselen energie uit. We zouden op basis van een botsingsmodel van deeltjes iets moeten kunnen leren over thermodynamica. In de thermodynamica gaat het dan om heel veel deeltjes. Maar laten we beginnen met twee botsende ‘deeltjes’.

Je raadt het misschien al... deze simulatie gaan we na bouwen! Daarbij maken we gebruik van Python classes uit het vorige hoofdstuk en de basics van simulaties geleerd in Q1. Zorg er dus voor dat je weet hoe dit werkt! We maken ook gebruik van plotten, en daarbovenop een animatie. Hoe de animatie precies werkt en hoe je die zelf maakt hoef je niet te weten. Je zou wel in staat moeten zijn de code te lezen.

In de onderstaande cell maken we de ParticleClass aan, en geven we enkele parameters van onze simulatie op.

# Importeren van libraries
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.animation import FuncAnimation

# Maken van de class
class ParticleClass:
    def __init__(self, m, v, r, R):
        self.m = m                         
        self.v = np.array(v, dtype=float)  
        self.r = np.array(r, dtype=float)  
        self.R = np.array(R, dtype=float)  

    def update_position(self):
        self.r += self.v * dt

# Simulation parameters
dt = 0.1                                                            # tijd stap
num_steps = 500                                                     # aantal te nemen stappen
particle = ParticleClass(m=1.0, v=[5.0, 0], r=[0.0, 0.0], R=1.0)    # het maken van ons deeltje

We hebben nu een deeltje met massa, een snelheid, een begin positie en een straal. We hebben ook al de stapgrootte bepaald!

We willen de beweging van dat deeltje straks bestuderen en moeten dus een plot maken:

# Creeer een figuur en de assen
fig, ax = plt.subplots()
ax.set_xlim(-10, 10)
ax.set_ylim(-10, 10)
ax.set_aspect('equal')
ax.set_title("Animatie")
ax.set_xlabel("x")
ax.set_ylabel("y")

# Toon het deeltje als een rode stip
dot, = ax.plot([], [], 'ro', markersize=10); # semicolon to suppress output
<Figure size 640x480 with 1 Axes>

Bij het aanmaken van ons deeltje hebben we het deeltje een beginpositie en snelheid mee gegeven. Als we dan per tijdstap de positie bepalen en deze laten plotten en die plots achter elkaar plakken, dan krijgen we een animatie van het deeltje. Met FuncAnimation wordt die animatie voor ons gedaan.

# Initialisieren van de functie voor de animatie
def init():
    dot.set_data([], [])
    return dot,

# Updaten van de functie voor elk frame
def update(frame):
    particle.update_position()
    dot.set_data([particle.r[0]], [particle.r[1]])
    return dot,

# Creeer de animatie
ani = FuncAnimation(fig, update, frames=range(200), init_func=init, blit=True, interval=50)

# Omdat we werken met Jup. Notebooks (en niet een .py file)
from IPython.display import HTML
HTML(ani.to_jshtml())
Loading...

Frames in FuncAnimation verwijst naar het totaal aantal frames dat gebruikt wordt. Interval naar de snelheid van de animatie, nl. 50 milliseconden ofwel 20 frames per seconde.

Merk op dat als we de laatste cel opnieuw runnen, het deeltje zich niet in de oorsprong bevindt. Dat is een ‘eigenaardigheid’ van Jupyter Notebooks. Het is nu beter om alle code in één cel te plaatsen (hieronder gedaan voor je), zodat we ervoor zorgen dat het deeltje altijd in de oorsprong begint.

#  waarden
dt = 0.1         # time step
num_steps = 500  # number of time steps
particle = ParticleClass(m=1.0, v=[5.0, 0], r=[0.0, 0.0],R=1.0)  

# figuur maken
fig, ax = plt.subplots()
ax.set_xlim(-10, 10)
ax.set_ylim(-10, 10)
ax.set_aspect('equal')
ax.set_title("Particle Animation")
ax.set_xlabel("x")
ax.set_ylabel("y")

# een rode bal maken
dot, = ax.plot([], [], 'ro', markersize=10)

# functie voor animatie
def init():
    dot.set_data([], [])
    return dot,

# positie updaten
def update(frame):
    particle.update_position()
    dot.set_data([particle.r[0]], [particle.r[1]])
    return dot,

# animatie
ani = FuncAnimation(fig, update, frames=range(200), init_func=init, blit=True, interval=50)

# For Jupyter notebook:
from IPython.display import HTML
HTML(ani.to_jshtml())
Loading...
<Figure size 640x480 with 1 Axes>

Een van de dingen die je kunt opmerken, is dat het deeltje niet in zijn doos blijft.


# doorlopende doos

#your code/answer
# Waarden
dt = 0.1         # time step
num_steps = 500  # number of time steps
particle = ParticleClass(m=1.0, v=[5.0, 0], r=[0.0, 0.0],R=1.0)  

# figuur maken
fig, ax = plt.subplots()
ax.set_xlim(-10, 10)
ax.set_ylim(-10, 10)
ax.set_aspect('equal')
ax.set_title("Particle Animation")
ax.set_xlabel("x")
ax.set_ylabel("y")

# een rode bal
dot, = ax.plot([], [], 'ro', markersize=10)

# animatie
def init():
    dot.set_data([], [])
    return dot,

# Update positie
def update(frame):
    particle.update_position()
    
    # Alleen periodieke randvoorwaarden in x-richting nodig
    # Y blijft altijd 0
    
    # X-richting (-10 tot 10)
    if particle.r[0] > 10:
        particle.r[0] = -10 + (particle.r[0] - 10)
    elif particle.r[0] < -10:
        particle.r[0] = 10 + (particle.r[0] + 10)
    
    dot.set_data([particle.r[0]], [particle.r[1]])
    return dot,
#  animatie
ani = FuncAnimation(fig, update, frames=range(200), init_func=init, blit=True, interval=50)

# For Jupyter notebook:
from IPython.display import HTML
HTML(ani.to_jshtml())
Loading...
<Figure size 640x480 with 1 Axes>
# doos met harde wanden

# doos met harde wanden

#your code/answer
# Waarden
dt = 0.1         # time step
num_steps = 500  # number of time steps
particle = ParticleClass(m=1.0, v=[5.0, 0], r=[0.0, 0.0],R=1.0)  

# figuur maken
fig, ax = plt.subplots()
ax.set_xlim(-10, 10)
ax.set_ylim(-10, 10)
ax.set_aspect('equal')
ax.set_title("Particle Animation")
ax.set_xlabel("x")
ax.set_ylabel("y")

# rode bal maken
dot, = ax.plot([], [], 'ro', markersize=10)

# animatie maken
def init():
    dot.set_data([], [])
    return dot,

# Update positie
def update(frame):
    particle.update_position()
    
    #botsen
    if particle.r[0] - particle.R < -10:
        particle.r[0] = -10 + particle.R       
        particle.v[0] *= -1                
    elif particle.r[0] + particle.R > 10:     
        particle.r[0] = 10 - particle.R
        particle.v[0] *= -1
    dot.set_data([particle.r[0]], [particle.r[1]])
    return dot,


# Create animatie
ani = FuncAnimation(fig, update, frames=range(200), init_func=init, blit=True, interval=50)

# For Jupyter notebook:
from IPython.display import HTML
HTML(ani.to_jshtml())
Loading...
<Figure size 640x480 with 1 Axes>
#your code/answer
# doos met harde wanden

#your code/answer
# Waarden
dt = 0.1         # time step
num_steps = 500  # number of time steps
particle = ParticleClass(m=1.0, v=[5.0, 3.0], r=[0.0, 0.0],R=1.0)  

# figuur maken
fig, ax = plt.subplots()
ax.set_xlim(-10, 10)
ax.set_ylim(-10, 10)
ax.set_aspect('equal')
ax.set_title("Particle Animation")
ax.set_xlabel("x")
ax.set_ylabel("y")

# rode bal maken
dot, = ax.plot([], [], 'ro', markersize=10)

# animatie maken
def init():
    dot.set_data([], [])
    return dot,

# Update positie
def update(frame):
    particle.update_position()
    
    #botsen verticalen
    if particle.r[0] - particle.R < -10:
        particle.r[0] = -10 + particle.R       
        particle.v[0] *= -1                
    elif particle.r[0] + particle.R > 10:     
        particle.r[0] = 10 - particle.R
        particle.v[0] *= -1
    #botsen horinzontalen
    if particle.r[1] - particle.R < -10:           
        particle.r[1] = -10 + particle.R
        particle.v[1] *= -1

    elif particle.r[1] + particle.R > 10:     
        particle.r[1] = 10 - particle.R
        particle.v[1] *= -1
    
    dot.set_data([particle.r[0]], [particle.r[1]])
    return dot,


#  animatie maken
ani = FuncAnimation(fig, update, frames=range(200), init_func=init, blit=True, interval=50)

# For Jupyter notebook:
from IPython.display import HTML
HTML(ani.to_jshtml())
Loading...
<Figure size 640x480 with 1 Axes>

Laten we teruggaan naar ons deeltje. Er is een functie om de positie bij te werken, hoewel de snelheid hetzelfde lijkt te blijven... kunnen we de snelheid veranderen door (bijvoorbeeld) de versnelling door de zwaartekracht?

# Maken van de class met versnelling
class ParticleClass:
    def __init__(self, m, v, r, R):
        self.m = m                  # mass of the particle
        self.v = np.array(v, dtype=float)  # velocity vector
        self.r = np.array(r, dtype=float)  # position vector
        self.R = np.array(R, dtype=float)  # radius of the particle

    def update_position(self):
        self.r += self.v * dt
    
    def update_velocity(self):
        self.v[1] -= 9.81 * dt


dt = 0.1         # time step
num_steps = 500  # number of time steps
particle = ParticleClass(m=1.0, v=[5.0, 3.0], r=[0.0, 0.0],R=1.0)  

# Create the figure and axis
fig, ax = plt.subplots()
ax.set_xlim(-10, 10)
ax.set_ylim(-10, 10)
ax.set_aspect('equal')
ax.set_title("Particle Animation")
ax.set_xlabel("x")
ax.set_ylabel("y")

# Create the particle as a red dot
dot, = ax.plot([], [], 'ro', markersize=10)

# Initialization function for animation
def init():
    dot.set_data([], [])
    return dot,

# Update function for each frame
def update(frame):
    particle.update_velocity()
    particle.update_position()
    
    #botsen verticalen
    if particle.r[0] - particle.R < -10:
        particle.r[0] = -10 + particle.R       
        particle.v[0] *= -1                
    elif particle.r[0] + particle.R > 10:     
        particle.r[0] = 10 - particle.R
        particle.v[0] *= -1
    #botsen horinzontalen
    if particle.r[1] - particle.R < -10:           
        particle.r[1] = -10 + particle.R
        particle.v[1] *= -1

    elif particle.r[1] + particle.R > 10:     
        particle.r[1] = 10 - particle.R
        particle.v[1] *= -1
    
    dot.set_data([particle.r[0]], [particle.r[1]])
    return dot,


# Create animation
ani = FuncAnimation(fig, update, frames=range(200), init_func=init, blit=True, interval=50)

# For Jupyter notebook:
from IPython.display import HTML
HTML(ani.to_jshtml())
Loading...
<Figure size 640x480 with 1 Axes>
#your code/answer

Een optie om de simulatie te verbeteren is de tijdstap Δt\Delta t kleiner te maken, maar dan hebben we ook meer geduld meer nodig - het aantal berekeningen schaalt met 1Δt\frac{1}{\Delta t}. Een tweede optie is een meer directe oplossing: We weten dat a=const.a = const. en daarom weten we ook de bewegingsvergelijking van het deeltje!

Daarnaast willen we graag weten waar het deeltje is geweest, dat is in onderstaande code toegevoegd.

import numpy as np
import matplotlib.pyplot as plt
from matplotlib.animation import FuncAnimation

# Maken van de class met versnelling
class ParticleClass:
    def __init__(self, m, v, r, R):
        self.m = m                  
        self.v = np.array(v, dtype=float)  
        self.r = np.array(r, dtype=float)  
        self.R = R  # straal als scalar

    def update_position(self, dt):
        self.r += self.v * dt
    
    def update_velocity(self, a, dt):
        self.v += a * dt

# Simulatie parameters
dt = 0.1         
num_steps = 500  
particle = ParticleClass(m=1.0, v=[0,0], r=[0.0, 0.0], R=1.0)  
a = np.array([0.0, -9.81])  

# creeeren van de plot en de assen
fig, ax = plt.subplots()
ax.set_xlim(-10, 10)
ax.set_ylim(-10, 10)
ax.set_aspect('equal')
ax.set_title("Particle Animation (Elastic Bounces)")
ax.set_xlabel("x")
ax.set_ylabel("y")

# creeeren van ons rode deeltje (zonder track line)
dot, = ax.plot([], [], 'ro', markersize=10)

# initializeren van onze functie voor de animatie
def init():
    dot.set_data([], [])
    return dot,

# Update function for each frame
def update(frame):
    # 1. Eerst versnelling toepassen op snelheid
    particle.update_velocity(a, dt)
    
    # 2. Dan positie updaten met de nieuwe snelheid
    particle.update_position(dt)
    
    # 3. Botsing detecteren en corrigeren (elastisch)
    
    # Linker- en rechterwand
    if particle.r[0] - particle.R < -10:
        particle.r[0] = -10 + particle.R
        particle.v[0] *= -1  # perfect elastisch
    elif particle.r[0] + particle.R > 10:
        particle.r[0] = 10 - particle.R
        particle.v[0] *= -1
    
    # Onder- en bovenwand
    if particle.r[1] - particle.R < -10:
        particle.r[1] = -10 + particle.R
        # Perfect elastische botsing: omkeren en grootte behouden
        particle.v[1] = abs(particle.v[1])
    elif particle.r[1] + particle.R > 10:
        particle.r[1] = 10 - particle.R
        particle.v[1] = -abs(particle.v[1])
    
    # Update dot positie
    dot.set_data([particle.r[0]], [particle.r[1]])
    
    return dot,

# Create animation
ani = FuncAnimation(fig, update, frames=range(num_steps), init_func=init, blit=True, interval=50)

# For Jupyter notebook:
from IPython.display import HTML
HTML(ani.to_jshtml())
Loading...
<Figure size 640x480 with 1 Axes>

We hebben steeds slechts gewerkt met een enkel deeltje. Maar om de simulatie uit het filmpje te maken, hebben we twee deeltjes nodig.

import numpy as np
import matplotlib.pyplot as plt
from matplotlib.animation import FuncAnimation

# Maken van de class met versnelling
class ParticleClass:
    def __init__(self, m, v, r, R, color='blue'):
        self.m = m                  
        self.v = np.array(v, dtype=float)  
        self.r = np.array(r, dtype=float)  
        self.R = float(R)  # straal als scalar
        self.color = color

    def update_position(self, a, dt):
        self.r += self.v * dt + 0.5 * a * dt**2  
    
    def update_velocity(self, a, dt):
        """Update the particle's velocity."""
        self.v += a * dt

# Simulation parameters
dt = 0.1         
num_steps = 500  

# Zwaartekracht
a = np.array([0.0, -9.81])  

# Deeltje 1: met horizontale beginsnelheid
particle1 = ParticleClass(m=1.0, v=[5.0, 0], r=[0.0, 0.0], R=1.0, color='blue')  
# Deeltje 2: zonder beginsnelheid
particle2 = ParticleClass(m=1.0, v=[0, 0], r=[0.0, 0.0], R=1.0, color='red')  

# creeeren van de plot en de assen
fig, ax = plt.subplots()
ax.set_xlim(-10, 10)
ax.set_ylim(-10, 10)
ax.set_aspect('equal')
ax.set_title("Two Particles - Same Start Point")
ax.set_xlabel("x")
ax.set_ylabel("y")

# creeeren van de deeltjes
dot1, = ax.plot([], [], 'bo', markersize=10, label='met horizontale snelheid')
dot2, = ax.plot([], [], 'ro', markersize=10, label='geen initiele snelheid')
ax.legend()

# initializeren van onze functie voor de animatie
def init():
    dot1.set_data([], [])
    dot2.set_data([], [])
    return dot1, dot2

# Update function for each frame
def update(frame):
    # Update beide deeltjes
    for particle in [particle1, particle2]:
        particle.update_position(a, dt)
        particle.update_velocity(a, dt)
    
    # Botsing detectie voor particle1
    if particle1.r[0] - particle1.R < -10:
        particle1.r[0] = -10 + particle1.R
        particle1.v[0] *= -1
    elif particle1.r[0] + particle1.R > 10:
        particle1.r[0] = 10 - particle1.R
        particle1.v[0] *= -1
    
    if particle1.r[1] - particle1.R < -10:
        particle1.r[1] = -10 + particle1.R
        particle1.v[1] *= -1
    elif particle1.r[1] + particle1.R > 10:
        particle1.r[1] = 10 - particle1.R
        particle1.v[1] *= -1
    
    # Botsing detectie voor particle2
    if particle2.r[0] - particle2.R < -10:
        particle2.r[0] = -10 + particle2.R
        particle2.v[0] *= -1
    elif particle2.r[0] + particle2.R > 10:
        particle2.r[0] = 10 - particle2.R
        particle2.v[0] *= -1
    
    if particle2.r[1] - particle2.R < -10:
        particle2.r[1] = -10 + particle2.R
        particle2.v[1] *= -1
    elif particle2.r[1] + particle2.R > 10:
        particle2.r[1] = 10 - particle2.R
        particle2.v[1] *= -1
    
    # Update dot posities
    dot1.set_data([particle1.r[0]], [particle1.r[1]])
    dot2.set_data([particle2.r[0]], [particle2.r[1]])
    
    return dot1, dot2

# Create animation
ani = FuncAnimation(fig, update, frames=range(num_steps), init_func=init, blit=True, interval=50)

# For Jupyter notebook:
from IPython.display import HTML
HTML(ani.to_jshtml())
Loading...
<Figure size 640x480 with 1 Axes>

We waren gebleven bij het maken van een botsingsmodel, waarbij we nu twee deeltjes hebben die onderhevig zijn aan zwaartekracht.

Laten we de zwaartekracht even vergeten en alleen 1D kijken.

import numpy as np
import matplotlib.pyplot as plt
from matplotlib.animation import FuncAnimation

# Maken van de class met botsing
class ParticleClass:
    def __init__(self, m, v, r, R):
        self.m = m                  
        self.v = np.array(v, dtype=float)  
        self.r = np.array(r, dtype=float)  
        self.R = float(R)  # straal als scalar

    def update_position(self):
        self.r += self.v * dt  

    def collide_detection(self, other):
        dx = self.r[0] - other.r[0]
        dy = self.r[1] - other.r[1]
        rr = self.R + other.R
        return dx**2 + dy**2 < rr**2

# Simulation parameters
dt = 0.1
num_steps = 200

particleA = ParticleClass(m=1.0, v=[2.5, 0], r=[-2.0, 0.0], R=0.45)  
particleB = ParticleClass(m=1.0, v=[-1, 0], r=[0.0, 0.0], R=0.45)  

fig, ax = plt.subplots()
ax.set_xlim(-10, 10)
ax.set_ylim(-10, 10)
ax.set_aspect('equal')
ax.set_title("Particle Animation")
ax.set_xlabel("x")
ax.set_ylabel("y")

dotA, = ax.plot([], [], 'ro', markersize=10)
dotB, = ax.plot([], [], 'bo', markersize=10)

def init():
    dotA.set_data([], [])
    dotB.set_data([], [])
    return dotA, dotB

def update(frame):
    particleA.update_position()
    particleB.update_position()
    
    dotA.set_data([particleA.r[0]], [particleA.r[1]])
    dotB.set_data([particleB.r[0]], [particleB.r[1]])

    # botsing tussen de deeltjes onderling
    if particleA.collide_detection(particleB):
        dx = particleA.r[0] - particleB.r[0]
        dy = particleA.r[1] - particleB.r[1]
        dist = np.sqrt(dx**2 + dy**2)
        nx = dx / dist
        ny = dy / dist
        
        # Relatieve snelheid
        dvx = particleA.v[0] - particleB.v[0]
        dvy = particleA.v[1] - particleB.v[1]
        dvn = dvx * nx + dvy * ny
        
        if dvn < 0:  # naderen
            # Elastische botsing (gelijke massa)
            particleA.v[0] -= dvn * nx
            particleA.v[1] -= dvn * ny
            particleB.v[0] += dvn * nx
            particleB.v[1] += dvn * ny
            
            # Uit elkaar duwen om overlap te voorkomen
            overlap = particleA.R + particleB.R - dist
            particleA.r[0] += overlap * nx * 0.5
            particleA.r[1] += overlap * ny * 0.5
            particleB.r[0] -= overlap * nx * 0.5
            particleB.r[1] -= overlap * ny * 0.5

    # botsing met de wand voor particleA (met positiecorrectie)
    if particleA.r[0] - particleA.R < -10:
        particleA.r[0] = -10 + particleA.R
        particleA.v[0] *= -1
    elif particleA.r[0] + particleA.R > 10:
        particleA.r[0] = 10 - particleA.R
        particleA.v[0] *= -1
    
    if particleA.r[1] - particleA.R < -10:
        particleA.r[1] = -10 + particleA.R
        particleA.v[1] *= -1
    elif particleA.r[1] + particleA.R > 10:
        particleA.r[1] = 10 - particleA.R
        particleA.v[1] *= -1

    # botsing met de wand voor particleB (met positiecorrectie)
    if particleB.r[0] - particleB.R < -10:
        particleB.r[0] = -10 + particleB.R
        particleB.v[0] *= -1
    elif particleB.r[0] + particleB.R > 10:
        particleB.r[0] = 10 - particleB.R
        particleB.v[0] *= -1
    
    if particleB.r[1] - particleB.R < -10:
        particleB.r[1] = -10 + particleB.R
        particleB.v[1] *= -1
    elif particleB.r[1] + particleB.R > 10:
        particleB.r[1] = 10 - particleB.R
        particleB.v[1] *= -1
    
    return dotA, dotB

# Creeer animatie
ani = FuncAnimation(fig, update, frames=range(num_steps), init_func=init, blit=True, interval=50)

# Voor Jupyter notebook:
from IPython.display import HTML
HTML(ani.to_jshtml())
Loading...
<Figure size 640x480 with 1 Axes>

Terug naar ons vraagstuk... we willen een simulatie waarbij een deeltje met een massa m1m_1 en snelheid v1v_1 op een andere stilstaand deeltje met massa m2m_2 botst. Deeltje twee beweegt naar een muur, botst tegen de muur en beweegt richting deeltje 1 en botst tegen dit deeltje. Hoe vaak vindt deze botsing plaats als functie van de massa verhouding m1m2\frac{m_1}{m_2}?

Daarvoor moeten we even terug naar het botsingsmodel zoals geleerd in Klassieke Mechanica. Bij elastische botsingen is zowel het impulsmoment als de kinetische energie behouden.

∑ip⃗ibefore=∑ip⃗iafter\sum_i \vec{p}_i^{before} = \sum_i \vec{p}_i^{after}
(1)
Ekin,before=Ekin,afterE_{kin, before} = E_{kin, after}
(2)

Voor een botsing met twee deeltjes levert dit een analytische oplossing:

v1′=v1−2m2m1+m2 ⟨v1−v2, x1−x2⟩∥x1−x2∥2 (x1−x2),v2′=v2−2m1m1+m2 ⟨v2−v1, x2−x1⟩∥x2−x1∥2 (x2−x1)\begin{align} \mathbf{v}'_1 &= \mathbf{v}_1-\frac{2 m_2}{m_1+m_2} \ \frac{\langle \mathbf{v}_1-\mathbf{v}_2,\,\mathbf{x}_1-\mathbf{x}_2\rangle}{\|\mathbf{x}_1-\mathbf{x}_2\|^2} \ (\mathbf{x}_1-\mathbf{x}_2), \\ \mathbf{v}'_2 &= \mathbf{v}_2-\frac{2 m_1}{m_1+m_2} \ \frac{\langle \mathbf{v}_2-\mathbf{v}_1,\,\mathbf{x}_2-\mathbf{x}_1\rangle}{\|\mathbf{x}_2-\mathbf{x}_1\|^2} \ (\mathbf{x}_2-\mathbf{x}_1) \end{align}
(3)

we vragen je natuurlijk niet om deze vergelijking zelf te schrijven in Python, maar die vergelijking operationaliseren we hieronder wel.

import numpy as np
import matplotlib.pyplot as plt
from matplotlib.animation import FuncAnimation

class ParticleClass:
    def __init__(self, m, v, r, R):
        self.m = m                  # mass of the particle
        self.v = np.array(v, dtype=float)  # velocity vector
        self.r = np.array(r, dtype=float)  # position vector
        self.R = float(R)  # radius as scalar (not array)

    def update_position(self):
        self.r += self.v * dt 

    def collide_detection(self, other):
        return np.linalg.norm(self.r - other.r) <= (self.R + other.R)

# Simulation parameters
dt = 0.1         # time step
num_steps = 530  # number of time steps
m1 = 1.0
m2 = 1.0

particleA = ParticleClass(m=m1, v=[0, 0], r=[-4.0, 0.0], R=0.45)  
particleB = ParticleClass(m=m2, v=[-1, 0], r=[-2.0, 0.0], R=0.45)  

fig, ax = plt.subplots()
ax.set_xlim(-10, 10)
ax.set_ylim(-10, 10)
ax.set_aspect('equal')
ax.set_title("Particle Animation")
ax.set_xlabel("x")
ax.set_ylabel("y")

dotA, = ax.plot([], [], 'ro', markersize=10)
dotB, = ax.plot([], [], 'bo', markersize=10)

counter = 0
counter_text = ax.text(-9.5, 9, "", fontsize=12, color='black')

# Initialization function for animation
def init():
    dotA.set_data([], [])
    dotB.set_data([], [])
    counter_text.set_text("")
    return dotA, dotB, counter_text

# Update function for each frame
def update(frame):
    global counter
    particleA.update_position()
    particleB.update_position()

    dotA.set_data([particleA.r[0]], [particleA.r[1]])
    dotB.set_data([particleB.r[0]], [particleB.r[1]])

    # collision detection and response between particles
    if particleA.collide_detection(particleB):
        vA, vB, mA, mB, rA, rB = particleA.v, particleB.v, particleA.m, particleB.m, particleA.r, particleB.r
        vA_new = vA - 2 * mB / (mA + mB) * np.dot(vA - vB, rA - rB) / (1e-12+np.linalg.norm(rA - rB))**2 * (rA - rB)
        vB_new = vB - 2 * mA / (mA + mB) * np.dot(vB - vA, rB - rA) / (1e-12+np.linalg.norm(rB - rA))**2 * (rB - rA)
        particleA.v = vA_new
        particleB.v = vB_new
        counter += 1  # particle-particle collision

    # wall collision detection and response (with position correction)
    for particle in [particleA, particleB]:
        # Horizontal walls
        if particle.r[0] - particle.R < -10:
            particle.r[0] = -10 + particle.R
            particle.v[0] = -particle.v[0]
            counter += 1  # wall collision counted
        elif particle.r[0] + particle.R > 10:
            particle.r[0] = 10 - particle.R
            particle.v[0] = -particle.v[0]
            counter += 1  # wall collision counted
        
        # Vertical walls
        if particle.r[1] - particle.R < -10:
            particle.r[1] = -10 + particle.R
            particle.v[1] = -particle.v[1]
            counter += 1  # wall collision counted
        elif particle.r[1] + particle.R > 10:
            particle.r[1] = 10 - particle.R
            particle.v[1] = -particle.v[1]
            counter += 1  # wall collision counted

    counter_text.set_text(f"Collisions: {counter}")
    return dotA, dotB, counter_text

# Create animation
ani = FuncAnimation(fig, update, frames=range(num_steps), init_func=init, blit=True, interval=50)

# For Jupyter notebook:
from IPython.display import HTML
HTML(ani.to_jshtml())
Loading...
<Figure size 640x480 with 1 Axes>
VSC
Introductie thermo simulaties (1h)
VSC
Recap