`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