There's probably a better way to display the code but I'll just dump it here:
from __future__ import division
from visual import *
from visual.graph import *
## INITIAL CONDITIONS
r0 = 3 ## meters, relaxed length of spring
theta0 = math.radians(106.35) ## free angle
h = 5 ## meters, initial height of ball
vx0 = 6 ## initial horizontal velocity
## DEFINE SPRING
initialaxis = vector(r0*math.cos(theta0),r0*math.sin(theta0),0)
spring = helix(pos=(initialaxis.x,h - initialaxis.y,0), axis=initialaxis, radius=0.6,color=color.yellow)
spring.constant = 90 ## N/m
## DEFINE BALL
ball = sphere(pos=(0,h,0), radius=1,make_trail = True)
ball.mass = 1 # kg
ball.velocity = vector(vx0,0,0)
ball.acceleration = vector(0,0,0)
ball.force = vector(0,0,0)
floor = box(size=(50,.1,2),pos=(0,0,0))
## OTHER DEFINITIONS
accelOfGrav = vector(0,-9.8,0) ## acceleration of gravity
K0 = 0.5 * ball.mass * vx0**2
Ugrav0 = ball.mass * 9.8 * h
energy0 = K0 + Ugrav0 ## joules
dt = 0.005
t = 0
## GRAPHING
graph1 = gdisplay()
energygraph = gcurve(color=color.red)
Kgraph = gcurve(color=color.green)
Ugravgraph = gcurve(color=color.white)
Uspringgraph = gcurve(color = color.yellow)
positiongraph = gcurve(color=color.cyan)
while True:
rate(100)
t += dt
ball.force = ball.mass * accelOfGrav
## STANCE PHASE
if ball.pos.y < initialaxis.y:## when spring touches floor
spring.pos.y = 0.1
spring.axis = ball.pos - spring.pos##update spring length
spring.displacement = initialaxis.mag*norm(spring.axis) - spring.axis
## ENERGY ADDITIONS/SUBTRACTIONS
dEnergy = energy0 - energy
springForce = spring.constant * spring.displacement
ball.force += springForce ## update force on ball
## FLIGHT PHASE
else:
spring.pos = ball.pos - initialaxis ##update spring position
spring.axis = initialaxis
spring.displacement = vector(0,0,0)
ball.acceleration = ball.force/ball.mass
ball.velocity += ball.acceleration * dt ##update ball velocity
ball.pos += ball.velocity * dt ##update ball position
## ENERGY ACCOUNTING
K = 0.5 * ball.mass * ball.velocity.mag**2
Ugrav = ball.mass * 9.8 * ball.pos.y
Uspring = 0.5 * spring.constant * spring.displacement.mag**2
energy = K + Ugrav + Uspring #calculate current energy
## UPDATE GRAPHS
energygraph.plot(pos=(ball.pos.x,energy))
positiongraph.plot(pos=(ball.pos.x, ball.pos.y))
Kgraph.plot(pos=(ball.pos.x,K))
Ugravgraph.plot(pos=(ball.pos.x,Ugrav))
Uspringgraph.plot(pos=(ball.pos.x,Ugrav))
I created graphs, did the energy accounting section, and made the graphs update. I also went and commented everything when I got slightly frustrated.
I had problems with spring potential energy. Like regarding placement of the code or whatnot. I suppose if I actually study algorithms or programming or something these kinds of problems won't appear.
Red is total energy (K+U). Yellow is spring potential energy. Green is kinetic energy. Gravitational potential energy is supposed to be plotted but it's not there. Blue is just position of the ball (x, y). Everything is plotted to the x-position of the ball, not time.
No comments:
Post a Comment