Search Home Members Contacts
About Us
Products
Downloads
Community
Support
Pages: [1]
  Print  
Author Topic: Basic intro to matrices  (Read 7401 times)
MPalenik
Community Member
*
Posts: 22


« on: October 29, 2004, 08:50:57 PM »

I haven't posted here for a long time - although about 3 years ago, I posted here quite a bit.  Back then I wrote some basic tutorials about physics and vectors that people seemed to like.  I know at that time there were a lot of people who considered matrices "advanced topics", although I was completely unable to get along without them, even using True Vision, especially since my game was a space based combat game, and I had to write completely new camera routines, since the camera in Truevision (at least at the time) was completely unsuited to that type of movement.

I was feeling sort of inspired today, so I thought I'd write a very light introduction on how matrices can be used for coordinate transformations.

Intro

I'll leave the actual generation of 4x4 matrices to Direct3D, since that's a slightly more complex subject.  The purpose of this tutorial is simply to make people who aren't using matrices for coordinate transformations, and simply relying on the built it properties of the engine, to become more comfortable doing so, since it allows for the implementation of many advanced effects (like mirrors, different types of billboards, and other things that you can do on your own, even when using another engine).

The purpose of this tutorial is more to give an introduction to linear algebra and its uses in describing vector spaces than to get into the specifics of how Direct3D specifically uses matrices.

Part I - Matrix operations.

If you already understand matrix multiplication and inversion, you can skip this section, if not, though, you should read through it.  Later on, I will expound upon the geometric implications of the rather mundane mathematics from this section, however, it will require an understanding of the concepts presented here, so bear with me, and hopefully the point of this section will be revealed later.

A matrix is just a collection of numbers, grouped together into a sort of rectangle, with a certain number of rows, and a certain number of columns.  The dimensions of a matrix are usually written as rowsxcolumns, so a 2x3 matrix, for example, would have 2 rows and three columns.

This means it would look something like this:

|m11 m12 m13|
|m21 m22 m23|

Notice how the rows are horizontal - from left to right, and the columns are vertical, up and down.

m11, m12, etc. are all just regular numbers that can be anything you want.

There are a couple of things we can do to this matrix - we can multiply it by a constant, add/subtract it from another matrix, multiply it by another matrix, or invert it.

Mutiplying a matrix by a constant is easy.  Every number in the matrix simply gets multiplied by a the constant.  So, for example, if we were to multiply our old 2x3 matrix (lets call it M) by a constant k, we would get

|k*m11 k*m12 k*m13|
|k*m21 k*m22 k*m23|

k, like the m's, is just a regular number that can be anything you want.

Adding one matrix to another is easy as well.  Just add the numbers from the same columb and row in the two matrices together.

in mathematical notation, where the matrices M1 and M2 are being added, to produce a new matrix, M3, this is:

M3ij = M1ij + M2ij, where i is the row, and j is the column.

To give you an example with two 2x2 matrices (this can be generalized to matrices of any size, though)

|m11 m12| + |n11 n12| = |m11+n11 m12+n12|
|m21 m22|     |n21 n22|    |m21+n21 m22+n22|

And that's all there is to it.  Note the end result is still a 2x2 matrix.

Multiplying two matrices, however, is a bit more difficult.  There are a few things you need to keep in mind.

1.  Matrix multiplication is not commutative.  So if we have 2 matrices, A and B, A*B is not equal to B*A.
2.  When multiplying two matrices the number of columns in the first matrix must match the number of rows in the second.  An easy way to do this is to write out the dimensions of your two matrices, mxn and oxp.  If the two matrices can be multiplied, n must equal o.  So for example a 3x2 and 2x4 matrix can be multiplied.  A 2x3 and 2x4 cannot, though.

It is perhaps easiest to first write the expression for matrix multiplication in Einstein notation, then explain what that means.

We have two matrices A and B, that we wish to multiply to get matric C.  In Einstein notation:

Cik = sum over i [sum over k[sum over j(Aij*Bjk)]] - note you'll see expressions like this in neural network learning methods the time.

Aij represents the i'th  row and j'th column of A.  B and C work similarly.

If you've had a lot of experience with math, you can probably decipher what that means for yourself, but if not, I'm going to explain in more detail.

We know that the number of columns in A equals the number of Rows in B.  So we can set the column we look at in A to equal the row we look at in B, and call this number n.  So, if n = 1, for example, we are looking at the first column of A an the first row of B.

Lets pick two other numbers also, o and p.  o is the row of matrix A, and p is the column of matrix B.

So, if we set n = 1, o = 2, and p = 2, we're looking at the coordinates (2,1) in matrix A, and (1,2) in matrix B.  If we set n = 1, o = 1, and p = 1, we're just looking at (1,1) in A, and (1,1) in B.

To multiply row o of A by a column p of B, pick a value of multiply Aon by Bnp, and take the sum for every value of n.  So in the most simple example, of 2x2 matrices:

if we want to multiply row 1 of A by column 1 of B, o = 1 and p=1, so

|a11 a12| * |b11 b12|
|a21 a22|    |b21 b22|

is what we're multiplying.

start with n=1, and sum between n=1 and n=2.

so we have

|a11 a12| * |b11 b12|             a11 * b11
|a21 a22|    |b21 b22|

plus (+)

|a11 a12| * |b11 b12|             a12 * b21
|a21 a22|    |b21 b22|


And that gives us a11*b11 + a12*b21

Simply multiply across matrix A and down matrix B.

This gives us element o, p of matrix C (where C = A*B)

so, this case, o, the row of the first matrix is 1, and p, the column of the second matrix is 1, so c11 = a11*b11 + a12*b21

To complete our multiplication, we need to find c12, c21, and c22.  Showing this graphically, we still have to multiply these rows and columns.

|a11 a12| * |b11 b12|
|a21 a22|    |b21 b22|

|a11 a12| * |b11 b12|
|a21 a22|    |b21 b22|

and

|a11 a12| * |b11 b12|
|a21 a22|    |b21 b22|

This can be generalized to larger matrices by applying the same rule.

Matrix inversion is slightly more complicated, so I won't cover that in very much depth right now.  I will, however, define it.

for any matrix with n rows, there is an nxn matrix I, the identity matrix, that you can multiply it by and get back the same matrix.  I is like the number 1.  If you multiply a number by one, you have the same number you started with.  If you multiply a matrix by I, you have the same matrix you started with.

