Verlet rope/snake sim

Here's an early version of my physics engine. Pre collisions-friction. The rope looks like a cylinder, but it is really just alot of sphere particles. Use F1 to pause and start, arrow keys as well as PgUp and PgDown for movement.

here it is:

[CODE]
Graphics3D 640,480,0,2

Global T.Vector=Vector()
Global T1.Vector = Vector()
Global T2.Vector = Vector()
Global dest_xang#,dest_yang#

Global Ntime,Otime,Ntsq,Lmilli,NdivO#

Const reps=4

Type structure
Field name$
Field Frst.entity,Lst.entity
Field FirstC.constraint,LastC.constraint
Field Nxt.structure
Field fric#
End Type

Global LastS.structure=Null

Type entity
Field name$,mesh
Field Nxt.entity,FirstP.verlet,LastP.verlet,FirstC.constraint,LastC.constraint,s.structure
Field center.verlet,x.verlet,y.verlet
End Type

Type verlet
Field Old.Vector,Current.Vector,Force.vector
Field ID,Nxt.verlet,e.entity
Field Anchor,Mass#,collidable
Field rep ;Del
End Type

Type constraint
Field ID,s.structure,e.entity,Nxt.constraint
Field p1.verlet,p2.verlet,Desired_Dist#
Field springiness#,collidable
Field rep
End Type

Type Vector
Field x#
Field y#
Field z#
End Type

ti=MilliSecs()
;ti=188221425
SeedRnd ti
DebugLog ti

Dim RetV.verlet(1000)

st.structure=Object.structure(VString(100,1))

l1=CreateLight(1)
RotateEntity l1,20,45,0
LightColor l1,255,255,255

Global cam=CreateCamera()
CameraRange cam,10,20000

PositionEntity cam,0,-50,-200

cube=CreateCube()
FlipMesh cube
ScaleEntity cube,200,200,200
rn=1
Global in =CreateBank(60)
While Not KeyHit(1)
l=MilliSecs()
If rn Then
update
dostructure st
positionreps
EndIf

If KeyDown(200) Then RetV(0)\Old\z#=RetV(0)\Old\z#+1 ElseIf KeyDown(208) Then RetV(0)\Old\z#=RetV(0)\Old\z#-1
If KeyDown(203) Then RetV(0)\Old\x#=RetV(0)\Old\x#-1 ElseIf KeyDown(205) Then RetV(0)\Old\x#=RetV(0)\Old\x#+1
If KeyDown(201) Then RetV(0)\Old\y#=RetV(0)\Old\y#+1 ElseIf KeyDown(209) Then RetV(0)\Old\y#=RetV(0)\Old\y#-1

If KeyHit(59) Then rn=Not rn
If Not MouseDown(3) Then
movecam
MoveEntity cam,0,0,MouseZSpeed()*2+MouseDown(1)*5+MouseDown(2)*-5
Else
MoveEntity cam,0,0,-MouseYSpeed()*2
MoveMouse GraphicsWidth()/2,GraphicsHeight()/2
EndIf

