*

Offline spoon

  • *
  • Posts: 1134
  • Foxy wins
    • View Profile
cool fun programming stuff
« on: April 27, 2015, 03:20:30 AM »
Code: [Select]
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
« Last Edit: April 27, 2015, 03:28:19 AM by spoon »
inb4 Blanko spoons a literally pizza

*

Offline xasop

  • Administrator
  • *****
  • Posts: 9776
  • Professional computer somebody
    • View Profile
Re: cool fun programming stuff
« Reply #1 on: April 27, 2015, 03:26:42 AM »
Wrap your code in [code] tags, and that won't happen.
when you try to mock anyone while also running the flat earth society. Lol

*

Offline spoon

  • *
  • Posts: 1134
  • Foxy wins
    • View Profile
Re: cool fun programming stuff
« Reply #2 on: May 04, 2015, 05:30:35 AM »
Code: [Select]
from __future__ import division
from visual import *
from random import *
from visual.graph import *

#define constants and other parameters
number = 9 #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
KE = 0 #kinetic energy of the system
t = 0 #time
f1 = gcurve(color=color.cyan) # a graphics curve

#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)

#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))
#charges[0].charge = *electrons*q

#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

#click to start moving charges
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:
        t = t + 1 #increment time for kinetiv energy graph
        if scene.mouse.clicked: #click to break from loop and draw lines between charges
            break
        for a in arange(number): #This for takes a charge carrier and compares it to...
            for b in arange(number): #... 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
                    #Is this the best way to update velocity?
                    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
       
                    #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
                       
    #displays force arrow with each update. Change it so old arrows disappear. also traces position with dots.
        #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 / 50, color=color.yellow, opacity = .1))

        #update position of charges 
        for a in arange(number):
            rate(300)
            charges[a].pos = charges[a].pos + charges[a].velocity
            for i in arange(number):
                KE = KE + mag(charges[i].velocity)**2   
            f1.plot(pos=(t,KE,0))
            KE = 0
           
#This loop draws lines between each charge carrier after the second click.
for a in arange(number):
    for b in arange(number):
        if a < b:
            curve(pos = [charges[a].pos, charges[b].pos])

#This loop creates a convex solid out of the charges' positions and makes charges transparent.   
scene.waitfor('keydown')
a = convex(color = color.white)
for i in arange(number):
    a.append(pos = charges[i].pos)
    charges[i].opacity = 0
   
#make other objects transparent
item1.opacity = 0
xaxis.pos = [zero, zero]
yaxis.pos = [zero, zero]
zaxis.pos = [zero, zero]

Basically finished product, with added functionality. Now, after viewing the beautiful symmetry created by the lines between each vertex, you can view the solid created by those vertices. I think it's gorgeous stuff.

With a high number of charge carriers, the solid created begins to resemble a geodisic sphere.

You also get a graph that shows the kinetic energy of the system with respect to time. woo
inb4 Blanko spoons a literally pizza

*

Offline Tom Bishop

  • Zetetic Council Member
  • **
  • Posts: 10637
  • Flat Earth Believer
    • View Profile
Re: cool fun programming stuff
« Reply #3 on: May 09, 2015, 04:06:13 PM »
How do we get the visual module?

PS C:\Users\tbishop\python> python .\spoon.py
Traceback (most recent call last):
  File ".\spoon.py", line 2, in <module>
    from visual import *
ImportError: No module named visual
PS C:\Users\tbishop\python>

*

Offline Tom Bishop

  • Zetetic Council Member
  • **
  • Posts: 10637
  • Flat Earth Believer
    • View Profile
Re: cool fun programming stuff
« Reply #4 on: May 09, 2015, 04:28:02 PM »
Okay, so I downloaded vpython for win32, and pasted the code into the vpython editor, but can't seem to run it. I selected my base python 27 wpython.exe when the program ran. This is what happens when I press F5:

« Last Edit: May 09, 2015, 06:28:00 PM by Tom Bishop »

*

Offline Pete Svarrior

  • e
  • Planar Moderator
  • *****
  • Posts: 16073
  • (◕˽ ◕ ✿)
    • View Profile
Re: cool fun programming stuff
« Reply #5 on: May 09, 2015, 05:49:05 PM »
This is a bit of a random guess, but: make sure you don't have any .py files in the directory where pythonw.exe is located
Read the FAQ before asking your question - chances are we already addressed it.
Follow the Flat Earth Society on Twitter and Facebook!

If we are not speculating then we must assume

*

Offline spoon

  • *
  • Posts: 1134
  • Foxy wins
    • View Profile
Re: cool fun programming stuff
« Reply #6 on: May 09, 2015, 09:11:21 PM »
Use Vidle, not Idle.
inb4 Blanko spoons a literally pizza