So:
A*I = A, where A and I are both matrices.

The identity matrix, for a 2 row matrix is just

|1 0|
|0 1|

And for a 3 row

|1 0 0|
|0 1 0|
|0 0 1|

and so on.  Simply set the diagonal elements to 1 and all other elements to zero.

Matrix inversion can, then, be likened to the inversion of regular numbers.  Any number multiplied by its inverse will give you 1. For example, 3 * 1/3 = 1.  Similarly, in matrices, any matrix multiplied by its inverse will give you I.  So A * 1/A = I.  A inverse, however, is more typically written as A^-1.  So

A * A^-1 = I.

and also

A^-1 * A = I, despite the fact that matrix multiplication is usually noncommutative

There are several different computational methods for finding the inverse of any nonsingular matrix - and by nonsingular matrix, I mean any matrix that has an inverse.  The inverse matrix will come in very useful later on.

Part II - solving systems of linear equations with matrices.

This may seem completely unrelated to the idea of coordinate transformations, but actually, it is quite relevant.  I'm sure just about everyone at some point in middle school, or possibly earlier, was taught to use matrices to solve systems of linear equations, even if you've blocked it out of your mind by now.

A linear equation is simply an equation of the form Y = A*x1 + B*x2 + . . . + N*xn, where x1, x2, etc are all variables, and A, B, C, etc. are coefficients, or constants that those variables are multiplied by.

An example of a linear equation would be 3x1 + 4x2 + 7x3 = 10.
A simpler example would be 5x = 10
An even simpler example would be x = 2.

Notice all the variables are only to the first power, although we can have as many variables in the equation as we want, from x1 to x2, to x3, to x4, to. . . well, xinfinity.  However, if we want to find out the value of one particular variable, say for example, we want to find the value of x1 in the first equation, we have to have as many equations as there are variables.  So, in the first equation, we have x1, x2, and x3 - that's three variables.  That means we need 3 equations to actually find solutions for the variables.  Lets write them as:

Ax1 + Bx2 + Cx3 = k1
Dx1 + Ex2 + Fx3 = k2
Gx1 + Hx2 + Ix3 = k3

Taking the coefficients of our 3 equations and three variables, we can create a 3x3 matrix.

|A B C|
|D E F|
|G H I|

where A B and C are the coefficients from the first equation, and so on.

There are actually several different (although related) ways to solve this system using this matrix, but the most important invovles expilcitly finding the inverse of this matrix we have created.

What happens when we multiply that matrix we just created by this matrix?:
|x1|
|x2|
|x3|

We get a new 3x1 matrix

|Ax1 + Bx2 + Cx3|
|Dx1 + Ex2 + Fx3|
|Gx1 + Hx2 + Ix3|

This is just a matrix where each row is the left hand side of one of our equations.  And since the left hand side of an equation must equal the right hand side, it stands to reason that the matrix we have created must equal

|k1|
|k2|
|k3|

Since each row of that matrix equals the right hand side of one of our equations.

we now have:

|A B C| * |x1| = |k1|
|D E F|    |x2|     |k2|
|G H I|    |x3|     |k3|


Lets call the matrix with the A, B, C, D, E, F, G, H, I matrix M.
Now, if we take M^-1, and multiply it by both sides

             |A B C|   |x1|                 |k1|
M^-1 * |D E F| * |x2| = M^-1 * |k2|
             |G H I|    |x3|                 |k3|

so:
     |x1|                  |k1|
I * |x2| = M^-1 * |k2|
     |x3|                  |k3|

But since any matrix multiplied by I is the original matrix,

|x1|                 |k1|
|x2| = M^-1 * |k2|
|x3|                 |k3|

so after we multiply M^-1 by the k matrix, we'll get a 3x1 matrix (only one number in each row), and all we have to do is set x1 equal to the first row of this matrix, x2 to the second row, and x3 to the third row, and we'll have our solutions.

Keep this in mind.  We have one more section to go before we tie everything together.

Part III - A geometric description

Lets look at matrices and matrix multiplication in a slightly different light.  In the previous section, we made use of several 3x1 matrices, with only one column, but three rows.  A matrix of this form can also be described as a column vector.

Where we might have a vector <x, y, z>, we can write that vector as the matrix
|x|
|y|       - note the similarity to the x1, x2, x3 matrix from the previous section
|z|

We could of course, also write that vector as a 1x3 matrix, which is more similar to our original vector format.  In this case we have the 1x3 matrix
|x y z|

It is interesting to note that a 1x3 matrix and a 3x1 matrix can be multiplied, either by

|x| * |x y z|
|y|
|z|

or

|x y z| * |x|
              |y|
              |z|

The first form of multiplication isn't very useful to us for our purposes.  The second, however is actually a different way of writing a very important operation from another part of mathematics.

Say we have two vectors <x1, y1, z1> and <x2, y2, z2>.  The dot product of these two vectors is x1*x2 + y1*y2 + z1*z2.  This number is equal to the length (aka the norm) of the first vector, times the length of the second vector, times the cosine of the angle between them.

If we write our first vector as a 1x3 matrix, and our second as a 3x1 matrix

|x1 y1 z1| * |x2|
                    |y2|
                    |z2|

gives us x1*x2 + y1*y2 + z1*z2.  Look familiar?  It should, its the dot product of those two vectors.  what we've done is use matrices to perform vector operations!

The dot product is especially usefull, since it allows us to find the projection of one vector onto another.  The word "projection" in this context means exactly what it sounds like it means - it's sort of like the shadow of one vector on another.  If you were to take a line perpendicular to one vector, and find the point where it intersects both that vector and a second vector, the length along the first vector to the point of intersection is the projection of the second onto the first.

since u dot v = |u|*|v|*cos(theta) (where |u| indicates norm u, and |v| indicates norm v), the projection of v onto u is (u dot v)/|u|.  If we want that distance in terms of lengths of u, simply divide by |u| again, and we're left with (u dot v)/|u|^2.

We can even do some more interesting things, like putting several column vectors together to form a 3x3 matrix.

|x1 x2 x3|
|y1 y2 y3|
|z1 z2  z3|

Assuming these vectors are all linearly independent, we've just defined a 3x3 vector space.  That means that each column vector is one of our axis in 3 dimensional space.

The first column vector, <x1, y1, z1> is our x axis.  The second is our y axis, and the third is our z axis.

