More verlet stuff

As I mentioned in my previous topic, I was going to mess around with learning verlets. I think I've got it, take a look at this little program I've just written:

Graphics 1280, 768, 32, 2
SetBuffer BackBuffer()

SeedRnd MilliSecs()

Const TIMESTEP# = .995

Type dat
	Field x#, y#
	Field h#, v#
	Field oldx#, oldy#
	Field fric#
	Field mass#
End Type

Global p.dat
Global temp.dat


Type node_dat
	Field x1#, y1#, x2#, y2#
End Type

Global node.node_dat



For n=0 To 19
	p.dat = New dat
	
	p\x = Rand(GraphicsWidth())
	p\y = Rand(GraphicsHeight())
	p\oldx = p\x
	p\oldy = p\y
	p\h = Rnd(-1, 1)
	p\v = Rnd(-1, 1)
	p\fric = Rnd(.9, 1)
	p\mass = Rnd(1, 30)
Next

; An extra particle for creating a central 'sun', etc.
p.dat = New dat
	p\x = GraphicsWidth() / 2
	p\y = GraphicsHeight() / 2
	p\oldx = p\x
	p\oldy = p\y
	p\h = 0
	p\v = 0
	p\fric = 0
	p\mass = 100


Repeat
	Cls
	
	For p.dat = Each dat
		p\v = 0
		p\h = 0
		
		; Calculate pull for each particle interaction
		For temp.dat = Each dat
			xdiff# = temp\x - p\x
			ydiff# = temp\y - p\y
			
			dist# = Sqr(xdiff * xdiff + ydiff * ydiff)

			If Abs(dist#) < 100 Then
				node.node_dat = New node_dat
				node\x1 = p\x: node\y1 = p\y
				node\x2 = temp\x: node\y2 = temp\y
			End If

			
			If Abs(dist#) > 0 Then
				p\h = p\h + (xdiff / dist) / (dist * p\mass / 10)
				p\v = p\v + (ydiff / dist) / (dist * p\mass / 10)
				;p\h = p\h + (xdiff / dist) / p\mass
				;p\v = p\v + (ydiff / dist) / p\mass
			End If
		Next


		; Verlet
		temp_x# = p\x
		temp_y# = p\y
		
		p\x = p\x + p\x * p\fric - p\oldx * p\fric + p\h * TIMESTEP
		p\y = p\y + p\y * p\fric - p\oldy * p\fric + p\v * TIMESTEP
		
		If p\x < 0 Then p\x = 0
		If p\x > GraphicsWidth() Then p\x = GraphicsWidth()
		If p\y < 0 Then p\y = 0
		If p\y > GraphicsHeight() Then p\y = GraphicsHeight()
		
		p\oldx = temp_x#
		p\oldy = temp_y#		
	Next
	
	; Draw node lines
	Color 32, 32, 64
	For node.node_dat = Each node_dat
		Line node\x1, node\y1, node\x2, node\y2
		Delete node
	Next
	
	; Draw particles
	For p.dat = Each dat
			Color 128, 128, 255
			radius# = p\mass * .5
			Oval p\x - radius, p\y - radius, p\mass, p\mass
	Next

	Flip
Until KeyHit(1)
End

It simply simulates gravity around a point and an attempt to put some node stuff in there as well. Have a play around with a values, it's kinda fun. :)

EDIT: Note to self: Use more comment lines. ;)

Thanks. this is very interesting stuff and something I'll look into when I finish my current project.

One thing that's confusing me after a little look is adding a counter for nodes, I'm making about 4 times as many nodes as lines visible. Your not making checks between objects and self though, or duplicate A<>B and B<>A, so I guess that's where the extra lines are coming from, and this also means node resolving is happening twice per particle, right?

As an extra note, you don't need Abs(dist#) following a square root(dist#) as a square root is always positive (unless you're dealing with imaginary numbers).

