Created
January 18, 2014 18:51
-
-
Save orzelc/8494567 to your computer and use it in GitHub Desktop.
VPython mock-up of charged tape problem as a mass-on-a-spring, because the equations are easier to deal with.
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
from visual import * | |
# Constants and other definitions | |
e=1.6e-19 | |
k=9e9 | |
g=9.80 | |
length=0.21 | |
height=0.3 | |
sep=0.1524 | |
mass=0.0003 | |
charge=2.0e-8 | |
scale=1 | |
spring=100 | |
drag=0.001 | |
# Pivot points for the pendulums | |
pivot1=sphere(pos=vector(sep/2,height,0), radius=0.02, color=color.white) | |
pivot2=sphere(pos=vector(-sep/2,height,0), radius=0.02, color=color.white) | |
# The suspended point charges that stand in for charged tapes | |
mass1=sphere(pos=vector(sep/2, height-length,0), radius=0.02, color=color.green) | |
mass1.m=mass | |
mass1.q=charge | |
mass1.p=vector(0,0,0) | |
mass1.pos.y=mass1.pos.y-mass1.m*g/spring #stretch to account for the mass | |
mass2=sphere(pos=vector(-sep/2,height-length,0), radius=0.02, color=color.yellow) | |
mass2.m=mass | |
mass2.q=charge | |
mass2.p=vector(0,0,0) | |
mass2.pos.y=mass2.pos.y-mass2.m*g/spring #stretch to account for the mass | |
# Strings for the pendulum, to make it look pretty | |
string1=cylinder(pos=pivot1.pos, axis=mass1.pos-pivot1.pos, radius=0.005, color=color.white) | |
string2=cylinder(pos=pivot2.pos, axis=mass2.pos-pivot2.pos, radius=0.005, color=color.white) | |
# Gravitational forces | |
fgrav1=mass1.m*g*vector(0,-1,0) | |
fgrav2=mass2.m*g*vector(0,-1,0) | |
# Treating the strings as springs to get the right tension force | |
stretch1=string1.axis.mag-length | |
stretch2=string2.axis.mag-length | |
ft1=-(spring*stretch1)*(string1.axis/string1.axis.mag) | |
ft2=-(spring*stretch2)*(string2.axis/string2.axis.mag) | |
# Electrostatic force between the charged balls | |
r=mass1.pos-mass2.pos | |
rhat=r/r.mag | |
felec1=-(k*mass1.q*mass2.q/r.mag2)*rhat | |
felec2=-felec1 | |
# Arrows showing the force | |
fgarr1=arrow(pos=mass1.pos, axis=scale*fgrav1, color=color.blue) | |
fgarr2=arrow(pos=mass2.pos, axis=scale*fgrav2, color=color.blue) | |
ftarr1=arrow(pos=mass1.pos, axis=scale*ft1, color=color.magenta) | |
ftarr2=arrow(pos=mass2.pos, axis=scale*ft2, color=color.magenta) | |
felarr1=arrow(pos=mass1.pos, axis=scale*felec1, color=color.cyan) | |
felarr2=arrow(pos=mass2.pos, axis=scale*felec2, color=color.cyan) | |
# Loop parameters | |
tmax=5 | |
t=0 | |
dt=0.001 | |
while t<tmax: | |
rate(1000) | |
#print(t) | |
#Recalculate forces | |
stretch1=string1.axis.mag-length | |
stretch2=string2.axis.mag-length | |
ft1=-(spring*stretch1)*(string1.axis/string1.axis.mag) | |
ft2=-(spring*stretch2)*(string2.axis/string2.axis.mag) | |
r=mass1.pos-mass2.pos | |
rhat=r/r.mag | |
felec1=-(k*mass1.q*mass2.q/r.mag2)*rhat | |
felec2=-felec1 | |
#Drag force to make the system settle near equlibirum point | |
fdrag1=-drag*mass1.p/mass1.m | |
fdrag2=-drag*mass2.p/mass2.m | |
fnet1=fgrav1+ft1+felec1+fdrag1 | |
fnet2=fgrav2+ft2+felec2+fdrag2 | |
# Update momentum and position of the charged balls | |
mass1.p=mass1.p+fnet1*dt | |
mass1.pos=mass1.pos+mass1.p*dt/mass1.m | |
mass2.p=mass2.p+fnet2*dt | |
mass2.pos=mass2.pos+mass2.p*dt/mass2.m | |
# Update all the miscellaneous graphics | |
string1.axis=mass1.pos-pivot1.pos | |
string2.axis=mass2.pos-pivot2.pos | |
fgarr1.pos=mass1.pos | |
fgarr2.pos=mass2.pos | |
ftarr1.pos=mass1.pos | |
ftarr1.axis=scale*ft1 | |
ftarr2.pos=mass2.pos | |
ftarr2.axis=scale*ft2 | |
felarr1.pos=mass1.pos | |
felarr1.axis=scale*felec1 | |
felarr2.pos=mass2.pos | |
felarr2.axis=scale*felec2 | |
t=t+dt | |
# If the balls are in contact, stop the loop. | |
r=mass1.pos-mass2.pos | |
if r.mag<(mass1.radius+mass2.radius): | |
print(t) | |
t=tmax+dt | |
# Print out some parameters for data recording | |
angle=abs(asin((mass2.pos.x-pivot2.pos.x)/string2.axis.mag)) | |
print(t/tmax,r.mag/(mass1.radius+mass1.radius),(r.x)/(sep),angle) |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment