Simulation anziehender und abstoßender Felder in Python 3 (Software)

Martin Vogel ⌂ @, Dortmund / Bochum, Sun, 01.01.2017, 17:54 (vor 2697 Tagen)

Eine Simulation ohne allzu konkreten physikalischen Bezug zeigt die Wirkung zweier entgegengesetzt wirkender Felder auf eine Anzahl Partikel. Ähnlich der Gravitation werden die Partikel mit einer Kraft angezogen, die sich umgekehrt quadratisch zum Abstand der einzelnen Teilchen verhält. Gleichzeitig werden sie von einer Kraft mit stärkerer Wirkung aber geringerer effektiver Reichweite abgestoßen. Ganz entfernt erinnert das Verhalten an die Wechselwirkung der Elementarkräfte im Atomkern.

Sinn der Übung ist weniger die Darstellung eines realen physikalischen Vorgangs (insofern muss ich alle enttäuschen, die das als Last-Minute-Idee für die aktuelle Programmieraufgabe aufgreifen wollen) als vielmehr, ein Programmierbeispiel für den Einfluss eines Reglerpaneels auf eine laufende Simulation zu sein. Diese wird zu Beginn des Programms einmal angestoßen und läuft dann kontinuierlich weiter, während alle möglichen Parameter bis hin zur Partikelanzahl lustig geändert werden können. Gelegentlich läuft die Simulation total aus dem Ruder, was sich in der Regel durch entschlossenes Hochziehen des Dämpfungsreglers schnell beheben lässt.

Wer kann mir das kuriose Verhalten der Partikel erklären, wenn man den Regler für die Abstoßung ganz oder fast ganz auf null herunterzieht?

[image]
[image]

