--- /dev/null
+import sys\r
+import os\r
+import random\r
+import pygame\r
+\r
+scale = 20 #Scale in pixels per angstrom\r
+\r
+mass = {"proton": 10**4, "electron":1} #Masses normalised to electron mass\r
+charge = {"proton":+1, "electron":-1} #Charges normalised to fundamental charge\r
+k = 10**-3 #Coulomb constant\r
+A = 10**-3 #Repulsion constant\r
+m = 6 #Inverse power for the repulsion law (Fr = A/r^m)\r
+\r
+class particle():\r
+ def __init__(self, mass, charge, r):\r
+ """Create a particle"""\r
+ self.mass = mass\r
+ self.charge = charge\r
+ self.r = r\r
+ self.v = [0,0]\r
+ self.f = [0,0]\r
+\r
+ def step(self):\r
+ """Do a single step"""\r
+ #Update velocity\r
+ if self.mass != 0:\r
+ for i in range(0,2):\r
+ self.v[i] += self.f[i] / self.mass\r
+ \r
+ #Update position \r
+ for i in range(0,2):\r
+ self.r[i] += self.v[i]\r
+\r
+ #Reflect off walls\r
+ if self.r[0] <= 0 or self.r[0] >= 640:\r
+ self.v[0] = -self.v[0]\r
+ if self.r[1] <= 0 or self.r[1] >= 480:\r
+ self.v[1] = -self.v[1]\r
+ \r
+ \r
+\r
+ \r
+ \r
+ \r
+ #sys.stdout.write("Position of " + str(self) + " is " + str(self.r) + "\n")\r
+ \r
+\r
+ def calcForce(a,b):\r
+ r = [0,0]\r
+ for i in range(0,2):\r
+ r[i] = b.r[i] - a.r[i]\r
+ if r[0] == 0 and r[1] == 0:\r
+ return [0,0]\r
+ #if r[0] > 10**-10 or r[1] > 10**-10:\r
+ # return [0,0]\r
+ f = [0,0]\r
+ #Coulomb force\r
+\r
+ distance = (r[0]**2 + r[1]**2)**(1/2)\r
+ #sys.stderr.write(str(r[i] * (k * a.charge * b.charge)) + " ")\r
+ #sys.stderr.write("Distance is " + str(distance) + "\n")\r
+ for i in range(0,2):\r
+ #Print out calculation\r
+ #sys.stderr.write("f["+str(i)+"] -= " + str(r[i]) + " * (" + str(k) + " * " + str(a.charge) + " * " + str(b.charge) + " / (" + str(r[0]) + "^2 + " + str(r[1]) + "^2)^(3/2)\n")\r
+ f[i] -= r[i] * (k * a.charge * b.charge / (distance**3))\r
+ #Repulsive force\r
+ #if distance > 10:\r
+ for i in range(0,2):\r
+ f[i] -= r[i] * (A / ((distance)**(m+1)))\r
+ #else:\r
+ # for i in range(0,2):\r
+ # f[] -= r[i] * (A / ((distance-10)**(m+1)))\r
+ return [f[0], f[1]]\r
+\r
+ def draw(self, window):\r
+ \r
+ if self.charge == 1:\r
+ colour = pygame.Color(255,0,0,1)\r
+ size = 4\r
+ elif self.charge == -1:\r
+ colour = pygame.Color(0,0,255,1)\r
+ size = 4\r
+ elif self.charge == 0:\r
+ colour = pygame.Color(255,255,255,1)\r
+ size = 4\r
+ else:\r
+ colour = pygame.Color(0,255,0,1)\r
+ size = 4\r
+\r
+ pygame.draw.circle(window,colour,[int(self.r[0] * scale), int(self.r[1] * scale)],size,0)\r
+ #pygame.draw.line(window,colour,[int((self.r[0] - self.v[0]) / scale), int((self.r[1] - self.v[1]) / scale)], [int(self.r[0] / scale), int(self.r[1] / scale)])\r
+ \r
+ \r
+\r
+ \r
+class universe():\r
+ def __init__(self):\r
+ self.particles = []\r
+ for i in range(0,0):\r
+ self.particles.append(particle(mass["proton"],charge["proton"], universe.getRandomPosition()))\r
+ self.particles.append(particle(mass["proton"] + 2*mass["electron"],charge["electron"],universe.getRandomPosition()))\r
+\r
+ def step(self):\r
+ for a in self.particles:\r
+ a.f = [0,0]\r
+ for b in self.particles:\r
+ if b != a:\r
+ fAB = particle.calcForce(a,b)\r
+ for i in range(0,2):\r
+ a.f[i] += fAB[i]\r
+\r
+ a.step()\r
+\r
+ def getRandomPosition():\r
+ return [float(random.randint(0,640) / scale), float(random.randint(0,480) / scale)]\r
+\r
+ def draw(self, window):\r
+ for a in self.particles:\r
+ a.draw(window)\r
+\r
+if __name__ == "__main__":\r
+ #create the screen\r
+ window = pygame.display.set_mode((640, 480))\r
+ theUniverse = universe()\r
+\r
+ run = True\r
+ while run:\r
+ #Clear screen\r
+ \r
+ \r
+ #Step\r
+ #sys.stdout.write("Stepping...\n")\r
+ theUniverse.step()\r
+ #Draw\r
+ window.fill(pygame.Color(0,0,0,0))\r
+ theUniverse.draw(window)\r
+ #Update draw events\r
+ pygame.display.flip()\r
+ #Wait a bit\r
+ #pygame.time.wait(1)\r
+ #Input handling\r
+ for event in pygame.event.get(): \r
+ if event.type == pygame.QUIT: \r
+ pygame.quit()\r
+ run = False\r
+ elif event.type == pygame.MOUSEBUTTONDOWN:\r
+ mousestate = pygame.mouse.get_pressed()\r
+ mousepos = pygame.mouse.get_pos()\r
+ if mousestate[0] and mousestate[2] == False:\r
+ theUniverse.particles.append(particle(mass["proton"],charge["proton"], [mousepos[0]/scale,mousepos[1]/scale]))\r
+ elif mousestate[2] and mousestate[0] == False:\r
+ theUniverse.particles.append(particle(mass["electron"],charge["electron"], [mousepos[0]/scale,mousepos[1]/scale]))\r
+ break\r
+ else: \r
+ #print(event)\r
+ pass\r
+ \r