In a regular, un-rotated, unscaled coordinate system, the x axis is <1,0,0>, the y axis is <0,1,0>, and the z axis is <0,0,1>.  This gives us the matrix

|1 0 0|
|0 1 0|
|0 0 1|

which is the identity matrix!  This is called the "standard basis".  I'll talk more about what this all means in the next section.

Another, more difficult to relate to way of looking at this involves calculus, but for those interested, we can define our coordinates X Y and Z by the functions

X = x1*x + x2*y + x3*z
Y = y1*x + y2*y + y3*z
Z = z1*x + z2*y + z3*z

we then create a linear map from x y z to X Y Z by the jacobian matrix (I'll use d to denote partial derivitive, rather than full derivitive, as it usually does)

|dX/dx dX/dy dX/dz|
|dY/dx dY/dy dY/dz|
|dZ/dx dZ/dy dZ/dz|

and we're back to

|x1 x2 x3|
|y1 y2 y3|
|z1 z2  z3|

Jacobian matrices are used a lot in back proppegation methods of neural net learning - although if you read any tutorials on it, they probably won't explicitly tell you that.  But learn to recognize Einstein notation and some more advanced topics in calculus, and you'll see it right away.

In DirectX, you'll typically be using 4x4 matrices, and 4 vectors, projecting the results back into 3 dimensional space.  Because of this dimensional reduction, the operations you're using aren't simple linear maps, but the idea is very similar.

Part IV - Defining vector spaces

As was pointed out in the last section, a 3x3 vector space can be described as a linear map from one coordinate system to the standard basis.

Since the three column vectors in our 3x3 matrix are our 3 x, y, and z axis, respectively, we are free to rotate or scale our x, y, and z axis in any way we want.  These axis make up our reference frame.  If we turn the axis, it's like we're turning our perspective, changing the directions that we call left and right, forward and backward, up and down.

For example, if I want to rotate my point of view 90 degrees counter clockwise, I have to reorient my coordinate system 90 degrees counter clockwise.  My new x axis, instead of being <1,0,0> is now <0,1,0>.  My new y axis, instead of being <0,1,0> is now <0,-1,0>

So whereas our original matrix defining our vector space was

|1 0 0|
|0 1 0|
|0 0 1|

our new vector space is defined by the matrix

|0 -1 0|
|1  0 0|
|0  0 1|

The z axis is unchanged.  Ths means that we have rotated around our z axis.

The general rule in 3 dimensional space is, that to rotate around any given axis, take the two vectors of the axis you're not rotating around, and the two components of those vectors that don't belong to the axis you're rotating around - so for example, if you're rotating around the z axis, you would want the x and y components of the x and y vectors, or if you were rotating around the x axis, you would want the y and z components of the y and z vectors - and change them to <cos(theta), sin(theta)>, <-sin(theta), cos(theta)>, where theta is the angle you're rotating by.

These linear maps can also be multiplied to create new vector spaces that combine the properties of the two.  If you multiply two rotation matrices, the resulting matrix will be describe a vector space where the axis have been rotated first by one matrix, then by the other.

Part V - putting it all together

So what's the point of all this stuff?  Now, finally, you're ready to find out.

A linear map (or a 4x4 translation/rotation matrix) allows us to take a vector in one coordinate frame and represent it in another.  Since our "natural" coordinate system is described by the identity matrix, for the most part, we're going to want to find how to map vectors from the identity matrix to vectors in the space of another matrix, and vice versa.

Suppose we have a vector and we want to rotate it 90 degrees around the z axis.  If we describe the new, rotated vector in a vector space that is rotated 90 degrees around the z axis, we will have our original vector back.  However, having the vector in this rotated coordinate frame doesn't do us much good.  We need a way to convert a vector from a different vector space into ours.

Suppose our vector space is represented by a matrix A, and the rotated vector space is represented by a matrix B.  What we first need to do, is find the coordinates of our vector from B's space to some natural, identity space, and from there into our space.

For the time being, lets just assume that our space is the identity space - this will actually be true most of the time, anyway.

Assume we have some vector in B's space <x,y,z>.  Since the x axis of B in identity space is Bcolumn1, the x component of our vector in identity space is simply  x*Bcolumn1.  Similarly, the y component of our vector in identity space is y*BColumn2, and the z component in identity space is z*BColumn3.

The total description of the vector in identity space, then, is the sum of the x component in identity space, the y component in identity space, and the z component in identity space.

B, as a matrix, is written as

|x1 x2 x3|
|y1 y2 y3|
|z1 z2  z3|

As a column vector, our vector in B's space is
|x|
|y|
|z|

That same vector in identity space, as a column vector, then becomes

|x1*x + x2*y + x3*z|
|y1*x + y2*y + y2*z|
|z1*x + z2*y + z2*z |

If you look closely, this matrix should look familiar.  Remember when we talked about matrix multiplication?  This column vector is simply the matrix produced by taking

|x1 x2 x3|    |x|
|y1 y2 y3| * |y|
|z1 z2  z3|    |z|

So, if R is the vector in B's space, and V is the vector in identity space, B*R = V.


Another problem that will come up quite a bit is converting from the identity space into a different vector space.

suppose we have a matrix A that describes a vector space.  The matrix itself is:

|x4 x5 x6|
|y4 y5 y6|
|z4 z5 z6 |

We know that our vector, R, in identity space is

|x|
|y|
|z|

We need to find what that vector is in A's space, though.  However, from the transformation we just derived derived, we know that, if R is the vector in A's space and A is the vector in identity space, A*R = V, so:

|x4 x5 x6|  |v1|     |x|
|y4 y5 y6|*|v2| = |y|
|z4 z5 z6 |  |v3|     |z|

V1, V2, and V3 are unknowns that we need to solve for.

This, however, should look very familiar as well.  What we really have is a system of 3 equations:

x4*v1 + x5*v2 + x6*v3 = x
y4*v1 + y5*v2 + y6*v3 = y
z4*v1 + z5*v2 + z6*v3 = z

And fortunately, the system is already in a matrix form, so we can use the equation solving techniques described in section II to find v1, v2, and v3.

Remember from II

A^-1 * A * V = A^-1*R
I * V = A^-1*R
V = A^-1*R

This means that to find R in A's space, simply multiply A^-1 by R.

A slightly more complicated (but not too much more) operation is converting from one vector space to another that's not the identity matrix.

However, we already know how to convert from one vector space to the identity space, and from the identity space to another space, so all we need to do is combine those two operations.

So to convert from A's space to B's

First:

Convert from A to I by multiplying A by R.  We now have V = A*R

Next, convert from I to B by multiplying B inverse by V.  And since V = A*R, the full expression to convert from A to B is

Va->b = B^-1*A*R, where R is the vector in B's space that you want to convert.

Although I have described everything in terms of 3x3 rotation matrices - linear maps - this all applies in the exact same manner to DirectX's 4x4 combined translation/rotation matrices.
Logged
e!chhoernchen
Community Member
*
Posts: 327


WWW
« Reply #1 on: October 30, 2004, 03:49:21 AM »

bufftata... ull get a medal for that Wink
no it is really good, but my english isnt  lol
Logged

MPalenik
Community Member
*
Posts: 22


« Reply #2 on: October 30, 2004, 05:28:01 AM »

I'm not sure what buftatta means, but thanks.  I see that the formatting has gotten screwed up a bit, especially in the places where I multiply row vectors by column vectors, as it seems leading spaces are deleted.  I'll have to see if there's something I can do about that.
Logged
newborn
Customers
Community Member
*****
Posts: 2437


WWW
« Reply #3 on: October 30, 2004, 07:35:31 AM »

Quote from: "MPalenik"
I'm not sure what buftatta means, but thanks.  I see that the formatting has gotten screwed up a bit, especially in the places where I multiply row vectors by column vectors, as it seems leading spaces are deleted.  I'll have to see if there's something I can do about that.


yes, use the BBCode [code ] and [/code ]
Logged

SylvainTV
Administrator
Community Member
*****
Posts: 4478


WWW
« Reply #4 on: October 30, 2004, 07:56:00 AM »

Thanks for this matrix explaination. I think this should go in the ressource library on the new site Smiley Very explainful
Logged

Regards

Sylvain Dupont
TrueVision3D Developer
sylvain@truevision3d.com

TV3D IRC at http://chat.truevision3d.com or on server irc.truevision3d.com #Truevision3D. Come talk with us !
Waterman
Customers
Community Member
*****
Posts: 1141


« Reply #5 on: October 30, 2004, 04:04:29 PM »

Quote from: "SylvainTV"
Thanks for this matrix explaination. I think this should go in the ressource library on the new site Smiley Very explainful

If that's still in distant future (?), perhaps stick it now so it won't be lost? There's still space for stickies in off-topic.
Logged

Things should be described as simply as possible - but not simpler [A. Einstein]
e!chhoernchen
Community Member
*
Posts: 327


WWW
« Reply #6 on: October 31, 2004, 03:45:04 AM »

hehe   'bufftata' = im impressed, overroled of n1 stuff and so on  Smiley
Logged

MPalenik
Community Member
*
Posts: 22


« Reply #7 on: October 31, 2004, 05:33:37 PM »

Well, if this thread does go sticky, and there are any more tutorials relating to math or physics (and their application to game programming), or how certain aspects of 3D engines work (I've developed my own, albiet incomplete ones, or a few occassions, so I can probably answer most questions), I'd be happy to add more tutorials to the thread.
Logged
gappy_
Community Member
*
Posts: 74


« Reply #8 on: November 28, 2004, 07:50:30 PM »

Very impressive ,it looks very useful for me,It answers my question about matrix.
Logged

Click to view my project. Smiley                                                                                        
  By Gappy
tido
Community Member
*
Posts: 384


« Reply #9 on: November 28, 2004, 11:16:50 PM »

man,
I'm too lazy to read all that.
I'm currently learning about matrices at school (course called algèbre et géomètrie vectorielle in french) so looking at this make me thing about school.lol I may come back later to read
Logged

-tido
houde
Community Member
*
Posts: 482


WWW
« Reply #10 on: November 28, 2004, 11:23:59 PM »

cours de cégep ça no?
le peu que je sais des matrices vient de mon cours de math pour ma technique "mathématiques de l'informatique" et c'était de la marde, mais toi je crois que c'est pas mal plus poussé.. tk good luck Cheesy

Zwi zwi.. I said something you don't need to know the meaning.. as always ^^
Logged

web designer/programmer
------------------------------
tido
Community Member
*
Posts: 384


« Reply #11 on: November 29, 2004, 11:01:00 AM »

oui, c'est au cégep.  Dernier cours de math de sciences de la nature.
Logged

-tido
MPalenik
Community Member
*
Posts: 22


« Reply #12 on: April 16, 2005, 07:32:10 PM »

This tutorial will show you how to implement the Lagrange formalism of mechanics into your programs, and explain what situations it might be useful in.

Intro

The first introduction that people are usually given to physics uses Newton's laws to derive equations of motion for particles.  These equations are simple to understand and use, however, for more complex problems, it can sometimes be difficult to find the equations of motion using this formalism.  Instead, we can turn to two newer formulations of physics, the Hamiltonian and Lagrangian.  These formulations are so powerful, in fact, that they can even be applied to situations where Newton's laws would typically be completely useless, like quantum mechanics and relativity (although some difficulties do arise with relativistic Hamiltonians).

Although the latter two cases are almost certain never to come up in game physics, there are many other times when the Lagrangian and Hamiltonian may come in useful.  I'm not aware of any specific cases relating to game programming, however, I did notice that the use of the Lagrangian was mentioned in a book on game physics that I briefly looked at, so whatever the applications are, I figure they must exist.  I will try to give a few examples of where it might be better to use the Lagrangian, but they may or may not be relevant to most games.

Even outside of physics, there may be some kind of application of the Lagrangian formalism to AI.  I'm not very familiar with A* algorithms, since AI isn't an area I've studied very intensively, but I do know that they revolve around minimizing some kind of path.  The Lagrangian makes use of the calculus of variations to do something similar - to find the path that minimizes T - U, so insight into this method of calculation may also lend some insight to people interested in A*.


Part I: Background Physics Information - Forces and Potentials

There's more than one way of deriving Lagrangian mechanics, and not all of them directly involve the calculus of variations to minimize the action (T-U), but unfortunately, that is the only derivation that I can remember.  The derivation that I'm going to present is the standard derivation, and I doubt that other derivations could have significantly different interpretations.  The idea of path minimization is very important in physics as well as other areas, so before I get into any math, I think it warrants further explanation as to what exactly is going on.

In physics, we have things called force fields.  They're not quite the same as the ones you see on star trek that knock you over when you run into them, like hitting a brick wall.  Instead, they're a collection of vectors in space that cause objects in them to move around.  When an object is in the force field, it will accellerate, accellerating in the direction that the little vector in space is pointing, by an amount equal to that vector's length.

Gravity, for example, is a force field.  Objects in a gravitational field accellerate toward the body giving off the gravitational force.  In the case of the earth, we can look at the space around the earth as being filled with lots of little vectors, pointing toward the surface of the earth that get smaller by a factor of 1/r^2 (where r is the distance from the center of the earth).

Of course, there are other forces as well, namely the electromagnetic force, the strong nuclear force, and the weak nuclear force.  The principle is the same for these forces, except that all objects won't accellerate equally in them.  The force on the objects gets multiplied by some value, like charge.

This is a very good description of mechanics in general, that forces act on objects making them move.  In some sense, we can say that forces are real things that actually inhabit points in space.

However, we can define something else whose reality is slightly less definite - something that we can't actually say is a real, physical thing, like force, but which is just as useful in physics.  It's called a potential.

To lift something up in a gravitational field, or any field that actually acts on the body, it requires a certain amount of work.  The amount of work required is dependant on the mass (or charge, as in the case of an electric field) of the object and the strength of the force field.  From now on, I'm only going to talk about gravitational fields, as they are the most relevant, and it will do away with other complications that may arise as in the case of the magnetic field (as well as getting rid of the need to mention charge).

The amount of work required to lift an object can simply be calculated from the method we use to calculate work:

Work = Integral(F*dx)

F is the force acting on the object, which is the gravitational field multiplied by the mass, and dx is an infintesimal component of the path.

The meaning of that equation might seem a little bit hard to decipher to some people at first.  But all it's really telling us is that the work done on an object is equal to the force on the object multiplied by the distance over which that force is applied.  Since the actual values of force and the direction/speed of our motion can change with time, we have to transform the simple quantity F*x, force times distance into an integral, F*dx, the dot product of force and dx.  If the math isn't too clear to you, don't worry.  The only thing you really need to remember is that work is equal to the force on an object times the distance over which it's applied.

Gravitational fields are nice, because they're what we call "conservative" fields.  That means that it takes work to move up and down in a gravitational field, but not to go in any other direction.  Of course, there are more rigerous mathematical definitions (e.g. curl(g) = 0), but basically, I have relayed to you the important information about conservative fields.  Since it doesn't take any work to move left or right (or to circle around the earth, if you want to be technical), we can say that an object has a specific amount of potential energy at a certain height, no matter how it got there, and no matter where on earth it is.

This means that we could define a potential energy for the earth for every point in space, which would be a function of distance.  Instead of being made up of vectors, like a force field, it would just be made up of scalers - that is to say, numbers, with no specific direction.  However, it would still contain all of the information that the force field had in it, and we can calculate the force field from it simply by taking the gradient (if you don't know what a gradient is, don't worry.  It doesn't really matter.  Just know that it's the change in the potential energy with respect to distance expressed as a vector).

When an object moves up in a gravitational field, the force and the motion are pointing in opposite directions.  However, we still want our potential energy to increase as we go up, instead of going down.  The integral(F*dx) will be negative because F and dx have opposite signs.  This means, in defining our potential "U", we need to flip the sign on our integral, so:

U = -Integral(F*dx)

The fact that all of the information in the force field is contained within the potential should tell us that there's some way of doing physics using scaler potentials, rather than vector fields, and that is the topic that I'm going to cover in the rest of this tutorial.

There's one more quantity we need to define: Kinetic energy.  Kinetic energy is basically the same as work.  Work is simply the change in kinetic energy after a force has acted on an object.  Kinetic energy, T, is calculated by:

T = Integral(F*dx)

So how, then, is this different from the potential?  Well, there are two differences.

1. The force F doesn't have to be solely from our gravitational field, like in the potential we calcluted.  F in this case could also be from somebody pushing on the object, or a rocket firing, or even some other kind of field.
2: The sign is reversed.  This means that as an object gets lower in a gravitational field, it's kinetic energy will go up (if no other forces are acting on it).  This makes sense, because kinetic energy is related to velocity, and we all know that things speed up as they fall.

Because our kinetic energy might have some source other than the gravitational field, we need a way of defining it independant of whatever forces are acting on it.  I hinted at this in #2 - the kinetic energy is related to the velocity.  It's defined as:

T = 1/2m*v^2

This is more or less all the background physics you need to know, however there are a few details that I'd like to clear up before moving on.  It's mainly semantics, so if you're not all that interested in the physics, skip to the next section.  However, to be completely rigerous, I need to point these things out.

When we defined work ealier, as the integral of F*dx, what I neglected to mention is the fact that we have to put limits on the integral.  The integral takes place between two points, a, and b.  So the work that we've defined is only the work that it takes to get from a to b.

When we define our potentials, we want to leave the limits off of our integral, to give us something more generic that can be used throughout all space.  To be mathematically correct, though, that means that we have to add a constant + C to our integral, which could, in theory, be any number.

However, to do the physics, we can leave off the + C, because the really important thing isn't the actual value of the potential at any point in space, but rather, the way the potential changes throughout space.  When looking at the change in the potential, the + C will drop out on its own, so it is never necessarry to write it in the first place.

Second of all, the quantity F is equal to some vector field multiplied by mass, in the case of gravity.  It is often usefull to write a potential V as integral(F/m*dx), because it is more generic than what I'll call U, integral(F*dx).  U is simply V*m, and thus V can be used for any object, regardless of mass.


Part II: The Lagrangian

There's a very interesting behavior that objects obey, from which we can derive their equations of motion.  Suppose you toss a ball up into the air at a certain speed, with a certain direction.  We could figure out how it's going to move and where it's going to land using Newton's laws.  But if we look it's motion, we notice there's a very special property about it's path that is different from every other path.

At any given point in time, an object on earth has a kinetic energy, which as described in the last section, depends on it's velocity, and a potential energy, which is usually just due to it's height above the ground (let's ignore any other force fields that might be present.  Generally speaking, we only have to worry about gravity).

We can sum up the kinetic energy that our particle has, over all of the points in it's path, and do the same thing with the potential energy.  Since all real paths have an infinite number of points in them, we need to make this sum into an integral.

To make things simple, let's say that we're throwing the ball straight up into the air, so we only have one coordinate to worry about - the y coordinate, height.

If we have our function y(t), which is y coordinate as a function of time, we can calcluate from it our potential and kinetic energies.

Our potential energy, is just some function of height, which we'll call U(y).
Our kinetic energy is a function of velocity.  Velocity is a function of the time derivitive of y.  The time derivitive of y(t), dy/dt, we'll call y'.  Our kinetic energy, then, is some function T(y')

When we throw our ball up into the air, we start out with a high velocity, but very low to the ground.  That means a high kinetic energy and low potential.  As the ball gets higher and higher, the potential energy goes up, and the kinetic energy goes down.  As it comes back down, the potential energy decreases again, and the kinetic energy goes up.  We want to sum up the potential and kinetic energies at all the different times along the entire path.

The sum of our kinetic energy over the entire path is:
Integral(T(y'))

The sum of our potential energy over the entire path is:
Integral(U(y))

Now, let's take the difference of those two quantities, and call it L

L = Integral(T(y')) - Integral(U(y))

Because of the properties of integrals, we can combine them into one:

L = Integral(T(y') - U(y))

It turns out, that for the path our ball took, L is smaller (closer to zero - we should say the absolute value of L is smaller) than it would have been over any other path.

Of course, the initial kinetic energy of the ball depends on how hard we throw it.  There's no way to get around that.  But to make L smaller, the ball has to go up, so that T gets smaller and U gets bigger.  However, as it goes up, it's kinetic energy will go down, until eventually U is larger than T, and |L| starts to get bigger again.  There's no way around that.  Eventually, the ball will completely stop, and U will be big and T will be small, so it will have to reverse direction, and fall back downward, so that U gets smaller and T gets larger again.

So how do we figure out what path does this?  There's an interesting kind of calculus involved, which I'll describe below.  If you're not interested in the math, skip this next section between the ---'s, and you can see the formula at the end.
----------------------

We want to find y(t), such that L(y, y') is minimized.  Assume that we have the path y(t) that minimizes L.  We know that y starts at point a and ends at point b.

If we add some other function to it, n(t), that doesn't change the start or end points, Integral(L(y+n, y'+n')dt) should be larger than Integral(L(y,y')dt).

If we multiply n(t) by alpha, we get Integral(L(y+alpha*n, y'+alpha*n')dt).  We'll call this integral J.  

J = Integral(L(y+alpha*n, y'+alpha*n')dt)

J is maximized when alpha = 0.  The problem is, we don't actually know what y+alpha*n and y'+alpha*n' are when alpha equals zero (because we don't know what y and y' are).

To find out, we take the derivitive of J with respect to alpha, and call it J':

J' = d/dalpha Integral(L(y+alpha*n, y'+alpha*n'))

Following from regular calculus, J is maximized when J' = 0.  So:

d/dalpha Integral(L(y+alpha*n, y'+alpha*n')) = 0

since d/dalpha and the integral commute, we can rewrite this as:

Integral(d/dalpha L(y+alpha*n, y'+alpha*n')) = 0

Let's define an new function Y, such that Y = y+alpha*n, so now we have:

Integral(d/dalpha L(Y,Y'))) = 0


And:

Integral(d/dalpha L(Y, Y'))) = Integral(dL/dY*dY/dalpha + dL/dY'*dY'/dalpha) = 0

since Y = y+alpha*n, dY/dalpha = n
and since Y' = y'+alpha*n', dY'/dalpha = n'

Integral(dL/dY*n + dL/dY'*n') = 0

Since the second term, dL/dY'*n' = dL/dY'*dn/dt, we can integrate this by parts.  Doing so, gives us [dL/dY'*n] - Integral(d/dt(dL/dY')*n).  So our full integral becomes

Integral(dL/dY*n - d/dt(dL/dY')*n) + [dL/dY'*n] = 0

I put dL/dY'*n in brackets because it must be evaluated from a to b.  Since we know that n(t) goes to zero at a and b, this term goes to zero, leaving us with

Integral(dL/dY*n - d/dt(dL/dY')*n) = Integral((dL/dY - d/dt(dL/dY'))*n) = 0


In order for the integral to go to zero, (dL/dY - d/dt(dL/dY'))*n must equal zero.  Since n is an arbitrary function, we know it can't go to zero.

That means that:

dL/dY - d/dt(dL/dY') = 0

Solving for Y, we get our function y(t), which minimizes L.

----------------------

So, now we have our equations of motion.  If your problem has more than one dimension, just apply these equations to all of your dimensions.  So, in the real world, where we have three dimensions, you would solve for:

dL/dX - d/dt(dL/dX') = 0
dL/dY - d/dt(dL/dY') = 0
dL/dZ - d/dt(dL/dZ') = 0


and from that, you would get, X(t), Y(t), and Z(t).  The x, y, and z coordinates of your object as a function of time.  X', Y', and Z' denote the derivitives of these functions with respect to time (the x, y, and z components of velocity).


Part III: An Example

Let's return to the problem of throwing a ball on earth.  This time, let's not assume that it's getting thrown straight up, though.  We'll do it in 3 dimensions, and solve the problem with the Legrangian.

We know that the force of gravity on earth falls off as 1/r^2.  However, the earth is really large, and we're already really far away from the center, so our change in r, if we go up a few feet isn't that big compared to our distance from the center of the earth.  This means, that to a good approximation, the force of gravity is constant with respect to height.

For most games, you'll want the force of gravity to be constant with respect to height.  There's absolutely no reason for any land based game to have gravity fall off as 1/r^2, since if you did it right, you'd never be able to tell the difference.

I know I used feet two paragraphs ago, but I'm going to switch to meters, because it's the standard.  The accelleration due to gravity at the surface of the earth is approximately constant equalling:

g = 9.8 m/s^2

That's 9.8 meters per second per second (or about 32 ft/s^2).  That means that when you're falling, every second you fall 9.8 meters per second faster than the last.  Your accelleration, a is equal to g.

Since from Newton's laws, we know that

F = m*a and we know
a = -g, then the force of gravity is
F = -m*g

there is a negative sign infront of g because g pulls us downward.

Now, here's where you need to remember what we learned in the first section about potentials.  We need to calculate the potential, U on earth, which is the integral of -F*dx.  Since F is constant, though, we can just replace this with -F*X, force times distance.

If we lift an object off the ground to a height, y, we've moved it a distance y, and our potential becomes:

U = m*g*y

As always, our kinetic energy is:

T = 1/2*m*v^2

So our Lagrangian, L is

L = 1/2*m*v^2 - m*g*y

Now, we need to find out what  v is in terms of x, y, and z.

Well, since velocity is the vector <x', y', z'> (where ' denotes the time derivitive, as before), the length of this vector, which we'll call v, is:

v = sqr(x'^2 + y'^2 + z'^2).

Fortunately, what we're really interested in is v^2, so that nasty square root will disappear, and we get

v^2 = x'^2 + y'^2 + z'^2

This is always true, for any problem, as long as you use rectangular (x,y,z) coordinates.
so, rewriting our L in terms of our formula for v^2, we get:

L = m*(1/2(x'^2 + y'^2 + z'^2) + g*y)

plugging this into the formula we derived in the last section, we get:

dL/dx = 0
dL/dy = m*g
dL/dz = 0


dL/dx' = m*x'
dL/dy' = m*y'
dL/dz' = m*z'


and
d/dt(dL/dx') = m*x''
d/dt(dL/dy') = m*y''
d/dt(dL/dz') = m*z''
where '' denotes the second derivitive with respect to time (accelleration)

So, putting it all together:

0 - mx'' = 0
m*(g - y'') = 0
0 - mz'' = 0


Which means that:

x'' = 0
y'' = g
z'' = 0


This means that our ball isn't accellerating in the x or z directions, but it is accellerating in the z direction, at a rate of g.  This makes sense, because the force of gravity is only acting in the y direction, and since there are no other forces, that should be the only direction of accelleration.

Solving the three equations we have, to find x, y, and z, we get:

x(t) = x_0 + v_x0*t
y(t) = y_0 + v_y0*t + 1/2*g*t^2
z(t) = z_0 + v_z0*t


<x_0, y_0, z_0> is the location of the person who threw the ball, which we can specify to be anything we want.  <v_x0, v_y0, v_z0> is the velocity that the ball was thrown at, which we can also specify to be anything we want.

This is the exact same answer as you would get using Newtons laws, but deriving it a different way.  Admittedly, in this situation, it would have been just as easy, or easier to use Newtons laws, but there are many situations where it is not.  I just used this is a simple example to show you how to go about solving problems this way.


Part IV: Applying it to a game

The Lagrangian formalism seems all well and good, but the problem is, it involves derivitives and integrals, which we can't do with our computers, so it would seem that it's not very useful for computer programming.  In reality, Newton's laws rely on derivitives and integrals as well, it's just that computer programmers use approximations in games to get around actually taking the derivitives and integrals.

What we need to remember is that we can find the path that minimizes the overall L by finding the direction to head at every point in time that minimizes L at that particular time.  So we just need to get our object to follow the lagrange equations (approximately) at every cycle of our game loop.

The first thing we need to do is write some code that will hold our masless potential, V, as a function of x, y, and z.  For a gravitational field, this is easy.  It's just

Code:
Function VofX(x as single, y as single, z as single) as single

'declare g as a constant somewhere in the program
VofX = g*y

End Function


Our potential, U,  that we use in the Lagrangian, is just

U = m*VofX

where m is the mass of our particle (or vehicle - whatever object it is that you want to apply the physics to), defined somewhere else in the program.

For other fields, our potential wouldn't be so trivial.  For example, if you wanted to write a function to calculate the potential of a planet in space, for some kind of space combat game, you would do something like

Code:
Function VofX(x as single, y as single, z as single) as single

Dim r as single
r = sqr(x*x + y*y + z*z)

'G_const is not the same as g.  It's a small constant that determines how much force
'each unit of mass will put out.  Also, M_planet should be defined somewhere.
'It is the mass of the planet.

VofX = -M_planet*G_const/r

End Function


But regardless of our potential, the process is the same.  First, calculate a vector that is the change in the potential in the x, y, and z directions.

You can get this by doing:

Code:
'Vel is a vector that holds the velocity of our object

DeltaU.x = (VofX(x + Vel.x, y, z) - VofX(x - Vel.x, y, z))/(2*Vel.x)
DeltaU.y = (VofX(x, y + Vel.y, z) - VofX(x, y - Vel.y, z))/(2*Vel.y)
DeltaU.z = (VofX(x, y + .1, z + Vel.z) - VofX(x, y, z - Vel.z))/(2*Vel.z)


In cartesian coordinates, since our potential, U, isn't dependant on velocity (at least for a gravitational field.  Things get more complicated with magnetic fields), d/dt(dL/dx'), d/dt(dL/dy'), and d/dt(dL/dz') are always the same.  They are:

m*x'', m*y'', and , m*z'', respectively.

We can use this fact to calculate our accellerations, which we can then use to change our velocities.

Pluggin in our calculations into the Lagrange equations, we get

Code:

Function Update Position()

'Accel is a vector holding our accelleration.  m is our mass, as usual.

Accel.x = 1/m * DeltaU.x
Accel.y = 1/m * DeltaU.y
Accel.z = 1/m * DeltaU.z

'Now, add the accelleration to the velocity.  We could have just added 1/m * DeltaU,
'but I chose to store our accelleration first, to illustrate the point about how it is
'calculated.

Vel.x = Vel.x + Accel.x
Vel.y = Vel.y + Accel.y
Vel.z = Vel.z + Accel.z

'Now, add our velocity to our position, to get our new position

Pos.x = Pos.x + Vel.x
Pos.y = Pos.y + Vel.y
Pos.z = Pos.z + Vel.z

End function


And there you have it - an implementation of the Lagrangian formalism into your program.  All you need to do is set the function VofX to whatever particular potential you're using in your game.

So far, however, we haven't seen how implementing this formalism would be more beneficial than using Newton's.  In fact, since all we're doing in the above example is recalculating our forces from the potential, and getting accellerations from them, if we just continue using the Lagrangian like this, it will be exactly the same as an implementation of Newton's laws, only slower.

In the next section, I will show you where the true power of the Lagrangian lies, and how it can provide much easier answers to certain problems than Newton's laws.


Part V: Advanced topics - generalized coordinates

As I said before, in the preceeding example, our implementation of the Lagrangian is actually less efficient and more difficult to figure out than an implementation of Newtonian physics, since it involves calculating a potential, then, in essence, recalculating forces from it.

However, there are certain situations where it is very difficult to use Newton's laws, but the Lagrangian formalism may give us insight into how to solve the problem.

Take for example, the problem of a swinging pendulum.  We'll make the same approximations as we did before, namely that the gravitational force is constant with respect to height.  The pendulum, however, doesn't simply move up and down, it swings from side to side, and is constrained to move along a circular path.

The force on the pendulum isn't just the force of gravity, it's the force of the rod, pulling or pushing on the weight at the end, forcing it to keep it's (semi)circular path.  To get the forces to use in Newton's laws, we have to project gravity onto a vector parallel to the velocity of the pendulum, and go through several difficult step to get things to come out right.

Such a situation is a situation of constrained motion.  In this case, our weight at the end of a pendulum is constrained to a circular path.  In situations like this, the Lagrangian can provide very easy solutions to our problems.

There are two ways to get a solution to this problem with the Lagrangian.  One is to do exactly what we did before, except throw in something called "Lagrange multipliers".  This method is nice because it will also spit out the forces that could be used with Newton's laws.  However, it is not the simplest, or most elegant solution, so I won't cover it here.

The other way to solve this problem is to use something called generalized coordinates.

What are generalized coordinates?  They're exactly what the name says.  Who says that our coordinate system has to be cartesian?  Why do we have to use x, y, and z coordinates?  Why not something else?

The answer is, we can.  There's no reason to use x, y, and z as our coordinates, if there's something more efficient for our purposes.  And the best part is, the Lagrangian will still work!

Why use a different set of coordinates?  Well, let's think.  What do we know about this situation with the pendulum?

First of all, the z coordinate can't change, because the pendulum only swings back and forth along the x axis, moving up and down.  So let's throw out z (Even if we want the pendulum to be at an agle with respect to the x axis, it's only moving in one horizontal dimension, so we can simply rotate the vector of motion at the end).

Now that that's out of the way, let's think about what else we know.  Since the pendulum has to move on a circle, a natural choice of coordinates would be polar coordinates.  We can describe the location of our pendulum in terms of r, the radius and theta, the angle it makes with the vertical axis.

But we also know that the radius of our pendulum can't change, so we can throw out r.  Now all we're left with is theta, the angle our pendulum makes with the vertical axis.  We've just narrowed the number of coordinates we have from three down to one!

What we need to do now is calculate T and U in terms of our new coordinate.

In our coordinate system

X = R*sin(theta)
Y = -R*cos(theta)


Remember, R is a constant, not a coordinate.

When the angle theta with the vertical is zero, our pendulum is at the bottom, and when the angle is pi/2, our pendulum is facing straight off to the side.

You might have noticed that Y becomes negative when theta is zero.  We could write Y as

Y = -R*cos(theta) + c

where c is the height of the pendulum off the ground.  But there's really no point, since the constant c will just drop out of our equations anyway, since our potential U is linear with height.

To calculate kinetic energy, we need X'^2 + Y'^2.  Differentiating X and Y with respect to time, we get:

X' = R*Cos(theta)*theta' and
Y' = R*Sin(theta)*theta'

so:

X'^2 + Y'^2 = v^2 = R^2*(Cos(theta)^2 + Sin(theta)^2)*theta'^2

And by trig identities:

v^2 = R^2*theta'^2

So T and U, and L are:

T = 1/2*m*R^2*theta'^2
U = -m*g*R*Cos(theta)
L = T - U = m*(1/2*R^2*theta'^2 + Cos(theta)

Plugging this into our Lagrange equations, we get:

-m*Sin(theta) - m*R^2*theta'' = 0 or
Sin(theta) + R^2*theta'' = 0
Providing that I haven't made any mistakes Smiley

The equation Sin(theta) + R^2*theta'' = 0 is impossible to solve exactly by hand.  However, we can calculate it approximately, which means that our physics game engine could hand it just fine.

Here's a function that will update the location of the end of our pendulum as a function of time:

Code:
Function UpdatePend()

'R is a constant defined somewhere else in the program.
'Theta is a single, and ThetaAccel is a single that holds the change in ThetaVel.
'ThetaVel is a singel that holds the change in Theta with each cycle.
'Pos is a vector that holds position

ThetaAccel = Sin(theta)/R^2
ThetaVel = ThetaVel + ThetaAccel
Theta = Theta + ThetaVel

Pos.x = R*Sin(theta)
Pos.Y = -R*Cos(theta)

End Function


As you can see, we have come up with a relatively simple set of equations, and a very simple function to describe the motion of a pendulum from our somewhat ugly and daunting Lagrangian.  You can apply this method of using generalized coordinates to just about any situation.


Conclusion

Although it was originally my intention to provide an introduction to the Hamiltonian as well, I feel like I've gone on long enough and will leave that for another day.  Hopefully the introduction I've given to the Lagrangian will be of some use to some programmers out there.

Remember that the true usefulness of the Lagrangian doesn't come from the way it describes the equations of motion, but rather, from the way that it allows you to find the equations of motion.  It is a usefull tool in game physics for finding out how an object is supposed to act when it would be difficult to program in all the different forces that are acting on it.  For situations of constrained motion, the Lagrangian is a very useful tool.

I can see that it might be useful to apply the Lagrangian to a game that uses a landscape.  By constraining your land vehicle to move along the ground, the Lagrangian could describe the motion of the vehicle, even if it would be difficult to find the exact force that the ground exherts on the vehicle.  I won't try to figure out a good method to put in this tutorial, but rather, will leave it to you as an interesting problem to solve, if you would like.
Logged
darqSHADOW
Administrator
Community Member
*****
Posts: 2734


« Reply #13 on: April 17, 2005, 12:22:51 AM »

I think this should move into the wiki...  Wink

DS
Logged

TrueVision3D Project Manager
The fast and simple way of 3D development.
Technosites
Community Member
*
Posts: 102


« Reply #14 on: June 05, 2006, 12:02:29 PM »

Hi,
About your first tutorial.
Firstly, its excellent. Excactly want I wanted. You really should get a medal. But, how does that help me do movement (camera and ship) in space? I'm not quite sure what Wector Transformations are and I'm not sure if im already using them. I've read this is the best way of doing a space camera, so how do I do it. I don't want you to just do it all for me, an explaintion of how with code examples would be nice. (Visual Basic 6 if you can Smiley )

Many Thanks, if you can't help me more you've still helped a lot with ur tutorial.
Technosites
Logged

There are 10 types of people in this world. Those who understand binary and those who dont.
Pages: [1]
  Print  
 
Jump to:  

Powered by SMF 1.1.3 | SMF © 2006-2007, Simple Machines LLC
Seo4Smf v0.2 © Webmaster's Talks