Bump! I wrote a super fancy renderer for a youtube video of this, so I thought I'd share it.
You need advanced 2D and... I think the matrix1utils to run it.
image of the program
type phyobj
x as float
y as float
vx as float
vy as float
fx as float
fy as float
mass as float
radius as float
obj
endtype
dim plot1height(-1) as float
dim plot1time(-1) as float
global numplot1 as integer
global numplot2 as integer
numplot1=-1
numplot2=-1
dim phyarr(-1) as phyobj
global phynum as integer
phynum=-1
global gaccel as float //acceleration due to gravity in meters per second. .1 is a good value, but 9.81 is realistic.
global chainnum as integer //number of links in the chain
global chainlength as float //length of the chain in meters
global chainmass as float //total mass of the chain in kilograms (mass of one link=chainmass/chainnum)
global chainK as float //the spring constant of each link (the k in F=-kx). To model a chain, this should be
//VERY high. I've had it at 5000 with a step size of .00005 and it still had significantly more give in it
//than a real chain. I think the spring constant of a steel ball is around 100,000 or so.
global chainDampening as float //dampen the velocity in the direction of chain collision. Should be very high.
//If this value were zero, the mass on the chain (bungee jumper) would return to its original position.
global objmass as float //mass of the falling objects in kilograms
global startheight as float //starting height of the two masses in meters
global separation as float //separation of the chain's first link and the falling object in meters
global floorK as float //floor spring constant
global floorDampening as float //floor dampening constant
global floorFriction as float
global elapsed as float //total time elapsed in seconds
global stepsize as float
global iterationsPerLoop as integer
global viewWidth as float //Since the motion along the vertical is greater than along the horizontal,
//the viewport has been stretched horizontally so that all the interesting parts are visible.
set display mode 1280,720,32,1,8,0
global viewTopBuffer as float
global viewBotBuffer as float //# of meters visible above and below the start point and the ground, respectively
global chainFriction as float
elapsed=0
gosub realistic_setup
drawGraph=1
createChain()
createFreefallObjs()
hold=1
sync on
sync rate 0
autocam off
position camera 0,viewwidth/2,50,-70
point camera 0,viewwidth/2,50,0
set camera view 0,0,0,screen width()/2,screen height()
set camera aspect screen height()*2.0/screen width()
tmpobj=find free object()
make object box tmpobj,viewwidth,.5,5
position object tmpobj,viewwidth/2,0,0
color object tmpobj,rgb(255,0,0)
global centerX as float`graph center, in units
global centerY as float
global centerxpix as integer
global centerypix as integer
global PPU as float `pixels per unit (so larger PPU means more zoomed in, smaller PPU means zoomed out)
global DPU as float `divisions per unit ("divisions" are the little lighter white lines dividing up a 1x1 unit)
PPU=45
DPU=4
centerX=0
centerY=0
centerYPix=Screen Height()-20
centerXPix=20
genGraphPaperImage(screen width()/2,screen height())
a2SetLineAA 1
set text size 20
do
elapsed=elapsed+stepsize*iterationsPerLoop
ink rgb(200,0,0),0
text 0,0,"Time Elapsed: "+str$(elapsed)
paste image 1,screen width()/2,0
ink rgb(0,0,0),0
text Screen Width()/2+30,0,"Acceleration Graph"
if spacekey() then hold=0
for n=0 to iterationsPerLoop
updateGravity()
calcFloorCollision()
updateChainPhy()
HoldParticle(0)
if hold
holdParticle(chainnum)
holdParticle(chainnum+1)
endif
if n=iterationsPerLoop
addPoint1(phyarr(chainnum).fy/phyarr(chainnum).mass,elapsed)
endif
integrate(stepsize)
next
if shiftkey()
drawgraph=1
else
drawgraph=0
endif
drawGraph1()
if drawgraph
drawTickmarkScale()
endif
drawPhyobjs()
sync
loop
exaggerated_setup:
//high number of links, high spring constants, low initial separation, chain length
//shorter than height, 100kg chain, 70kg mass, very small step size. View width smaller
//than view height so the horizontal movement is greatly exaggerated.
floorFriction=.01
gaccel=.2
chainnum=100
chainlength=100
chainmass=400
chainK=200
chainDampening=100
chainFriction=.01
objmass=100
startheight=90
separation=4
floorK=chainK
floorDampening=10
viewTopBuffer=10
viewBotBuffer=10
viewWidth=20
stepsize=.0001
iterationsPerLoop=1000
return
realistic_setup:
//high number of links, high spring constants, low initial separation, chain length
//shorter than height, 100kg chain, 70kg mass, very small step size. View width equal to
//view height, so the view is realistic
floorFriction=.01
gaccel=9.81
chainnum=200
chainlength=70
chainmass=100
chainK=1000000 //spring coefficient of a million!
chainDampening=5000
chainFriction=.01
objmass=100
startheight=90
separation=1
floorK=chainK
floorDampening=1000
viewTopBuffer=10
viewBotBuffer=10
viewWidth=120
stepsize=.00005
iterationsPerLoop=1000
return
accurate_setup:
//high number of links, high spring constants, low initial separation, chain length
//shorter than height, 100kg chain, 70kg mass, very small step size.View width equal to
//view height, so the view is realistic
floorFriction=.01
gaccel=.1
chainnum=200
chainlength=70
chainmass=100
chainK=5000
chainDampening=5000
chainFriction=.01
objmass=70
startheight=90
separation=1
floorK=chainK
floorDampening=1000
viewTopBuffer=10
viewBotBuffer=10
viewWidth=120
stepsize=.00005
iterationsPerLoop=1000
return
function drawTickmarkScale()
ink rgb(255,255,255),0
for n=-5 to 15
h=100*(1+n/gaccel)
line 0,h,10,h
next n
endfunction
function drawGraph1()
ink rgb(255,0,0)
if numplot1<1 then exitfunction
ylast as integer
xlast as integer
ylast=yUnitToPixel(plot1height(0))
xlast=xUnitToPixel(0)+screen width()/2
for x=1 to numplot1
ynew=yUnitToPixel(-plot1height(x))
xnew=xUnitToPixel(plot1time(x))+screen width()/2
for n=-1 to 1
for m=-1 to 1
a2line xlast+n,ylast+m,xnew+n,ynew+m,rgb(0,0,0)
next m
next n
xlast=xnew
ylast=ynew
next x
endfunction
function addPoint1(v as float, t as float)
if t<14
array insert at bottom plot1height()
array insert at bottom plot1time()
numplot1=numplot1+1
plot1height(numplot1)=v
plot1time(numplot1)=t
endif
endfunction
function HoldTopParticles()
HoldParticle(0)
endfunction
function HoldParticle(n)
increaseForce(n,-phyarr(n).fx,-phyarr(n).fy)
endfunction
function updateChainPhy()
for n=0 to chainnum-1
m=n+1
diffx#=phyarr(m).x-phyarr(n).x
diffy#=phyarr(m).y-phyarr(n).y
d#=sqrt(diffx#*diffx#+diffy#*diffy#) //distance from n to m
sumrad#=2*phyarr(n).radius
if d#>sumrad#
diffx#=diffx#/d# //now diff is the direction from n to m
diffy#=diffy#/d#
vrelx#=phyarr(m).vx-phyarr(n).vx //velocity of m relative to n
vrely#=phyarr(m).vy-phyarr(n).vy
vrel#=diffx#*vrelx#+diffy#*vrely# //magnitude of relative velocity in the direction of collision
//magnitude of relative velocity in the direction 90 degrees CCW from the direction of collision
vrel2#=-diffy#*vrelx#+diffx#*vrely#
sep#=d#-sumrad#
kforce#=-sep#*chainK
dforce#=-vrel#*chainDampening
frictionforce#=-sgn(vrel2#)*chainFriction*(kforce#+dforce#)
increaseForce(n,-diffx#*(kforce#+dforce#)-frictionforce#*diffy#,-diffy#*(kforce#+dforce#)+frictionforce#*diffx#)
increaseForce(m, diffx#*(kforce#+dforce#)+frictionforce#*diffy#, diffy#*(kforce#+dforce#)-frictionforce#*diffx#)
remstart
dforcex#=-vrelx#*chainDampening
dforcey#=-vrely#*chainDampening
increaseForce(n,-diffx#*kforce#-dforcex#,-diffy#*kforce#-dforcey#)
increaseForce(m,diffx#*kforce#+dforcex#,diffy#*kforce#+dforcey#)
remend
endif
next n
endfunction
function createFreefallObjs()
addPhyObj(objmass,viewWidth/2+separation/2.0,startheight,20)
addPhyObj(objmass,viewWidth/2+separation/2.0,startheight,20)
endfunction
function createChain()
startx#=viewWidth/2-separation/2.0
starty#=startheight
midx#=startx#+separation/2.0
midy#=startheight-sqrt(chainlength*chainlength-separation*separation)/2.0
endx#=startx#+separation
endy#=starty#
massper#=chainmass/chainnum
radius#=.5*chainlength/chainnum
for n=0 to chainnum-1
t#=n*1.0/chainnum
if t#<=.5
addPhyObj(massper#,(midx#-startx#)*2*t#+startx#,(midy#-starty#)*2*t#+starty#,radius#)
else
addPhyObj(massper#,(endx#-midx#)*(t#-.5)*2+midx#,(endy#-midy#)*(t#-.5)*2+midy#,radius#)
endif
next n
endfunction
function drawPhyobjs()
remstart
ink rgb(0,255,0),0
for n=0 to phynum
x=xWorldToScreen(phyarr(n).x)
y=yWorldToScreen(phyarr(n).y)
if n=chainnum then ink rgb(255,0,0),0
if n=chainnum+1 then ink rgb(0,0,255),0
circle x,y,phyarr(n).radius
next n
remend
endfunction
function addPhyobj(mass as float, x as float, y as float,radius as float)
array insert at bottom phyarr()
phynum=phynum+1
phyarr(phynum).mass=mass
phyarr(phynum).x=x
phyarr(phynum).y=y
phyarr(phynum).vx=0
phyarr(phynum).vy=0
phyarr(phynum).fx=0
phyarr(phynum).fy=0
phyarr(phynum).radius=radius
phyarr(phynum).obj=find free object()
make object sphere phyarr(phynum).obj,radius^.1
position object phyarr(phynum).obj,(x-viewwidth/2)*8+viewwidth/2,y,0
color object phyarr(phynum).obj,rgb(0,255,0)
endfunction
function xWorldToScreen(x as float)
x=x*screen width()/viewWidth
endfunction x
function yWorldToScreen(y as float)
y=screen height()-(y+viewTopBuffer)/(startheight+viewTopBuffer+viewBotBuffer)*screen height()
endfunction y
function integrate(dt as float)
for n=0 to phynum
phyarr(n).x=phyarr(n).x+phyarr(n).vx*dt+phyarr(n).fx/phyarr(n).mass*.5*dt*dt
phyarr(n).y=phyarr(n).y+phyarr(n).vy*dt+phyarr(n).fy/phyarr(n).mass*.5*dt*dt
phyarr(n).vx=phyarr(n).vx+phyarr(n).fx/phyarr(n).mass*dt
phyarr(n).vy=phyarr(n).vy+phyarr(n).fy/phyarr(n).mass*dt
zeroForce(n)
position object phyarr(n).obj,(phyarr(n).x-viewwidth/2)*7+viewwidth/2,phyarr(n).y,0
next n
endfunction
function calcFloorCollision()
for n=0 to phynum
if phyarr(n).y<0
normf#=phyarr(n).y*floorK+phyarr(n).vy*floorDampening
fricf#=sgn(phyarr(n).vx)*floorFriction*normf#
increaseForce(n,fricf#,-normf#)
endif
next n
endfunction
function updateGravity()
for n=0 to phynum
increaseForce(n,0,-gaccel*phyarr(n).mass)
next n
endfunction
remstart
a special function to do the equivalent of:
phyarr(n).fx=phyarr(n).fx+xamt
phyarr(n).fy=phyarr(n).fy+yamt
The reason I'm using a special function for this is so that it will be easier to calculate
the stress on an individual particle if you're not modifying forces directly.
(a particle with 0 accel and therefore 0 net force can still have huge stresses on it, so you
have to look at all the components of the force, not just the net force of .fx and .fy)
remend
function increaseForce(n as integer, xamt as float, yamt as float)
phyarr(n).fx=phyarr(n).fx+xamt
phyarr(n).fy=phyarr(n).fy+yamt
endfunction
function zeroForce(n)
phyarr(n).fx=0
phyarr(n).fy=0
endfunction
function ready_GraphUtils(graphimg)
centerX=0
centerY=0
PPU=100
DPU=4
graphimage=graphimg
endfunction
function xUnitToPixel(x as float)
x=(x-centerX)*PPU+centerxpix
endfunction x
function xPixelToUnit(x as float)
x=(x-centerxpix)/PPU+centerX
endfunction x
function yUnitToPixel(y as float)
y=-(y-centerY)*PPU+centerypix
endfunction y
function yPixelToUnit(y as float)
y=-(y-centerypix)/PPU+centerY
endfunction y
`generate the graph paper image
function genGraphPaperImage(imgw as integer, imgh as integer)
make image 1,imgw,imgh
set draw target 2,1
`advanced 2D batch drawing commands, because the built-in command is slow.
a2StartDotBatch imgw*imgh
for x=0 to imgw-1
for y=0 to imgh-1
a2dot x,y,rgb(225+rnd(10),225+rnd(10),230+rnd(10))
next y
next x
a2endbatch
`get the leftmost and rightmost pixels on the screen
leftUnit as float
rightUnit as float
upUnit as float
downUnit as float
leftUnit=floor(xPixelToUnit(0))
rightUnit=ceil(xPixelToUnit(imgw))
upUnit=ceil(yPixelToUnit(0))
downUnit=floor(yPixelToUnit(imgh))
`handle vertical bars:
pos#=leftUnit
num=0
while pos#<=rightUnit `for every vertical bar...
xpix=xUnitToPixel(pos#) `x pixel position
if (num mod DPU)=0 `regular unit divisions are one color
ink rgb(180,180,190),0
line xpix,0,xpix,imgh
else `subunits are lighter
ink rgb(210,210,230),0
line xpix,0,xpix,imgh
endif
pos#=pos#+1.0/DPU `increase the counters by a constant.
num=num+1
endwhile
`handle horizontal bars:
pos#=downUnit
num=0
`for every horizontal vertical bar
while pos#<=UpUnit
ypix=yUnitToPixel(pos#)
if (num mod DPU)=0
ink rgb(180,180,190),0
line 0,ypix,imgw,ypix
else
ink rgb(210,210,230),0
line 0,ypix,imgw,ypix
endif
pos#=pos#+1.0/DPU
num=num+1
endwhile
cx#=xUnitToPixel(0)
cy#=yUnitToPixel(0)
if 0<cx# and cx#<imgw
ink rgb(0,0,0),0
line cx#,0,cx#,imgh
endif
if 0<cy# and cy#<imgh
ink rgb(0,0,0),0
line 0,cy#,imgw,cy#
endif
set draw target 1,0
endfunction