Download: [image]ZIP-Archiv_2ABK27ET1.zip

  1. #!/usr/bin/python3
  2.  
  3. # Ein Haufen Partikel, der sich auf große Entfernung anzieht
  4. # und auf geringe Entfernung abstößt
  5.  
  6. #
  7. # Version vom 1. Januar 2017
  8. # Autor: Martin Vogel, martinvogel.de
  9. # Lizenz: cc-by-3.0 https://creativecommons.org/licenses/by/3.0/de/
  10. #
  11.  
  12. from tkinter import *
  13. from random import randint, choice, random
  14. from time import perf_counter
  15. from math import hypot
  16.  
  17.  
  18. class Einstellungen:
  19. """Einige programmweit gültige Einstellungen, die meisten sind über
  20. # Schieberegler beeinflussbar."""
  21.  
  22. def __init__(self):
  23. # Die Anfangsgröße der Canvas
  24. self.Höhe = 600
  25. self.Breite = 800
  26.  
  27. # Die Anfangszahl der Partikel
  28. self.Anzahl = 200
  29.  
  30. # Die Elastizität der Wand:
  31. # 1: Impuls wird voll zurückgegeben, 0: Partikel "klebt" an der Wand.
  32. self.Wandelastizität = 0.5
  33.  
  34. # Anteil, um den die Geschwindigkeit pro Rechenschritt vermindert wird.
  35. self.Dämpfung = 0.01
  36. # Anmerkung: da die Zahl der Rechenschritte pro Sekunde von der Anzahl
  37. # der Partikel abhängt (der Rechenaufwand wächst quadratisch), wirkt
  38. # sich die Dämpfung bei geringer Partikelzahl viel stärker aus als bei
  39. # großer Partikelzahl. Wer mag, kann das gerne mal verbessern.
  40.  
  41. # Die Stärke von Anziehung und Abstoßung.
  42. self.Anziehungsfaktor = 9
  43. self.Abstoßungsfaktor = 18
  44.  
  45. # Bezugswert für Zeitmessungen:
  46. self.Zeit = perf_counter()
  47.  
  48. def setze_dämpfung(self, event):
  49. # Lies den Schieberwert!
  50. self.Dämpfung = SD.get()
  51.  
  52. def setze_anzahl(self, event=None):
  53. # Wo steht der Schieber?
  54. self.Anzahl = SA.get()
  55. while len(Partikel) < self.Anzahl:
  56. # Brauchen wir noch Partikel? Neuen erzeugen!
  57. C_Partikel()
  58. while len(Partikel) > self.Anzahl:
  59. # Haben wir zu viele Partikel? Irgendeins löschen!
  60. choice(Partikel).lösche()
  61.  
  62. def setze_anziehungsfaktor(self, event):
  63. # Hole den Schieberwert!
  64. self.Anziehungsfaktor = SAn.get()
  65. # Die Grafik im Kontrollpaneel aktualisieren
  66. C.after_idle(ZeigeVerlauf)
  67.  
  68. def setze_abstoßungsfaktor(self, event):
  69. # Welcher Schieberwert wurde eingestellt?
  70. self.Abstoßungsfaktor = SAb.get()
  71. # Die Grafik im Kontrollpaneel aktualisieren
  72. C.after_idle(ZeigeVerlauf)
  73.  
  74. def setze_wandelastizität(self, event):
  75. # Lies den Schieber, setze das Attribut:
  76. self.Wandelastizität = SW.get()
  77.  
  78. def setze_canvasgröße(self, event):
  79. # Falls die Größe der Canvas verändert wird, neue Abmessungen auslesen
  80. # und merken:
  81. self.Breite = C.winfo_width()
  82. self.Höhe = C.winfo_height()
  83.  
  84.  
  85. # Die einzige Instanz dieser Klasse:
  86. E = Einstellungen()
  87.  
  88. # Alle Partikel werden in dieser Liste verwaltet.
  89. Partikel = []
  90.  
  91.  
  92. def zufallsfarben():
  93. "Erzeugt einen nicht zu dunklen RGB-Farbwert und eine dunklere Variante"
  94. # Ein Farbpaar für die Partikel: Heller Kern, dunkler Rand
  95. r, g, b = randint(64, 255), randint(64, 255), randint(64, 255)
  96. return "#%02x%02x%02x" % (r, g, b), "#%02x%02x%02x" % (r//2, g//2, b//2)
  97.  
  98.  
  99. class C_Partikel:
  100. "Massepunkt mit Ort und Geschwindigkeit."
  101. def __init__(self, x=None, y=None):
  102. # Ort
  103. if x is None:
  104. x = randint(0, E.Breite)
  105. if y is None:
  106. y = randint(0, E.Höhe)
  107. self.x = x
  108. self.y = y
  109. # Geschwindigkeit
  110. self.vx = 0
  111. self.vy = 0
  112. # Farbe
  113. hell, dunkel = zufallsfarben()
  114. # Canvasobjekt
  115. self.ID = C.create_oval(self.x-3, self.y-3,
  116. self.x+3, self.y+3,
  117. fill=hell,
  118. outline=dunkel)
  119. # Eintragen in die Liste aller Partikel
  120. Partikel.append(self)
  121. # Der Punkt kann mit der Maus verschoben werden.
  122. C.tag_bind(self.ID, "<B1-Motion>", self.bewegezu)
  123.  
  124. def lösche(self):
  125. "Punkt wird entfernt"
  126. # Lösche Punkt aus der Partikelliste
  127. Partikel.remove(self)
  128. # Lösche Canvas-Objekt
  129. C.delete(self.ID)
  130.  
  131. def bewegezu(self, event):
  132. "Punkt wird verschoben"
  133. C.move(self.ID, event.x-self.x, event.y-self.y)
  134. self.x, self.y = event.x, event.y
  135.  
  136.  
  137. def Animationsschleife(event=None):
  138. """Einsammeln aller Kräfte, die auf jedes Partikel wirken und anschließende
  139. Berechnung der Kraftwirkung"""
  140.  
  141. # Berechnung der seit dem letzten Aufruf verstrichenen Zeit delta_t
  142. jetzt = perf_counter()
  143. delta_t = jetzt-E.Zeit
  144. E.Zeit = jetzt
  145.  
  146. Anz, Abs = E.Anziehungsfaktor, E.Abstoßungsfaktor
  147. Dmp = E.Dämpfung
  148. Wez = E.Wandelastizität
  149.  
  150. # Kopie der Partikelliste erzeugen, damit es keine Verfälschungen durch
  151. # Neuberechnungen gibt.
  152. P = Partikel[:]
  153. # Für alle Partikel:
  154. for i, Pi in enumerate(P):
  155. # Aktuelle Koordinaten
  156. x, y = Pi.x, Pi.y
  157. # Beschleunigungssumme in x- und y-Richtung
  158. ax = ay = 0
  159. # Einwirkungen aller anderen Partikel einsammeln:
  160. for Pj in P:
  161. # Abstände in x- und y-Richung
  162. dx = Pj.x-Pi.x
  163. dy = Pj.y-Pi.y
  164. # Abstand diagonal
  165. d = hypot(dx, dy)
  166. if d:
  167. # Falls Pj nicht dasselbe Partikel wie Pi ist:
  168. # Die Beschleunigungsgleichung:
  169. # Die Anziehung lässt mit dem Quadrat der Entfernung nach
  170. # (ähnlich der Gravitation).
  171. # Die Abstoßung soll hier mit der dritten Potenz der
  172. # Entfernung schwächer werden (ohne physikalische Grundlage).
  173. a = (Anz/d)**2 - (Abs/d)**3
  174. # Beschleunigungskomponenten aufaddieren
  175. ax += dx*a
  176. ay += dy*a
  177. # Nachdem die resultierende Beschleunigung aufaddiert wurde, kann
  178. # die Geschwindigkeitsänderung berechnet werden:
  179. Partikel[i].vx = (ax*delta_t+Pi.vx)*(1-Dmp)
  180. Partikel[i].vy = (ay*delta_t+Pi.vy)*(1-Dmp)
  181. # Daraus ergibt sich der neue Ort des Partikels.
  182. x += Partikel[i].vx*delta_t
  183. y += Partikel[i].vy*delta_t
  184. # Bei den folgenden Kollisionen wird das Partikel um eine zufällige
  185. # Entfernung von der Wand abgesetzt, um kuriose Effekte zu vermeiden,
  186. # die entstehen, wenn ein Haufen benachbarter Partikel dieselben
  187. # x- oder y- Werte aufweisen.
  188.  
  189. if x < 0:
  190. # Kollision mit linker Wand?
  191. x = random()/10
  192. Partikel[i].vx *= -Wez
  193. Partikel[i].vy *= Wez
  194.  
  195. elif x > E.Breite:
  196. # Kollision mit rechter Wand?
  197. x = E.Breite-random()/10
  198. Partikel[i].vx *= -Wez
  199. Partikel[i].vy *= Wez
  200.  
  201. if y < 0:
  202. # Kollision mit oberer Wand?
  203. y = random()/10
  204. Partikel[i].vx *= Wez
  205. Partikel[i].vy *= -Wez
  206.  
  207. elif y > E.Höhe:
  208. # Kollision mit unterer Wand?
  209. y = E.Höhe-random()/10
  210. Partikel[i].vx *= Wez
  211. Partikel[i].vy *= -Wez
  212.  
  213. # Canvasobjekt des Partikels auf neue Koordinaten setzen
  214. C.move(Partikel[i].ID, x-Partikel[i].x, y-Partikel[i].y)
  215. Partikel[i].x, Partikel[i].y = x, y
  216.  
  217. # Sobald die CPU wieder Luft geholt hat, geht’s von vorne los.
  218. C.after_idle(Animationsschleife)
  219.  
  220. # Das Hauptfenster mit der Canvas
  221. T = Tk()
  222. T.title("Wechselwirkung nah und fern")
  223.  
  224. # Die Zeichenfläche. Sie passt sich Größenänderungen des Fensters an.
  225. C = Canvas(width=E.Breite, height=E.Höhe, bg="black", highlightthickness=0)
  226. C.pack(expand=True, fill="both")
  227. C.bind("<Configure>", E.setze_canvasgröße)
  228.  
  229. # Das Kontrollfenster
  230. KF = Toplevel()
  231. KF.title("Reglerpult")
  232. # Das Reglerpult
  233. SLF = LabelFrame(KF, text="Parameter", padx=5, pady=5)
  234. SLF.grid(row=0, column=0, columnspan=2, sticky="wens", padx=5, pady=5)
  235.  
  236. # Zeile 1 im Labelframe (die Regler) soll dessen Größenanpassung übernehmen,
  237. # außerdem soll sie mindestens 200 Pixel hoch sein, damit die
  238. # Skalenzahlen gut zu erkennen sind.
  239. KF.rowconfigure(0, weight=1)
  240. SLF.rowconfigure(1, weight=1, minsize=200)
  241. # Die 5 Spalten sollen sich gleichmäßig verteilen können.
  242. SLF.columnconfigure(0, weight=1)
  243. SLF.columnconfigure(1, weight=1)
  244. SLF.columnconfigure(2, weight=1)
  245. SLF.columnconfigure(3, weight=1)
  246. SLF.columnconfigure(4, weight=1)
  247.  
  248. # Schieberegler für die Anzahl der Partikel
  249. Label(SLF, text="Anzahl").grid(row=0, column=0)
  250. SA = Scale(SLF, from_=500, to=0, width=20, orient="vertical",
  251. tickinterval=100, resolution=1, command=E.setze_anzahl)
  252. SA.set(E.Anzahl) # Anfangswert des Reglers
  253. SA.grid(row=1, column=0, sticky="nsew")
  254.  
  255. # Schieberegler für den Anziehungsfaktor
  256. Label(SLF, text="Anziehung").grid(row=0, column=1)
  257. SAn = Scale(SLF, from_=50, to=0, width=20, orient="vertical",
  258. tickinterval=10, resolution=0.1, command=E.setze_anziehungsfaktor)
  259. SAn.set(E.Anziehungsfaktor) # Anfangswert des Reglers
  260. SAn.grid(row=1, column=1, sticky="nsew")
  261.  
  262. # Schieberegler für den Abstoßungsfaktor
  263. Label(SLF, text="Abstoßung").grid(row=0, column=2)
  264. SAb = Scale(SLF, from_=50, to=0, width=20, orient="vertical",
  265. tickinterval=10, resolution=0.1, command=E.setze_abstoßungsfaktor)
  266. SAb.set(E.Abstoßungsfaktor) # Anfangswert des Reglers
  267. SAb.grid(row=1, column=2, sticky="nsew")
  268.  
  269. # Schieberegler für die Wand-Elastizität
  270. Label(SLF, text="Wand-nelastizität").grid(row=0, column=3)
  271. SW = Scale(SLF, from_=1, to=0, width=20, orient="vertical",
  272. tickinterval=0.2, resolution=0.01, command=E.setze_wandelastizität)
  273. SW.set(E.Wandelastizität) # Anfangswert des Reglers
  274. SW.grid(row=1, column=3, sticky="nsew")
  275.  
  276. # Schieberegler für die Dämpfung
  277. Label(SLF, text="Dämpfung").grid(row=0, column=4)
  278. SD = Scale(SLF, from_=1, to=0, width=20, orient="vertical",
  279. tickinterval=0.2, resolution=0.001, command=E.setze_dämpfung)
  280. SD.set(E.Dämpfung) # Anfangswert des Reglers
  281. SD.grid(row=1, column=4, sticky="nsew")
  282.  
  283.  
  284. def ZeigeVerlauf(event=None):
  285. """In einer Canvasgrafik unterhalb der Regler wird der Wirkungsradius
  286. und die Stärke von Anziehung und Abstoßung (innerhalb gewisser Grenzen)
  287. als Diagramm angezeigt."""
  288. # Einige Werte zur Darstellung ausrechnen und zwischenspeichern
  289. Werte = [0]
  290. for d in range(1, 600):
  291. Anziehung = (E.Anziehungsfaktor/d)**2
  292. Abstoßung = (E.Abstoßungsfaktor/d)**3
  293. Werte.append(Anziehung-Abstoßung)
  294.  
  295. # Extrema?
  296. Wmin, Wmax = min(Werte)-0.001, max(Werte)+0.001
  297. # Der betragsmäßig kleinere der beiden Werte bestimmt den Maßstab
  298. Wmin, Wmax = max(Wmin, -2*Wmax), min(Wmax, -2*Wmin)
  299. # … noch etwas Randabstand
  300. Wmax += 0.05*(Wmax-Wmin)
  301. Wmin -= 0.05*(Wmax-Wmin)
  302.  
  303. # Lage der Nulllinie in y-Richtung
  304. Null = 100*(Wmax)/(Wmax-Wmin)
  305.  
  306. # Punkte des Funktionsgraphen
  307. Punkte = []
  308.  
  309. # Skalierung der Werte auf den Darstellungsbereich
  310. for x in range(1, 600):
  311. Punkte.append(x)
  312. y = min(101, max(-1, 100*(Wmax-Werte[x])/(Wmax-Wmin)))
  313. Punkte.append(y)
  314.  
  315. # Radius der Grenze zwischen Anziehung und Abstoßung
  316. try:
  317. r0 = E.Abstoßungsfaktor**3/E.Anziehungsfaktor**2
  318. except:
  319. r0 = float("inf")
  320.  
  321. # Radius wird für die Darstellung auf Anzeigebereich begrenzt
  322. x0 = min(600, r0)
  323.  
  324. # Grafik neu aufbauen
  325. SC.delete("all")
  326. # Nullinie
  327. SC.create_line(0, Null, 600, Null, fill="#008000")
  328. # An-Ab-Grenze
  329. SC.create_line(x0, 0, x0, 100, fill="#008000")
  330. # untenliegende Beschriftung
  331. SC.create_text(x0+2, 2, text="Anziehung", anchor="nw", fill="#40C000")
  332. SC.create_text(x0-2, 2, text="Abstoßung", anchor="ne", fill="#40C000")
  333. # Funktionsgraph
  334. SC.create_line(Punkte, fill="#E0E000")
  335. # Beschleunigungsgleichung
  336. SC.create_text(598, 98, text="a = (%.1f/x)² - (%.1f/x)³" % (
  337. E.Anziehungsfaktor, E.Abstoßungsfaktor),
  338. anchor="se", fill="#E0E000")
  339. # Radius als gestrichelter Kreis mit symbolischem Partikel im Zentrum
  340. SC.create_oval(-3, Null-3, 3, Null+3, fill="red", outline="red")
  341. SC.create_oval(-x0, Null-x0, x0, Null+x0, outline="red", dash=(3, ))
  342. SC.create_text(x0/2, Null-1, text="r0 = %.1f" % r0, fill="red", anchor="s")
  343.  
  344. # Die Canvas im Kontrollfenster
  345. SC = Canvas(KF, width=600, height=100, bg="black", highlightthickness=0)
  346. SC.grid(row=1, column=0, padx=5, pady=5)
  347.  
  348. # Anzeige des Diagramms
  349. ZeigeVerlauf()
  350.  
  351. # Canvas anzeigen
  352. C.update()
  353.  
  354. # Erstellung der Anfangsbestückung an Partikeln
  355. E.setze_anzahl()
  356.  
  357. # Animation starten
  358. Animationsschleife()
  359.  
  360. # Kontrollfenster in den Vordergrund holen
  361. KF.lift()
  362.  
  363. # Lift off!
  364. T.mainloop()

--
Dipl.-Ing. Martin Vogel
Leiter des Bauforums

Bücher:
CAD mit BricsCAD
Bauinformatik mit Python

RSS-Feed dieser Diskussion
powered by my little forum