You're correct. I'm not doing any checking for duplicate nodes at the moment, the nodes thing was added in afterwards as just an experiment. I'm wondering how I could do it actually. I'm thinking perhaps the node type could store the handles of both the particles, and when creating a new node it checks to see if a node connecting those two particles already exists. If it works it might be an interesting way of making AI follow paths or something similar.

The square-root thingy, that's just a bad habit. Ignore it. :D

Yup, work by the look of things - take a look:

Graphics 1280, 768, 32, 2
SetBuffer BackBuffer()

SeedRnd MilliSecs()

Const TIMESTEP# = .995

Type dat
	Field x#, y#
	Field h#, v#
	Field oldx#, oldy#
	Field fric#
	Field mass#
	Field rad#
End Type

Global p.dat
Global temp.dat


Type node_dat
	Field x1#, y1#, x2#, y2#
	Field link[2]
End Type

Global node.node_dat


For n=0 To 29
	p.dat = New dat
	
	p\x = Rand(GraphicsWidth())
	p\y = Rand(GraphicsHeight())
	p\oldx = p\x
	p\oldy = p\y
	p\h = Rnd(-1, 1)
	p\v = Rnd(-1, 1)
	p\fric = 1;Rnd(.9, 1)
	p\mass = Rnd(5, 80)
	p\rad = p\mass * .1
	If p\rad < 1 Then p\rad = 1
Next

Repeat
	Cls
	
	For p.dat = Each dat
		p\v = 0
		p\h = 0
		
		; Calculate pull for each particle interaction
		For temp.dat = Each dat
			xdiff# = temp\x - p\x
			ydiff# = temp\y - p\y
			
			dist# = Sqr(xdiff * xdiff + ydiff * ydiff)

			If dist# < 100 And p <> temp Then
				node_duplicated = False
				For node.node_dat = Each node_dat
					If node\link[0] = Handle(p.dat) And node\link[1] = Handle(temp.dat) Then node_duplicated = True
					If node\link[0] = Handle(temp.dat) And node\link[1] = Handle(p.dat) Then node_duplicated = True
				Next
				
				If node_duplicated = False Then
					node.node_dat = New node_dat
					node\x1 = p\x: node\y1 = p\y
					node\x2 = temp\x: node\y2 = temp\y
					node\link[0] = Handle(p.dat)
					node\link[1] = Handle(temp.dat)
				End If
			End If

			
			If dist# > 0 Then
				;p\h = p\h + (xdiff / dist) / (dist * p\mass / 10)
				;p\v = p\v + (ydiff / dist) / (dist * p\mass / 10)
				p\h = p\h + (xdiff / dist) / p\mass
				p\v = p\v + (ydiff / dist) / p\mass
			End If
		Next


		; Verlet
		temp_x# = p\x
		temp_y# = p\y
		
		p\x = p\x + p\x * p\fric - p\oldx * p\fric + p\h * TIMESTEP * TIMESTEP
		p\y = p\y + p\y * p\fric - p\oldy * p\fric + p\v * TIMESTEP * TIMESTEP
		
		If p\x < 0 Then p\x = 0
		If p\x > GraphicsWidth() Then p\x = GraphicsWidth()
		If p\y < 0 Then p\y = 0
		If p\y > GraphicsHeight() Then p\y = GraphicsHeight()
		
		p\oldx = temp_x#
		p\oldy = temp_y#		
	Next
	
	; Draw node lines
	Color 32, 32, 64
	nodes = 0
	For node.node_dat = Each node_dat
		nodes = nodes + 1
		Line node\x1, node\y1, node\x2, node\y2
		Delete node
	Next
	
	; Draw particles
	For p.dat = Each dat
			Color 128, 128, 255
			radius# = p\rad * .5
			Oval p\x - radius, p\y - radius, p\rad, p\rad
	Next

	Text 10, 10, "Nodes: " + nodes

	Flip
Until KeyHit(1)
End

