from __future__ import division
from visual import *
from random import *
#define constants and other parameters
number = 8 #number of charges
q = 1.6e-19 #1 unit of charge (charge of an electron)
k = 9e9 # k (1/4piEnaught)
rad = 100 #radius of spherical region
shrink = 1.75 #used to make sure all charges begin within the sphere
zero = vector(0,0,0) #used for declaration purposes
electrons = 7e11*(1/number) #puts a lot of charge on each charge carrier
d = .5 #damping effect the inside of the sphere has on the charges (not used)
#generate translucent spherical region
item1 = sphere(pos = (0,0,0), radius = rad, color = color.white, opacity = .3)
#generate axes for clarity in charge position and depth within the sphere
xaxis = curve(pos = [(-1.5*rad,0,0),(1.5*rad,0,0)], color = color.blue, radius = rad/100)
yaxis = curve(pos = [(0,-1.5*rad,0),(0,1.5*rad,0)], color = color.blue, radius = rad/100)
zaxis = curve(pos = [(0,0,-1.5*rad),(0,0,1.5*rad)], color = color.blue, radius = rad/200)
#label axes
text(text = 'x', pos = (1.6*rad,-rad/30,rad/50), height = rad/15, depth= -rad/25, color=color.white)
text(text = 'y', pos = (-rad/20,1.6*rad,rad/50), height = rad/10, depth= -rad/25, color=color.white)
text(text = 'z', pos = (-rad/20,-rad/20,1.6*rad), height = rad/10, depth= -rad/25, color=color.white)
#Create a list and fill it with charge carriers in random locations in the spherical region.
charges = []
for i in arange(number):
charges.append(sphere(pos = ((randrange(-rad*100,rad*100)/(shrink*100)),(randrange(-rad*100,rad*100)/(shrink*100)),(randrange(-rad*100,rad*100)/(shrink*100))), radius = rad/33, charge = electrons*q, force = zero, velocity = zero, color = color.red))
#This list will be used to make arrows and spheres and what have you.
objects = []
#This function is used to find the distance between points, which is needed to calculate electric field. a and b are charge carriers from 'charges'.
def getDistance(a,b):
distance = mag(charges[a].pos-charges[b].pos)
return distance
#This function calculates net force. a = which charge is having calculations done on it. b = the charge whose field is contributing a force to a. c = distance from getDistance function.
def getNetForce(a,b,c):
direction = charges[a].pos - charges[b].pos
charges[a].force = (((k*charges[b].charge)/(c**2))*norm(direction))
return charges[a].force
scene.waitfor('click')
#This chunk updates force, velocity, and position for each charge carrier.
for i in arange(number):
while mag(charges[i].pos) < 2 * rad:
for a in arange(number): #This for takes a charge carrier and compares it to...
for b in arange(number): #... to a charge carrier in this loop.
if a < b: #This if statement gets rid of overlap so the same calculation is not made twice.
distance = getDistance(a, b) #get distance between charge carriers
charges[a].force = getNetForce(a, b, distance) #get force on a due to b...
charges[b].force = getNetForce(b, a, distance) #... and on b due to a
charges[a].velocity = charges[a].velocity + charges[a].force #update...
charges[b].velocity = charges[b].velocity + charges[b].force #...velocity
sizeA = mag(charges[a].pos) #distance of first charge carrier from center...
sizeB = mag(charges[b].pos) #... and the second
#Charges near boundary, take two: Newtonian approach. This makes the charges "bounce off the inside of the wall by using the Rodrigues rotation formula.
if sizeA >= rad:
charges[a].pos = .999 * rad * norm(charges[a].pos) #moves charge carrier just inside sphere to avoid the charge getting "stuck"
posUnit = norm(charges[a].pos) #unit vector in direction of charge position
projection = proj(charges[a].velocity, charges[a].pos) #projection of velocity vector onto position vector
charges[a].velocity = -.9 *((charges[a].velocity * cos(pi)) + (cross(posUnit, charges[a].velocity) * sin(pi)) + (projection * (1 - cos(pi)))) #Rodrigues' rotation formula
if sizeB >= rad:
charges[b].pos = .999 * rad * norm(charges[b].pos) #moves charge carrier just inside sphere to avoid the charge getting "stuck"
posUnit = norm(charges[b].pos) #unit vector in direction of charge position
projection = proj(charges[b].velocity, charges[b].pos) #projection of velocity vector onto position vector
charges[b].velocity = -.9 *((charges[b].velocity * cos(pi)) + (cross(posUnit, charges[b].velocity) * sin(pi)) + (projection * (1 - cos(pi)))) #Rodrigues' rotation formula
## Clunky clode that deals with charges near boundary. Doesn't work well at all, sometimes glitchy. physics inconsistent
##
## if sizeA > rad:
## radialDirectionA = norm(charges[a].pos)
## charges[a].pos = radialDirectionA * rad
## crossA = cross(charges[a].pos, charges[a].force)
## newForceDirA = cross(crossA, charges[a].pos)
## angleA = diff_angle(charges[a].pos,charges[a].force)
## #fixed so it doesn't get glitchy, but doesn't redistribute charges evenly on surface.
## charges[a].force = mag(charges[a].force) * sin(2 * (pi - angleA)) * norm(newForceDirA) * d
## charges[a].velocity = charges[a].force
## if sizeB > rad:
## radialDirectionB = norm(charges[b].pos)
## charges[b].pos = radialDirectionB * rad
## crossB = cross(charges[a].pos, charges[a].force)
## newForceDirB = cross(crossB, charges[a].pos)
## angleB = diff_angle(charges[a].pos,charges[a].force)
## #fixed so it doesn't get glitchy, but doesn't redistribute charges evenly on surface.
## charges[b].force = mag(charges[b].force) * sin(2 * (pi - angleB)) * norm(newForceDirB) * d
## charges[b].velocity = charges[b].force
#Optional section that displays force arrow with each update. Change it so old arrows disappear. also traces position with translucent balls.
#for i in arange(number):
#objects.append(arrow(pos=charges[i].pos, axis=charges[i].force, color=color.yellow)
#objects.append(sphere(pos=charges[i].pos, radius = rad / 33, color=color.red, opacity = .1))
for a in arange(number):
rate(100)
charges[a].pos = charges[a].pos + charges[a].velocity #update position of charges...
#print mag(charges[0].force)
#print mag(charges[0].velocity)
Above is a program I am writing in vpython to study geometry by simulating the electric force! Placing
n like charge carriers in a hollow, neutral, insulating sphere, one would expect the to disperse among the surface of the sphere in such a way that equilibrium is reached. My supposition was that for charge carriers > 3, equilibrium would only be reached if the number of charges corresponded to the number of vertices in a
platonic polyhedral. This supposition was based on intuition and a limited understanding of geometry in 3 space.
I was pleasantly surprised. For every number of charge carriers I have tried, equilibrium is reached, sometimes in fascinating ways. My favorite is 8. Instead of forming a cube, the charges form a polyhedral with two parallel square bases, rotated such that it has 8 other faces that are equilateral triangles. Obviously, this is the orientation that results in a net force on each charge of zero, but it also seems correct to say that the same would be true for a cube. Each vertex on a cube is the same average distance from any other vertex.
This is still a work in progress, I need to streamline it and clean it up, but it is functional. Over the next few days, I am going to do some more research into this sort of geometry and see what I can find. In the mean time, feel free to copy/paste this into the most recent version of vpython and give it a whirl! Other cool things to do would be changing the charge on one or more charge carriers. If you have any ideas or insights, let me know.
Also, this can be a general thread for fun programming adventures.
Also, lol at everything being bold from calling "b" from a list "charges[]"
*thanks parsifal