RenderWorld
Flip
ttl#=ttl+(MilliSecs()-l)
c#=c#+1
DebugLog 1000/(ttl#/c#);1000/(MilliSecs()-l)
Wend

Function VString(num,Stp#)
st.structure=Object.structure(cst("String"))
e1.entity=Object.entity(ce("STE",st))
For a=0 To num
If a=0 Then
RetV.verlet (0)=Object.verlet(crv(e1,0,0,0,0,0,0,1))
Else
RetV.verlet (a)=Object.verlet(crv(e1,0,a*-Stp#,0,Rnd(-1,1),a*-stp#+Rnd(-1,1),Rnd(-1,1),0))
c1.constraint=Object.constraint(cce(e1,RetV(a),RetV(a-1),1))
If a>1 Then c2.constraint=Object.constraint(cce(e1,RetV(a),RetV(a-2),1))
EndIf
Next
Return Handle(st.structure)
End Function

Function cube()
st.structure=Object.structure(cst("Box"))

e1.entity=Object.entity(ce("Box entity",st))
v1.verlet=Object.verlet(crv(e1,-10,-10,-10,-10+Rnd(-1,1),-10+Rnd(-1,1),-10+Rnd(-1,1)))
v2.verlet=Object.verlet(crv(e1,10,-10,-10,10+Rnd(-1,1),-10+Rnd(-1,1),-10+Rnd(-1,1)))
v3.verlet=Object.verlet(crv(e1,-10,10,-10,-10+Rnd(-1,1),10+Rnd(-1,1),-10+Rnd(-1,1)))
v4.verlet=Object.verlet(crv(e1,10,10,-10,10+Rnd(-1,1),10+Rnd(-1,1),-10+Rnd(-1,1)))
v5.verlet=Object.verlet(crv(e1,-10,-10,10,-10+Rnd(-1,1),-10+Rnd(-1,1),10+Rnd(-1,1)))
v6.verlet=Object.verlet(crv(e1,10,-10,10,10+Rnd(-1,1),-10+Rnd(-1,1),10+Rnd(-1,1)))
v7.verlet=Object.verlet(crv(e1,-10,10,10,-10+Rnd(-1,1),10+Rnd(-1,1),10+Rnd(-1,1)))
v8.verlet=Object.verlet(crv(e1,10,10,10,10+Rnd(-1,1),10+Rnd(-1,1),10+Rnd(-1,1)))
cv.verlet=Object.verlet(csrv(e1,0,0,0,Rnd(-1,1),Rnd(-1,1),Rnd(-1,1),1))
xv.verlet=Object.verlet(csrv(e1,-10,0,0,-10+Rnd(-1,1),Rnd(-1,1),Rnd(-1,1),2))
yv.verlet=Object.verlet(csrv(e1,0,10,0,Rnd(-1,1),10+Rnd(-1,1),Rnd(-1,1),3,1))

c1a.constraint=Object.constraint(cce(e1,v1,v2,1))
c2a.constraint=Object.constraint(cce(e1,v2,v4,1))
c3a.constraint=Object.constraint(cce(e1,v3,v4,1))
c4a.constraint=Object.constraint(cce(e1,v1,v3,1))
c5a.constraint=Object.constraint(cce(e1,v2,v3,1))

c1b.constraint=Object.constraint(cce(e1,v5,v6,1))
c2b.constraint=Object.constraint(cce(e1,v6,v8,1))
c3b.constraint=Object.constraint(cce(e1,v7,v8,1))
c4b.constraint=Object.constraint(cce(e1,v5,v7,1))
c5b.constraint=Object.constraint(cce(e1,v6,v7,1))

c1c.constraint=Object.constraint(cce(e1,v1,v5,1))
c2c.constraint=Object.constraint(cce(e1,v2,v6,1))
c3c.constraint=Object.constraint(cce(e1,v3,v7,1))
c4c.constraint=Object.constraint(cce(e1,v4,v8,1))

c1d.constraint=Object.constraint(cce(e1,v1,v6,1))
c2d.constraint=Object.constraint(cce(e1,v2,v8,1))
c3d.constraint=Object.constraint(cce(e1,v3,v5,1))
c4d.constraint=Object.constraint(cce(e1,v4,v7,1))

c1e.constraint=Object.constraint(cce(e1,v1,cv,1))
c2e.constraint=Object.constraint(cce(e1,v5,cv,1))
c3e.constraint=Object.constraint(cce(e1,v2,cv,1))


c1f.constraint=Object.constraint(cce(e1,cv,xv,1))
c2f.constraint=Object.constraint(cce(e1,v4,xv,1))
c4f.constraint=Object.constraint(cce(e1,v7,xv,1))

c1g.constraint=Object.constraint(cce(e1,cv,yv,1))
c2g.constraint=Object.constraint(cce(e1,v3,yv,1))
c4g.constraint=Object.constraint(cce(e1,v7,yv,1))

e1\mesh=CreateCube()
ScaleEntity e1\mesh,10,10,10
EntityColor e1\mesh,255,0,0

Return Handle(st.structure)
End Function

Function update()
If Lmilli=0 Then
Ntime=10
Otime=10
NdivO=1
Else
Otime=Ntime
Ntime=MilliSecs()-LMilli
NdivO#=Float(Ntime)/Otime
EndIf
lmilli=MilliSecs()
Ntsq=Ntime*Ntime
End Function

Function positionreps()
For a.verlet=Each verlet
PositionEntity a\rep,a\Current\x#,a\Current\y#,a\Current\z#
Next
For b.entity=Each entity
If b\mesh<>0 Then
PositionEntity b\mesh,b\center\Current\x#,b\center\Current\y#,b\center\Current\z#
AlignToVector b\mesh,b\x\Current\x#-b\center\Current\x#,b\x\Current\y#-b\center\Current\y#,b\x\Current\z#-b\center\Current\z#,1
AlignToVector b\mesh,b\y\Current\x#-b\center\Current\x#,b\y\Current\y#-b\center\Current\y#,b\y\Current\z#-b\center\Current\z#,2
EndIf
Next
End Function

Function cce(ent.entity,tp1.verlet,tp2.verlet,sp#)
c.constraint=New constraint
c\e.entity=ent.entity
c\p1.verlet=tp1.verlet
c\p2.verlet=tp2.verlet
c\springiness#=sp#
c\desired_dist#=Sqr((tp1\current\x#-tp2\current\x#)^2+(tp1\current\y#-tp2\current\y#)^2+(tp1\current\z#-tp2\current\z#)^2)
If ent\firstC.constraint=Null Then ent\firstC.constraint=c.constraint
If ent\LastC.constraint<>Null Then ent\LastC\Nxt.constraint=c.constraint
ent\LastC.constraint=c.constraint
Return Handle c.constraint
End Function

Function ccs(st.structure,tp1.verlet,tp2.verlet,sp#)
c.constraint=New constraint
c\s.structure=st.structure
c\p1.verlet=tp1.verlet
c\p2.verlet=tp2.verlet
c\springiness#=sp#
c\desired_dist#=Sqr((tp1\current\x#-tp2\current\x#)^2+(tp1\current\y#-tp2\current\y#)^2+(tp1\current\z#-tp2\current\z#)^2)
If st\firstC.constraint=Null Then st\firstC.constraint=c.constraint
If st\LastC.constraint<>Null Then st\LastC\Nxt.constraint=c.constraint
st\LastC.constraint=c.constraint
Return Handle c.constraint
End Function

Function cst(Name$)
news.structure=New Structure
news\name$=name$
news\fric#=.99
If LastS.structure<>Null Then LastS\Nxt.structure=news.structure
LastS.structure=news.structure
Return Handle news.structure
End Function

Function ce(Name$,struct.structure)
newe.entity=New entity
newe\name$=name$
newe\s.structure=struct.structure
If struct\Frst.entity=Null Then struct\Frst.entity=newe.entity
If struct\Lst.entity<>Null Then struct\Lst\NXt.entity=newe.entity
struct\Lst.entity=newe.entity
Return Handle newe.entity
End Function

Function crv(ent.entity,CX#,CY#,CZ#,OX#,OY#,OZ#,Anch=0)
v.verlet=New verlet
v\current.vector=New vector
v\current\x#=CX#
v\current\y#=CY#
v\current\z#=CZ#
v\Old.Vector=New Vector
v\Old\x#=OX#
v\Old\y#=OY#
v\Old\z#=OZ#
v\mass#=1
v\Anchor=Anch
v\rep=CreateSphere()
ScaleEntity v\rep,5,5,5
EntityColor v\rep,255,0,0
PositionEntity v\rep,CX#,CY#,CZ#
EntityAlpha v\rep,1
v\force.vector=New vector
v\force\y#=-.0001
v\e.entity=ent.entity
If ent\FirstP.verlet=Null Then ent\FirstP.verlet=v.verlet
If ent\LastP.verlet<>Null Then ent\LastP\NXt.verlet=v.verlet
ent\LastP.verlet=v.verlet
Return Handle(v.verlet)
End Function

Function csrv(ent.entity,CX#,CY#,CZ#,OX#,OY#,OZ#,SP,Anch=0)
v.verlet=New verlet
v\current.vector=New vector
v\current\x#=CX#
v\current\y#=CY#
v\current\z#=CZ#
v\Old.Vector=New Vector
v\Old\x#=OX#
v\Old\y#=OY#
v\Old\z#=OZ#
v\mass#=1
v\Anchor=Anch
v\rep=CreateSphere()
ScaleEntity v\rep,5,5,5
EntityAlpha v\rep,1
If SP=3 Then EntityColor v\rep,0,0,255 Else EntityColor v\rep,255,0,0
PositionEntity v\rep,CX#,CY#,CZ#
v\force.vector=New vector
v\force\y#=-.0001
v\e.entity=ent.entity
Select sp
Case 1:
ent\center.verlet=v.verlet
Case 2:
ent\x.verlet=v.verlet
Case 3:
ent\y.verlet=v.verlet
End Select
Return Handle(v.verlet)
End Function

Function dostructure(st.structure)
e.entity=st\Frst.entity
If e.entity<>Null Then
Repeat
doverlets e
If e.entity=st\lst.entity Then Exit Else e.entity=e\Nxt.entity
Forever
For r=1 To reps
c.constraint=st\firstC.constraint
e.entity=st\Frst.entity
If c.constraint<>Null Then
Repeat
Setdistance c\p1.verlet,c\p2.verlet,c\desired_dist#,c\springiness
If c.constraint=st\LastC.constraint Then Exit Else c.constraint=c\Nxt.constraint
Forever
EndIf
Repeat
doconstraints e
If e.entity=st\Lst.entity Then Exit Else e.entity=e\Nxt.entity
Forever
Next
EndIf
End Function

Function doentity(e.entity)
If e.entity<>Null Then
doverlets e
doconstraints e
EndIf
End Function

Function doverlets(e.entity)
v.verlet=e\firstP.verlet
If v.verlet<>Null Then
Repeat
indivVerlet v.verlet
If v.verlet=e\LastP.verlet Then Exit Else v.verlet=v\Nxt.verlet
Forever
EndIf
If e\center.verlet<>Null Then indivVerlet e\center.verlet
If e\x.verlet<>Null Then indivVerlet e\x.verlet
If e\y.verlet<>Null Then indivVerlet e\y.verlet
End Function

Function doconstraints(e.entity)
c.constraint=e\firstC.constraint
If c.constraint<>Null Then
Repeat
Setdistance c\p1.verlet,c\p2.verlet,c\desired_dist#,c\springiness#
If c.constraint=e\LastC.constraint Then Exit Else c.constraint=c\Nxt.constraint
Forever
EndIf
End Function

Function indivVerlet(v.verlet)
If v.verlet<>Null Then
If v\Anchor=0 Then
CloneVector(T,v\current)
MulVecScalar(T2,v\Force,Ntsq) ; a*timestep*timestep (Acceleration offset)
Addvector(T1,v\current,T2) ; N+a*timestep*timestep (Acceleration with Current Position)
Subvector(T2,v\Current,v\Old) ; N-O (Uncorrected Velocity)
MulvecScalar2(T2,NdivO#*v\e\s\fric#) ; (Newtime/Oldtime)(N-O) (Time Corrected Velocity)
Addvector(v\current,T1,T2) ; N+Newtime/Oldtime(N-O)+a*Newtime*Oldtime (Complete, time corrected physics)
CloneVector(v\old,T)
Else
clonevector(v\current,v\old)
EndIf
EndIf
End Function

Function SetDistance(p1.verlet,p2.verlet,dist#,sp#)
SubVector(T1,p1\current,p2\current)
deltalength# = Sqr((T1\x * T1\x) + (T1\y * T1\y) + (T1\z * T1\z))
If deltalength# <= 0.0 deltalength# = 0.0001
diff# = (deltalength# - dist#) / deltalength#
tmass# = p1\mass + p2\mass
MulVecScalar(T2,T1,diff# * (p2\mass / tmass)*sp#)
SubVector2(p1\current,T2)
MulVecScalar(T2,T1,diff# * (p1\mass / tmass)*sp#)
AddVector2(p2\current,T2)
End Function

Function Vector.Vector(x#=0.0,y#=0.0,z#=0.0)
v.Vector = New Vector
v\x=x
v\y=y
v\z=z
Return v
End Function

;// Vector 1 is set to Vector 2
Function CloneVector(v1.Vector,v2.Vector)
v1\x = v2\x
v1\y = v2\y
v1\z = v2\z
End Function

;// Vector Scalar Multiplication
;// Form of Vector1 = Vector2 * Scalar
Function MulVecScalar(v1.Vector,v2.Vector,s#)
v1\x = v2\x * s
v1\y = v2\y * s
v1\z = v2\z * s
End Function

;// Form of Vector1 = Vector1 + Vector2
Function AddVector2(v1.Vector,v2.Vector)
v1\x = v1\x + v2\x
v1\y = v1\y + v2\y
v1\z = v1\z + v2\z
End Function

;// Form of Vector1 = Vector1 - Vector2
Function SubVector2(v1.Vector,v2.Vector)
v1\x = v1\x - v2\x
v1\y = v1\y - v2\y
v1\z = v1\z - v2\z
End Function

;// Form of Vector1 = Vector1 + Vector2
Function AddVector(v1.Vector,v2.Vector,v3.vector)
v1\x = v2\x + v3\x
v1\y = v2\y + v3\y
v1\z = v2\z + v3\z
End Function

;// Vector Subtraction
;// Form of Vector1 = Vector2 - Vector3
Function SubVector(v1.Vector,v2.Vector,v3.Vector)
v1\x = v2\x - v3\x
v1\y = v2\y - v3\y
v1\z = v2\z - v3\z
End Function

;// Form of Vector1 = Vector1 * Scalar
Function MulVecScalar2(v1.Vector,s#)
v1\x = v1\x * s
v1\y = v1\y * s
v1\z = v1\z * s
End Function

; the rest is from Rob hutchison
Function movecam()
mxs=MouseXSpeed()/2
mys=MouseYSpeed()/2

; Update destination camera angle x and y values
dest_xang#=dest_xang#+mys
dest_yang#=dest_yang#-mxs

; Curve camera angle values towards destination values
xang#=CurveValue#(xang#,dest_xang#,5)
yang#=CurveValue#(yang#,dest_yang#,5)

RotateEntity cam,xang#,yang#,0
MoveMouse GraphicsWidth()/2,GraphicsHeight()/2
End Function

Function CurveValue#(current#,destination#,curve)
current#=current#+((destination#-current#)/curve)
Return current#
End Function
[/CODE]

If you use it, please credit bot builder(myself).

[CODE]
;Change

st.structure=Object.structure(VString(100,1))

;to

st.structure=Object.structure(cube())

;and take out
If KeyDown(200) Then RetV(0)\Old\z#=RetV(0)\Old\z#+1 ElseIf KeyDown(208) Then RetV(0)\Old\z#=RetV(0)\Old\z#-1
If KeyDown(203) Then RetV(0)\Old\x#=RetV(0)\Old\x#-1 ElseIf KeyDown(205) Then RetV(0)\Old\x#=RetV(0)\Old\x#+1
If KeyDown(201) Then RetV(0)\Old\y#=RetV(0)\Old\y#+1 ElseIf KeyDown(209) Then RetV(0)\Old\y#=RetV(0)\Old\y#-1
[/CODE]

Also, to change accuracy, increase or decrease reps.

Fantastic stuff!!!!

This is an option to the actual ODE wrapper indeed :)

BTW, where are collisions ?

Thank you. Heh. Collisions aren't here yet. that's why I said pre collisions-friction. I could use blitze's built in collision detection, but they don't have line to poly collisions :( . I guess I could do with a long cylinder but blitz doesn't have that either. linepicks to slow. I need somthing that runs in less that 1 tenth of a millisecond, as I expect to run it hundreds of times per-loop, yet the physics should run at 60-100 fps, allowing the game some proccessing time to. I also plan to run diffferent objects at different rates. So far away objects have less physics detail. I also plan on having the far away objects update less accuratly. Working on the collisions now......

thats pretty cool stuff. good work. look forward to seeing it evlolve.

Heh. Thanks. Right now my main focus is on collisions, as I said. I've come up with a giant system
that should hopefully run pretty fast. It will basically be a whole custom collision system that
sadly won't be very user-friendly unless I create a custom editer. basically, it will be a large
hierchy of collision checks, all box checks, that contain quad checks, triangle checks, more
box checks, cylinder checks, and sphere checks. The hierarchy of collision checks is
similar to a portal system(I think). It checks whether or not an entity could be in an area,
and if it is then it goes on to another check that involves a smaller area, until it has such
few colision possibilities that they run quickly. I've also decided to specialize the collision
system to deal not only with triangles, but with quads, as it is fast to check if a quad is
being collided with. Especially in the case that the quad is horizontal or vertical. In fact,
I have seperate type structures to deal with quads aligned to the different axis. I haven't
done any work with the other collision checks yet, as they are quite a bit harder.

Once I get collisions and friction sorted out, evolving constraints are in, then
possibly evolving Verlets, then several editers to edit the collision structures,
verlet positions, etc. So basically, I have my work cut out :) .

my god man... you do have your work cut out for you. im especially interested in what you come up with for collision detection AND response. i havent had much luck trying to make my own.