I found that easier than what I thought it would be!

I might try and simulate accretion next. That should be pretty easy to do I think, when two particles are at a certain distance and a slow enough velocity so that they cannot escape each others pull they could be considered as 'one' particle.

Here's a challenge that I haven't solved with a look at this verlet shennanigans - have a 'rope' stretched across the screen. It should sags a bit due to gravity. And then jiggle one end so that the jiggles move through the length of the rope.

By editing out the x component to get a 2D verlet engine in only the y direction, the transmission of the impulse on one end along the length of the rope dies out quick, and I get creepage as the rope drifts, and generally it doesn't work. I also can't see where to specify a range constraint so a point can't move further than n away from it's neighbours. So all in all I'm not sure I get how this works! ;) eg. It seems increasing the value of mass decreases the 'weight' of a particle.

Well I've just solved my challenge myself, by returning to the original example posted in the rope thread. Obviously it's doing something (ConstrainMass ?) that your above example isn't, so when I changed your above code I couldn't get the same ropey behaviour. I'm not sure what you're doing above is 'true' verlet-ness because you haven't got a per-particle constraint thing. Aren't you only summing and averaging per particle?

I really ought to look harder at the papers, but it'd be nice if they gave more code examples rather than maths equations!

Hmm... I've used the formulas given in the Gamasutra article. I even did a little test in Excel to see how the 'Verlet' system I programmed above works different to the simpler Euler system I used in the previous thread.
This is how it works as far as I can gather; Using verlets the velocity variable isn't cumulative. It's just a sum of all the forces acting on the particle during that cycle. After that it resets the velocity to 0.

Say in the first cycle you set the Y velocity to 0.1; the particle will move 0.1 down the screen. In the next cycle, the velocity is set to 0, but the particle should still move 0.1 units down the screen. This is because it's using the previous coordinates to calculate the next position: the two positions should be exactly 0.1 units apart.


The rope example I posted in the other thread I don't think is 'correct' - it works to a degree, but I don't think it's the typical way of doing it. Tomorrow I'm going to have a go at doing the rope program using the 'verlet' system I used above, I have an idea of how I might do it.

I also can't see where to specify a range constraint so a point can't move further than n away from it's neighbours.


Ah - I did that using unit vectors. I don't know if you know how you use them, but I'll go over it anyway.

A unit vector is a vector with a length of exactly 1.

Take a vector v(4, 5) for example. That's 4 units up, and 5 units across. Using pythags theorem, you can work out the length to be:

vector length = sqr(4 * 4 + 5 * 5)

...Which equals approx. 6.40 units.

We want that to be a length of 1, but still point in the same direction.

To work out the X & Y coords of the unit vector you use the following:
X = X / dist
Y = Y / dist

This gives you the unit vector. To set the length of the vector, you simple multiply the coords by whatever amount you want it to be. E.g., if your vector needs to be 20 units long, then it's simply:

X = (X / dist) * 20
Y = (Y / dist) * 20

When you're calculating each location for the rope joints, you use the previous coordinates to work out where the next joint should be. You make the next joint 'n' units away from the previous point in whatever direction it's currently pointing.

I hope that makes some kind of sense! I pretty much taught myself all that over a weekend, so I'm still learning. ;)

I think you mean acceleration, rather than velocity. It's set to 0 after each cycle because acceleration generally isn't accumulative (for example, gravity, which accelerates objects downward, is constant).

Both Verlet and Euler are "integrators" (I think that's the word used...). Euler is an approximate method, while Verlet uses methods derived from the simple calculus problem of integrating the force acting on an object:

velocity = (integral of) gravity (dt)
velocity = acceleration * t + initial_velocity
position = (integral of) velocity (dt)
position = (acceleration/2) * t * t + initial_velocity * t + initial_position

(Where t = TIMESTEP)


It really starts to get fun when you connect particles together with rigid constraints, which form rigid bodys (like described in the Gamasutra article).