User:Timmy Foxclaw/About Coordinate Systems and Rotations

About Coordinate Systems and Rotations

The four coordinate systems

There are four coordinate systems that are related to LSL programming:

1. World coordinates
2. Region coordinates
3. Object coordinates (root prim coordinates)
4. Child prim coordinates (coordinates relative to the prim that the script is in, or refered to by link number)

The World coordinates refer to the map, and allow to include the sim in the coordinates, or refer to void water.

Region coordinates are relative to a given sim. The origin is in the South/West corner at height 0. The North/East corner then is 256, 256 and a Z coordinate for the height up to 4096 meter (on opensim you can go even higher).

Object coordinates are relative to the root prim. Hence, if the object is moved or rotated then the orientation of a child prim, when given in object coordinates, doesn't change. In LSL "local position" and "local rotation" refer to this coordinate system. "local" means relative to the root prim.

Prim Coordinates are relative to a given prim. If this prim is the root prim then the Prim Coordinates are the same as the Object Coordinates. For example, if a child prim is a cube with a size of 1,1,1 and one red surface where the center of that surface is at 0.5,0,0 then it will still be at 0.5,0,0 no matter how you move or rotate that child prim (relative to the other linked prims).

Positions

A different position of the origin of a coordinate system is easy to understand: You can think of positions as vectors that start in the origin of the coordinate system that they are given in and end in the point that they refer to. While the length of the vector is independent of the rotation of the coordinate system, the three coordinates are not; but a mental picture of an arrow doesn't have little numbers for three coordinates, so that picture works independent of the rotation too.

Since the rotation of the World Coordinate system and the Region Coordinate system is the same (X, Y and Z axis are parallel of both), and since World Coordinates aren't used in many LSL functions to begin with, we will ignore World Coordinates for now and only refer to Region Coordinates, or say "global" when we mean Region Coordinates.

Rotations

An LSL rotation internally stores a vector that is the axis around which to rotate and the angle of the rotation around that axis. Let V = <u, v, w> be the normalized vector (that is, with length 1) around which we rotate and let α (alpha) be the angle around which we have to rotate. Then the LSL rotation is a quaternion stored as r = <x, y, z, s> = < V * sin(α/2), cos(α/2) >. Thus, r.x = x * sin(α/2), and r.s = cos(α/2) etc. Note that the quaternion is also normalized. Also note that there is a duality here because inverting V (making it point the opposite way) and inverting the angle gives the same rotation; r and -r have different values but are the same rotation. Also note that if you don't rotate at all (α == 0) then it doesn't matter what axis V you pick, which is apparent because V drops out since sin(0) = 0. The quaternion <0, 0, 0, 1> is the ZERO_ROTATION quaternion.

The point of this technical story is to show that for an LSL rotation to make sense in terms of orientation, you need to be able to express a vector in three coordinates (u, v, w above): the axis around which we rotate is expressed relative to the X-, Y- and Z-axes of the coordinate system. Hence, it is the orientation of the X-, Y- and Z-axes that defines the meaning of a rotation in LSL.

In terms of a mental picture the origin with the (orientation of the) three axis, the red X-axis, the green Y-axis and the blue Z-axis is all the reference we need, combined with a vector for position or a quaternion for rotation.

When you edit an object, the viewer shows either 'World' (meaning Region here) or 'Local' axes, but really the 'World' (region) axes show the wrong origin, shifted to an averaged center of the object, because if the origin was drawn at (0, 0, 0) you'd most likely not see it. The 'Local' ruler shows the correct coordinate system for the selected prim as its Prim Coordinate System. Selecting the root prim with 'Local' ruler on then shows the Object Coordinate System.

The dimension of rotations

As said before, given some coordinate system, any point in space can be represented with a vector. Obviously space is three dimensional, and thus vectors exists of three real values: one needs three distinct floating point numbers, the x coordinate, the y coordinate and the z coordinate to uniquely identify a position in a given coordinate system.

However, if one limits oneself to only normalized vectors, vectors with a length of one, then those represent all points on the surface of a sphere with radius 1. A surface is obviously two dimensional, so it should be possible to uniquely identify any point on the surface of such a sphere with only two floating point numbers.

One might think that selecting just two coordinates of the three of the vector will suffice because the third is fixed by the length requirement, but that only works for half spheres; for example, if x and y are known then z can still be either plus the square root of x squared plus y squared, or minus that value.

Instead, a better choice would be the spherical coordinate system and express the unit vectors with the two angular coordinates φ (phi) and θ (theta), where φ is the angle between the vector and the positive Z axis, and θ the angle that the projection of the vector on to the X,Y plane makes with the positive X axis. In other words, starting with a unit vector along the positive Z-axis, one can obtain the required point on the surface of the sphere by first rotating this vector around the Y axis (towards the positive X axis) by an angle of φ and then rotating the result around the Z-axis by an angle of θ. Note that in both cases the rotations are counter-clockwise when one looks at it from the positive side of the axis that one rotates around (towards the origin).

What we just did was expressing a unit vector in terms of two rotations: one around the Z axis and one around the Y axis; and indeed every rotation can be expressed as two rotations around those two axis (or really, any two arbitrary axes, as long as they are independent (not parallel)).

It is therefore possible to represent unit vectors (aka, a direction) with a rotation, as both are two dimensional. The mental picture here is that any unit vector can be expressed as some fixed unit vector (for example the one pointing along the positive Z axis, but any other would do as well) and the rotation needed to turn that 'starting vector' into the unit vector that one wants to represent. The other way around almost works as well, with the exception of the starting vector itself, as that can be expressed by any arbitrary rotation around itself, so that information was lost and it is not possible to know which of those rotations was used. In fact, any vector close to the starting vector would give inaccurate results in the light of floating point round off errors and is not a good way to represent a rotation. Of course, the two dimensional vector (θ, φ) would be excellent to represent both: the unit vector, as well as any rotation; but LSL stores vectors in x, y, z coordinates, not in spherical coordinates.

Converting between coordinate systems

Consider a child prim of an object that is not the root prim and express its position with the vector pos_oc, and its orientation with the rotation rot_oc, both relative to the Object Coordinates system.

A script inside that prim can easily obtain those values:

vector pos_oc = llGetLocalPos();
rotation rot_oc = llGetLocalRot();

The postfix '_oc' stands for Object Coordinates and is used to make clear relative to which coordinate system the x, y and z values of both variables are. Remember that also the rotation contains a vector (the axis around which we rotate), so from now on we'll talk about 'coordinates' relative to a coordinate system for both, the position as well as the rotation, where the coordinates refer to the x, y and z components of the position and the rotation. For example, the first sentence of this paragraph will read "...express its position with the vector pos_oc, and its orientation with the rotation rot_oc, both in Object Coordinates."

Note that one has to be careful here, because a script in the root prim of the object will not obtain Object Coordinates but instead Region Coordinates! The reason is that the 'Local' functions return coordinates relative to the 'parent' coordinate system, where the parent of a non-root prim is the root prim, but the parent of the root prim is the region (or avatar in case of an attachment). The position and rotation of the root prim in Object Coordinates is trivial (ZERO_VECTOR and ZERO_ROTATION respectively).

Back to the script in the child prim. It can also easily obtain its position and rotation in Region Coordinates:

vector pos_rc = llGetPos();
rotation rot_rc = llGetRot();

Both sets, (pos_oc, rot_oc) and (pos_rc, rot_rc) refer to the same thing, so it should be possible to convert them into each other.

One can do this by looking at the orientation (position and rotation) of one coordinate system relative to the other and expressed in the coordinates of that other. The position and rotation of the Object Coordinate System relative to the Region Coordinate System, and in Region Coordinates, could be obtained as follows:

vector oc_pos_rc = llGetRootPosition();
rotation oc_rot_rc = llGetRootRotation();

As indicated by the '_rc' postfix, these are again in Region Coordinates.

Having the relative orientation of the Object Coordinate System in Region Coordinates, it is possible to convert Object Coordinates into Region Coordinates. More abstractly put, if you have two coordinate systems A and B, and you have their relative orientation expressed in A then you can convert orientations expressed in B into A.

The formula is as follows:

vector   pos_A = pos_B * B_rot_A + B_pos_A;            (1)
rotation rot_A = rot_B * B_rot_A;

The order of the multiplications is important; swapping the B_rot_A to the front will not work!

Note that the B_*_A are the position and rotation of the coordinate system B expressed in terms of the coordinate system A. The reason that I use the prefix and postfix in this order is directly related to the order enforced for rotation multiplications. Written this way, the left-hand postfix of a multiplication must match the right-hand prefix and conviently short-circuits the notation:

something_A = something_B * B_rot_A

the two B's 'connect' and drop out (along with the '_rot_') to give something_ ... A. This only works when the postfix (_A or _B) refers to the coordinate system that the variables are in, and the prefix is whatever is described by the variable, possibly another coordinate system (B_ in this example).

A complete example script follows:

default
{
    touch(integer num_detected)
    {
        llSay(0, "---------");

        // Get the pos and rot in Object Coordinates.
        vector pos_oc = llGetLocalPos();
        rotation rot_oc = llGetLocalRot();
        
        // Fix these values when this script is in the root prim!
        if (llGetLinkNumber() <= 1)
        {
            pos_oc = ZERO_VECTOR;
            rot_oc = ZERO_ROTATION;
        }
        
        // Get the pos and rot in Region Coordinates.
        vector pos_rc = llGetPos();
        rotation rot_rc = llGetRot();
        
        // Get the pos and rot of the Object Coordinate System in Region Coordinates.
        vector oc_pos_rc = llGetRootPosition();
        rotation oc_rot_rc = llGetRootRotation();
        
        // Print the values of the positions (the rotations are too hard to understand).
        llSay(0, "oc_pos_rc = " + (string)oc_pos_rc);
        llSay(0, "pos_os = " + (string)pos_os);
        llSay(0, "pos_rc = " + (string)pos_rc);
        
        // Calculate pos_rc and rot_rc from pos_oc and rot_oc : converting between _oc and _rc.
        vector   pos2_rc = pos_oc * oc_rot_rc + oc_pos_rc;
        rotation rot2_rc = rot_oc * oc_rot_rc;
        
        // Print the result for the position, and show the difference with what llGetPos() and llGetRot() returned.
        llSay(0, "pos2_rc = " + (string)pos2_rc);
        // These values should be ZERO_VECTOR and ZERO_ROTATION (or -ZERO_ROTATION).
        llSay(0, "pos_rc - pos2_rc = " + (string)(pos_rc - pos2_rc));
        llSay(0, "rot_rc / rot2_rc = " + (string)(rot_rc / rot2_rc));
    }
}

Obviously there is no preference for one coordinate system over the other for this conversion; one can equally as well convert orientations in A into B:

vector   pos_B = pos_A * A_rot_B + A_pos_B;
rotation rot_B = rot_A * A_rot_B;

The question arises then how to obtain the A_*_B values; there is no LSL function that returns directly the position or rotation of the region relative to a prim (aka, rc_*_oc).

Intuitively one would say that it has to be possible to express A_*_B in B_*_A on grounds of symmetry, and that is indeed the case.

To show how this works, lets just take the formula for pos_A and rot_A (see (1) above) and rework them to get pos_B and rot_B on the left-hand side. This results first in

pos_B * B_rot_A = pos_A - B_pos_A
rot_B * B_rot_A = rot_A

Then divide both sides by B_rot_A to get:

pos_B = pos_A / B_rot_A - B_pos_A / B_rot_A
rot_B = rot_A / B_rot_A

from which we can conclude that

A_pos_B = - B_pos_A / B_rot_A                           (2)
A_rot_B = ZERO_ROTATION / B_rot_A

Note how A_rot_B is the inverse of B_rot_A: A_rot_B * B_rot_A = ZERO_ROTATION. As per our notation convention you'd expect A_rot_B * B_rot_A to be A_rot_A, and that is actually the case as A_rot_A means the relative orientation of coordinate system A expressed in A's coordinates, which is the trivial ZERO_VECTOR and ZERO_ROTATION as we saw before for the root prim orientation expressed in Object Coordinates.

And for fun, note that if we do this inversion again we get:

B_pos_A = - A_pos_B / A_rot_B = - (- B_pos_A / B_rot_A) / (ZERO_ROTATION / B_rot_A) = B_pos_A / B_rot_A * B_rot_A = B_pos_A
B_rot_A = ZERO_ROTATION / A_rot_B = ZERO_ROTATION / (ZERO_ROTATION / B_rot_A) = B_rot_A

as expected.

We can now also convert Region Coordinates to Object Coordinates.

One could first calculate,

vector rc_pos_oc = - llGetRootPosition() / llGetRootRotation();
rotation rc_rot_oc = ZERO_ROTATION / llGetRootRotation();

and then convert as usual:

vector   pos_oc = pos_rc * rc_rot_oc + rc_pos_oc;
rotation rot_oc = rot_rc * rc_rot_oc;

or if you don't want to calculate the intermediate rc_*_oc values, you could do immediately:

vector   pos_oc = (pos_rc - llGetRootPosition()) / llGetRootRotation();
rotation rot_oc = rot_rc / llGetRootRotation();

To verify this, you could extend the above script with

        // Calculate pos_oc and rot_oc from pos_rc and rot_rc.
        vector   pos2_oc = (pos_rc - llGetRootPosition()) / llGetRootRotation();
        rotation rot2_oc = rot_rc / llGetRootRotation();
        
        llSay(0, "pos2_oc = " + (string)pos2_oc);
        // These values should again be ZERO_VECTOR and +/- ZERO_ROTATION.
        llSay(0, "pos_oc - pos2_oc = " + (string)(pos_oc - pos2_oc));
        llSay(0, "rot_oc / rot2_oc = " + (string)(rot_oc / rot2_oc));

More in general, if you want to convert from B to A, but you only have A_*_B instead of B_*_A, then (combining (1) and (2)) you can do:

vector   pos_A = (pos_B - A_pos_B) / A_rot_B;           (3)
rotation rot_A = rot_B / A_rot_B;

Cascading coordinate systems

Imagine you have a castle object where the floor is the root prim. Everything is linked because you want it to be easy to move the castle around. One of the child prims is a sculpty representing the hinges around which you want to rotate a large gate. The gate has a rotating wheel on it.

You want that the object (castle) to keep working when it is moved and/or rotated, but also when it is edited and the hinges are moved and/or rotated. This means that the state of the gate prims must be stored relative to the hinges, which is by far the easiest thing to do to begin with, after all, the gate rotates around the hinges and has no relation with the rest of the castle.

Hence, we have the Region Coordinate System, the Object Coordinate System (the castle) and the Hinges Coordinate System. The latter is our first example of a Prim Coordinates System, but I'll use the postfix '_hc' for "Hinges Coordinates" instead of '_pc'. Finally we have the Spill Coordinate System (or Gate Coordinates) where the spill is a prim around which we want to rotate the wheel, and I'll use '_sc' as postfix for that.

The variables used then would be an orientation for the wheel in *_sc coordinates, the orientation of the spill, and other gate prims, in *_hc coordinates and the orientation of the hinges in *_oc coordinates. Finally we could have variables for the orientation of the castle (root prim) in *_rc coordinates of course.

vector   castle_pos_rc; // The position of the castle in Region Coordinates.
rotation castle_rot_rc; // The rotation of the castle in Region Coordinates.

vector   hinges_pos_oc; // The position of the hinges in Object Coordinates.
rotation hinges_rot_oc; // The rotation of the hinges in Object Coordinates.

vector   spill_pos_hc;  // The position of the spill in Hinges Coordinates.
rotation spill_rot_hc;  // The rotation of the spill in Hinges Coordinates.

vector   wheel_pos_sc;  // The position of the wheel in Spill Coordinates.
rotation wheel_rot_sc;  // The rotation of the wheel in Spill Coordinates.

Now because we used existing prims orientations as coordinate system for the next prim, the position and rotation of those prims are the position and rotation of that coordinate system. Just like before we saw that llGetPos/llGetRot returns the global orientation of the root prim when used in a script in the root prim, while llGetRootPos/llGetRootRot gave us the relative orientation of the Object Coordinate System in Region Coordinates, but is in fact the same as what llGetPos/llGetRot in the root prim returns: the position/rotation of the root prim is the orientation of the Object Coordinate System relative to the Region.

Likewise, in this case, castle_*_rc is the relative orientation of the Object Coordinate System in Region Coordinates as well as the orientation of the root prim in Region Coordinates; they are the same thing.

The same then applies to the *_oc variables, which we defined as the orientation of the hinges in Object Coordinates, but at the same time are the relative orientation of the Hinges Coordinate System in Object Coordinates. The *_hc variables are the orientation of spill in Hinges Coordinates, but also the relative orientation of the Spill Coordinates in Hinges Coordinates.

Lets write it out using the same notation convention as used above, then we have:

vector   oc_pos_rc = castle_pos_rc; // The Object Coordinate System orientation in Region Coordinates.
rotation oc_rot_rc = castle_rot_rc;

vector   hc_pos_oc = hinges_pos_oc; // The Hinges Coordinate System orientation in Object Coordinates.
rotation hc_rot_oc = hinges_rot_oc;

vector   sc_pos_hc = spill_pos_hc;  // The Spill Coordinate System orientation in Hinges Coordinates.
rotation sc_rot_hc = spill_rot_hc;

It will therefore probably not surprise you that

vector   hc_pos_rc = hc_pos_oc * oc_rot_rc + oc_pos_rc; // The Hinges Coordinate System orientation in Region Coordinates.
rotation hc_rot_rc = hc_rot_oc * oc_rot_rc;

vector   sc_pos_rc = sc_pos_hc * hc_rot_rc + hc_pos_rc; // The Spill Coordinate System orientation in Region Coordinates.
rotation sc_rot_rc = sc_rot_hc * hc_rot_rc;

And finally, using the latter, we can express the wheel orientation in region coordinates:

vector   wheel_pos_rc = wheel_pos_sc * sc_rot_rc + sc_pos_rc;
rotation wheel_rot_rc = wheel_rot_sc * sc_rot_rc;

Although this was nothing else but applying the coordinate system conversion formula ((1)) three times, it might give some insight to get rid of the intermediate variables and express the wheel orientation in region coordinates without them, giving:

wheel_pos_rc = ((wheel_pos_sc * sc_rot_hc + sc_pos_hc) * hc_rot_oc + hc_pos_oc) * oc_rot_rc + oc_pos_rc
wheel_rot_rc = wheel_rot_sc * sc_rot_hc * hc_rot_oc * oc_rot_rc

Applying rotations and translations

Applying rotations and translations (moving) is very much the same as coordinate system transformations described above. To show this consider the following: imagine you have a prim 'pa' with its own coordinate system A. Then it trivially has position ZERO_VECTOR and rotation ZERO_ROTATION relative to A, per definition.

vector   pa_pos_A = ZERO_VECTOR;
rotation pa_rot_A = ZERO_ROTATION;

Now consider a prim 'pb' that has some relative position pb_pos_A and rotation pb_rot_A given in A coordinates. For example,

vector   pb_pos_A = <0.5, 0.6, 0.7>;
rotation pb_rot_A = <0.38, 0.26, 0.34, 0.82>;

Note that 0.38² + 0.26² + 0.34² + 0.82² = 1, as it must be since rotations are normalized.

As seen before, the orientation of pb relative to pa is the same as the orientation of pb relative to A, and the same as the orientation of B relative to A. pa and A are interchangeable as long as pa has ZERO_VECTOR and ZERO_ROTATION relative to A. Thus, with B being the coordinate system where pb currently is, we also have:

vector   B_pos_A = pb_pos_A; // The orientation of coordinate system B, in coordinates of A.
rotation B_rot_A = pb_rot_A;

Now suppose we want move the prim pa (without moving the coordinate system A), by applying the translation 'B_pos_A', after all a translation is stored in a vector and B_pos_A is a vector. Likewise, we want to rotate pa (again without influencing A) by applying the rotation B_rot_A. Then the most logical order to this in is by first applying the rotation and then applying the translation, because then the meaning of the translation is preserved and will match the meaning of B_pos_A in that pa will end up where pb is now. If you'd first apply the translation so that pa ends up where pb is, and then apply the rotation (around A!) then pa would swoop around the origin of A and move away from its desired position.

Doing this correctly should leave us with pa_pos_A == pb_pos_A and pa_rot_A == pb_rot_A, so that pa ended up precisely where pb is now. We can't just assign those values however, because we're looking for a formula that works in general; a way to apply B_rot_A and then B_pos_A to any prim with -say- orientation pc_pos_A, pc_rot_A.

Thus, first we apply the rotation:

pc_rot_A *= B_rot_A;
pc_pos_A *= B_rot_A;

And then apply the translation:

pc_pos_A += B_pos_A;

Now lets check if indeed this would cause pa to end up where pb is. The effect on pa_*_A would be:

pa_rot_A = pa_rot_A * B_rot_A = ZERO_ROTATION * B_rot_A = B_rot_A = pb_rot_A;
pa_pos_A = pa_pos_A * B_rot_A = ZERO_VECTOR * B_rot_A = ZERO_VECTOR;

and then the translation

pa_pos_A = pa_pos_A + B_pos_A = ZERO_VECTOR + B_pos_A = B_pos_A = pb_pos_A;

So yes, pa_*_A ends up as pb_*_A!

The real test however would be if we first applied some transformation to pa to bring it to some arbitrary orientation (aka, pc), then apply the transformation B_rot_A,B_pos_A and then reverse the first transformation but now relative to B, and then find that we indeed end up at pb.

Let us first define a 'transformation', which exists of a vector for the translation and a rotation for the rotation. We will have two transformations during this test: t1: pa -> pc, and t2: pa -> pb.

As we just saw (or we hope this is going to be the case), t2 can be defined as:

vector   t2_translation = B_pos_A;
rotation t2_rotation = B_rot_A;

where as t1 is arbitrary, just some values t1_translation and t1_rotation.

Applying t1 to pa, we get (this is thus the same as above, but now with the 't1' notation):

pa_rot_A *= t1_rotation;
pa_pos_A *= t1_rotation;
pa_pos_A += t1_translation;

Next we apply translation t2 to that:

pa_rot_A *= t2_rotation;
pa_pos_A *= t2_rotation;
pa_pos_A += t2_translation;

Finally we want to reverse t1, but in the coordinate system of B. Therefore, we first convert our coordinates to that of B (formula (1)):

pa_pos_B = pa_pos_A * A_rot_B + A_pos_B;
pa_rot_B = pa_rot_A * A_rot_B;

where (formula (2))

A_pos_B = - B_pos_A / B_rot_A;
A_rot_B = ZERO_ROTATION / B_rot_A;

and then reverse t1. Note now we have to do the translation first!

pa_pos_B -= t1_translation;
pa_pos_B /= t1_rotation;
pa_rot_B /= t1_rotation;

And if everything worked out then pa now should be at ZERO_VECTOR, ZERO_ROTATION relative to B!

As pa_pos_A,pa_rot_A started as ZERO_VECTOR,ZERO_ROTATION, the first step of applying t1 resulted in pa_rot_A = t1_rotation and pa_pos_A = t1_translation.

After applying t2 we get pa_rot_A = t1_rotation * t2_rotation, and pa_pos_A = t1_translation * t2_rotation + t2_translation.

Converting that to B gives

pa_pos_B = (t1_translation * t2_rotation + t2_translation) * A_rot_B + A_pos_B
pa_rot_B = (t1_rotation * t2_rotation) * A_rot_B

where, remembering that t2_translation = B_pos_A and t2_rotation = B_rot_A,

A_pos_B = - B_pos_A / B_rot_A = - t2_translation / t2_rotation
A_rot_B = ZERO_ROTATION / B_rot_A = ZERO_ROTATION / t2_rotation

so that we get

pa_pos_B = (t1_translation * t2_rotation + t2_translation) / t2_rotation - t2_translation / t2_rotation = t1_translation
pa_rot_B = (t1_rotation * t2_rotation) / t2_rotation = t1_rotation

and finally reversing t1 gives:

pa_pos_B = (t1_translation - t1_translation) / t1_rotation = ZERO_VECTOR
pa_rot_B = t1_rotation / t1_rotation = ZERO_ROTATION

So that we can conclude that if you have two prims pa and pb, each representing their own coordinate system A and B respectively, and you express the orientation of B relative to A as the pair (B_pos_A, B_rot_A) — then the following holds.

The transformation pair (t_translation, t_rotation) equals (B_pos_A, B_rot_A) as long as you first apply the rotation and then the translation while working in the coordinate system of A.

The latter is also very practical, because you normally won't know anything about the target coordinate system B until you actually did the transformation.

Now remember that we found that in order to convert variables from being relative to a coordinate system B to a coordinate system A, we did (formula (1)):

pos_A = pos_B * B_rot_A + B_pos_A
rot_A = rot_B * B_rot_A

which is thus exactly the same as first applying the rotation t1_rotation and then the translation t1_translation on the pair (pos_B, rot_B) where that pair is relative to A! That isn't too weird because if you consider a prim in the origin of A with no rotation, aka (ZERO_VECTOR, ZERO_ROTATION) relative to A and then apply the transformation then you expect to end up in the origin of B with no rotation relative to B. Aka pos_B = ZERO_VECTOR and rot_B = ZERO_ROTATION. So now work backwards: assume you already have the prim in B like that, then what are its coordinates in A? That would be the above conversion thus, with pos_B = ZERO_VECTOR and rot_B = ZERO_ROTATION.

Hence, doing the conversion with ZERO as input gives the same results as doing the t transformation with ZERO as input! And the final formula for a transformation just looks exactly the same:

target_pos_A = source_pos_A * t_rotation + t_translation        (4)
target_rot_A = source_rot_A * t_rotation

Scaling coordinate systems

Apart from rotation and translation, there is one more linear operator that can be used to convert between coordinate systems: scaling.

Let A_pos_B be the position of the origin of A in the coordinates of B. Let A_rot_B be the orientation of A within B, and let the vector A_B_scale be the scaling factors for each axes. Scaling has no influence on the origin of a coordinate system, but it matters for the rotation when the scaling factors per axis are different.

I think that a good choice is scale first and rotate afterwards. Because that leaves the values of the coordinates of A_rot_B alone: it remains the rotation as defined in B, which makes sense since it is written in the coordinates of B. However, it's just a matter of choice really, which in this case will deform B_rot_A instead of A_rot_B. I can imagine that the choice might depend on the application.

Chosing that we first scale, a point at pos_A is mapped onto the point pos_B as follows:

pos_B = (S * pos_A) * A_rot_B + A_pos_B

where S is a 3x3 matrix with as diagonal the vector A_B_scale (and the rest zeroes). Note that when each element of A_B_scale is the same then you can just use a scalar for S, equal to the scale of each axis.

Whether or not S is a matrix or a scalar, lets just replace it with 'A_B_scale' here. Just keep in mind that if it isn't a scalar then this won't work in LSL directly: then you'll have to split it up in three different code lines because LSL doesn't support scaling or matrices.

 pos_B = pos_A;
 pos_B.x *= A_B_scale.x;  // Scale
 pos_B.y *= A_B_scale.y;
 pos_B.z *= A_B_scale.z;
 pos_B *= A_rot_B;        // Rotate
 pos_B += A_pos_B;        // Translate

Note that the inverse of S, 1/S, is a 3x3 matrix with as diagonal the inverse of each element of A_B_scale. We will either write 'B_A_scale', or divide by 'A_B_scale' in our notation.

Therefore our notation becomes:

pos_B = ('A_B_scale' * pos_A) * A_rot_B + A_pos_B

Since A and B are arbitrary, we would expect the same formula to hold when swapping A and B:

pos_A = ('B_A_scale' * pos_B) * B_rot_A + B_pos_A

Solving pos_B from that by subtracting B_pos_A from both sides, then dividing both sides by B_rot_A and finally left multiplying both sides with the inverse of B_A_scale, gives:

pos_B = 'A_B_scale' * ((pos_A - B_pos_A) / B_rot_A) = 'A_B_scale' * (pos_A / B_rot_A) - 'A_B_scale' * (B_pos_A / B_rot_A)

from which we can conclude (by setting pos_A to ZERO_VECTOR) that

A_pos_B = -'A_B_scale' * (B_pos_A / B_rot_A)

Solving B_pos_A from that by left multiplying both sides with -B_A_scale and then right multiplying both sides with B_rot_A, gives:

B_pos_A = - ('B_A_scale' * A_pos_B) * B_rot_A

By setting A_pos_B to ZERO_VECTOR, and this B_pos_A too, we can also conclude that

('A_B_scale' * pos_A) * A_rot_B = 'A_B_scale' * (pos_A / B_rot_A)

If 'A_B_scale' is a scalar, so that the order of multiplication doesn't matter, then the scale just cancels on both sides and we find the conversion that we found before, that B_rot_A is the inverse of A_rot_B. However, if A_B_scale is a diagonal matrix then the order matters and the brackets are significant.

There is no way to let pos_A drop out: we cannot write B_rot_A as function of A_B_scale and A_rot_B.

The conclusion is that a scale is asymmetrical: if you choose that a scale A_B_scale means that you first apply the scale on a vector in A and then rotate, then that automatically implies that when going from B to A you first have to rotate and then apply the scale.

To make the choice clear in the name, I suggest to put an 'r' in front of the coordinate system that need rotation to be applied first.

In summary, if A_rB_scale is defined as

pos_B = ('A_rB_scale' * pos_A) * A_rot_B + A_pos_B      (5)

then that implies that rB_A_scale has to be applied after rotation. Note this means braces, NOT changing the left- or right- side of multiplications! Matrices (scales) are always on the left side of a vector, while rotations are always on the right side. Thus we get:

pos_A = 'rB_A_scale' * (pos_B * B_rot_A) + B_pos_A      (6)

If now we solve pos_B from the latter, we get:

pos_B = ('A_rB_scale' * (pos_A - B_pos_A)) / B_rot_A = ('A_rB_scale' * pos_A) / B_rot_A - ('A_rB_scale' * B_pos_A) / B_rot_A

Hence

B_pos_A = - 'rB_A_scale' * (A_pos_B / A_rot_B)          (7)
B_rot_A = ZERO_ROTATION / A_rot_B

and

A_pos_B = - ('A_rB_scale' * B_posA) / B_rot_A           (8)
A_rot_B = ZERO_ROTATION / B_rot_A

Note how these two sets condense into the same thing when the scale is a scalar, but also when A_rot_B and B_rot_A are ZERO_ROTATION. In fact, when the rotation involved is of the special case of one or more rotations of 90 degrees around any of the axes, so that the rotation merely swaps axes, one could define the inverse of A_B_scale by applying that rotation to that vector (swapping the scaling factors), next to inverting each factor individually. That is a different way of how we defined B_A_scale, but it would collapse the two sets into one for that special set of rotations. I'm using that trick in the example script below.

Direction vectors

A 'direction' is a normalized vector (its length is 1), or unit vector, where only the direction it points to is important.

Of course, in the light of rotations between coordinate systems, also the directions are relative to some coordinate system.

As stated before, using some 'starting' unit vector a rotation can be expressed as a direction by applying the rotation to a starting vector. Reversing this process runs into floating point round off errors however, or might even have lost all information of the rotation when we rotate around the starting vector itself.

There is a better way to convert directions (back) to a rotation however: by using two directions that are perpendicular, each having a unit vector along one of the axes as starting vector. That way there is always at least one vector enough rotated away from the actual rotation axis in order not to suffer from (too much) information loss.

For example, let us use the unit vector along the positive X axis and the unit vector along the positive Z axis. In LSL those are often called fwd and up respectively. Furthermore, lets work in a coordinate system A. Then we have the 'starting' vectors:

vector A_fwd_A  = <1, 0, 0>;            (9)
vector A_left_A = <0, 1, 0>;
vector A_up_A   = <0, 0, 1>;

We said we only knew the fwd and up directions, but the third can always be calculated from the other two with a cross product:

A_fwd_B  = A_left_B % A_up_B;           (10)
A_left_B = A_up_B % A_fwd_B;
A_up_B   = A_fwd_B % A_left_B;

Hence, the cross product exists of the other two in the order left to right where after 'up' you start at 'fwd' again; they rotate: ... left up fwd left up fwd ... Put a '=' at any space in that sequence and a '%' at the next space and you have the right formula. Also note that x % y = - y % x, in case you were wondering.

Now assume we have the target directions A_*_B, then we can find the rotation between A and B with:

A_rot_B = llAxes2Rot(A_fwd_B, A_left_B, A_up_B);       (11)

Where A_rot_B is the rotation that turns A_*_A into A_*_B:

A_fwd_B  = A_fwd_A  * A_rot_B
A_left_B = A_left_A * A_rot_B
A_up_B   = A_up_A   * A_rot_B

Finding coordinate systems

Assume you already have a coordinate system B (for example, region or object coordinates, the orientations of which you can find directly with #LSL functions).

Finding coordinate system A means that you find A_rot_B (or B_rot_A which is just the inverse), and A_pos_B or B_pos_A.

Finding the rotation is usually the easiest since it doesn't depend on where the origin is and can most likely be derived somehow. For example, in many cases you will have two orthogonal direction vectors in B that you know to be parallel with the x, y or z axes of A. Lets say that you have AsX_dir_B and a AsY_dir_B that make an angle of 90 degrees and have a length of 1, where AsX_dir_B is known to be parallel with the X axis of A and AsY_dir_B is known to be parallel with the Y axis of A. Note that in this case we can also use the notation A_fwd_B for AsX_dir_B since 'fwd' is in the X direction (AsX_dir_B means "A's X-axis direction in coordinates of B", and A_fwd_B means "A's forwards direction in coordinates of B", so they are the same thing). Then we find:

A_rot_B = llAxes2Rot(AsX_dir_B, AsY_dir_B, AsX_dir_B % AsY_dir_B);

as discussed in the previous paragraph.

Next, if you have a vector from_B_to_A pointing from the origin of B to the origin of A, then you have the translation between A and B too:

If from_B_to_A is in coordinates of A, then

B_pos_A = -from_B_to_A

while if from_B_to_A is in coordinates of B, then

A_pos_B = from_B_to_A

Slight counter intuitive you have to read 'A_pos_B' as 'The vector pointing from the origin of B to the origin of A', but plausible since it is in the coordinates of B and thus you should expect the vector to start in the origin of B.

If you have a series of vectors that are stacked, all the way from B to A, then of course you just add them all up to get A_pos_B. Likewise, if you have a series of vectors in the coordinates of A that are stacked to add up to go from A to B, then add them all up to get B_pos_A.

The interesting case therefore is where you have a vector C_pos_A that points from A to C and is in the coordinates of A, and a vector C_pos_B that points from B to C and is in the coordinates of B.

To solve this puzzle, lets first write the vectors slightly different so we're not confused about what they mean, even after converting between the coordinate systems A and B.

Let AC_vec_A be C_pos_A: the vector (arrow) from the origin of A to C (in coordinates of A) and let BC_vec_B be C_pos_B: the vector from the origin of B to C (in the coordinates of B). Note how the notation of _vec_ is similar to that of _dir_, except that a direction vector has a length of 1, and this vector has both a direction and a length: it describes how to get from one point to another. Like directions, when converting them between A and B you ONLY apply the rotation A_rot_B! Also note that a vector with a minus sign just means "in the opposite direction", hence CA_vec_A = -AC_vec_A.

Thus, we find that

AC_vec_B = AC_vec_A * A_rot_B

and find

A_pos_B = BA_vec_B = BC_vec_B + CA_vec_B = BC_vec_B - AC_vec_B = BC_vec_B - AC_vec_A * A_rot_B

then dropping the extra notation again

A_pos_B = C_pos_B - C_pos_A * A_rot_B         (12)

Which is rather logical when you read it in English as "To go from B to A in coordinates of B (A_pos_B), you first go from B to C in coordinates of B (C_pos_B) and then in the opposite direction (the minus sign) of going from A to C (namely from C to A) converted to the coordinates of B.

And oh wonder - if we write this result a little different we find:

C_pos_A = (C_pos_B - A_pos_B) / A_rot_B

which is the normal conversion (formula (3)) for the position of C from A to B coordinates!

LSL functions

What all of the above taught you is mainly a coding style: use prefixes and postfixes for your rotations (and translations)! It's all in the names of your variables. A 'rotation' type can be the rotation part of a transformation (not tied to a particular coordinate system, but rather rotating around whatever coordinate system you apply it to), it can be the representation of a unit vector, it can be the rotation part of the orientation of a prim relative to a given coordinate system, or it can represent the rotation part of a coordinate system conversion. If you don't use a consistent way to reflect all that in your variable names then you will get confused.

The LSL functions mostly deal with two types: orientations (position + rotation) of prims relative to some coordinate system, and simply transformations (rotation and/or translation).

The table below shows the coordinate systems involved.

		Coordinates of the prim containing the script (sc_)
Coordinate System	Abbreviation	script in child prim	script in root prim
World/Map Coordinates	wc	sc_pos_wc = sc_pos_rc + llGetRegionCorner() sc_rot_wc = sc_rot_rc
Region Coordinates	rc	sc_pos_rc = llGetPos() sc_rot_rc = llGetRot()
Object/Root Coordinates	oc	sc_pos_oc = llGetLocalPos() sc_rot_oc = llGetLocalRot()	sc_pos_oc = ZERO_VECTOR sc_rot_oc = ZERO_ROTATION
Self Coordinates	sc	sc_pos_sc = ZERO_VECTOR sc_rot_sc = ZERO_ROTATION

Where 'in the root prim' can be detect by testing that llGetLinkNumber() returns a value <= 1.

Alternatively, you can set link = llGetLinkNumber() and use the following table which shows how to obtain the coordinates of a prim with link number 'link' from a script anywhere in the same object. Note that this doesn't work when the object exist of a single prim: then llGetLinkNumber() returns 0, while 'link' must be larger than 0 here. If you need to get Object Coordinates (_oc) then do not use the negative value LINK_THIS. You can use LINK_THIS for the other coordinate systems, provided the object exists of at least two prims.

		Coordinates of a linked prim (link_)
Coordinate System	Abbreviation	link > 1	link = 1
World/Map Coordinates	wc	link_pos_wc = link_pos_rc + llGetRegionCorner() link_rot_wc = link_rot_rc
Region Coordinates	rc	link_pos_rc = llList2Vector(llGetLinkPrimitiveParams(link, [PRIM_POSITION]), 0) link_rot_rc = llList2Rot(llGetLinkPrimitiveParams(link, [PRIM_ROTATION]), 0)
Object/Root Coordinates	oc	link_pos_oc = llList2Vector(llGetLinkPrimitiveParams(link, [PRIM_POS_LOCAL]), 0) link_rot_oc = llList2Rot(llGetLinkPrimitiveParams(link, [PRIM_ROT_LOCAL]), 0)	link_pos_oc = ZERO_VECTOR link_rot_oc = ZERO_ROTATION
Self Coordinates	sc	link_pos_sc = (link_pos_oc - sc_pos_oc) / sc_rot_oc link_rot_sc = link_rot_oc / sc_rot_oc

While the above table works for link = LINK_ROOT, it is probably more convenient to use the following table:

		Coordinates of the root prim (oc_)
Coordinate System	Abbreviation	script in any prim
World/Map Coordinates	wc	oc_pos_wc = oc_pos_rc + llGetRegionCorner() oc_rot_wc = oc_rot_rc
Region Coordinates	rc	oc_pos_rc = llGetRootPosition() oc_rot_rc = llGetRootRotation()
Object/Root Coordinates	oc	oc_pos_oc = ZERO_VECTOR oc_rot_oc = ZERO_ROTATION
Self Coordinates	sc	oc_pos_sc = - sc_pos_oc / sc_rot_oc oc_rot_sc = ZERO_ROTATION / sc_rot_oc

As shows from the above table, LSL does not directly support _sc coordinates, and I used conversions to fill in the tables. All in all, we listed four coordinate systems (wc, rc, oc and sc; not including the 'link_' stuff, which isn't really a separate coordinate system) and thus could write out four times four is sixteen ways to convert one into the other. Four of those would be trivial, and of the remaining twelve we can express six as the inverse of the other six.

However, lets ignore the World Coordinates for now, which are just the same as the Region Coordinates with an offset for the position. Then we have three coordinate systems (rc, oc and sc), leading to nine conversion, three of which are trivial, and of the remaining six three can be expressed as the inverse of the other three. Those three are: oc_*_rc, sc_*_rc and sc_*_oc, which are all listed above. The trivial ones, rc_*_rc, oc_*_oc and sc_*_sc, are all ZERO_* (the latter two also are listed in the tables above). Finally, the inverse values, rc_*_oc, rc_*_sc and oc_*_sc can be expressed with the inversion formula given in #Converting between coordinate systems. So this is really all there is.

The fourth dimension

Picture a coordinate system as three perpendicular lines, like an object in-world perhaps existing of three really long prims in the colors red, green and blue. You can toss them into the region where they randomly bounce around like dice – a whole bunch of them! And then you freeze time - make a snapshot.

With the knowledge so far you can now express the position and rotation of an object or prim relative to any of those coordinate systems – and convert between them at will. Apply rotations and translations relative to any chosen system. All of the coordinate systems are virtually the same! But in practice this isn't the case, why is that?

The reason is that in practice coordinate systems change in time, in a way that you cannot predict. It's like they are still tumbling and bouncing around and you have no idea where they are.

This noise comes mainly from one thing: objects are being moved. Taken into repositories and rezzed back, rotated and possibly even scaled! In fact, not just objects, but also prims inside objects can be changed at will at any moment and a good script must take that into account.

The art of programming therefore becomes design. You must choose which coordinate systems you want to use; what your global variables store, relative to which coordinate system.

Unfortunately, a script can not detect when this might have happened and only perform extra computations or re-computation, of otherwise constant variables, when a prim or object has actually been changed. There is an LSL event moving_start but that only works half on SL and not at all on opensim. For rotation detection there is nothing.

As such, it is necessary to call the LSL functions listed in the tables above every time you want to do anything that depends on them; which hopefully won't be too often.

The following example script can be put in the root prim of an object of any number of cubes and allows one to sit on each cube, using llLinkSitTarget, even when you move or rotate child prims.

// A script written to show how to use coordinate systems in LSL.
//
// Copyright(c) Timmy Songbird @DreamNation
//              Timmy Foxclaw  @SL,
//              2014
//
// This program is free software; you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published by
// the Free Software Foundation; either version 2 of the License, or
// (at your option) any later version.
//
// This program is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
// GNU General Public License for more details.
//
// You should have received a copy of the GNU General Public License
// along with this program; if not, write to the Free Software
// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307, USA.')

//----------------------------------------------------------------------------
//#line "/opt/secondlife/scripts/scripts/tests/llLinkSitTarget_Test.m4s:5"
// GitLocation: tests/llLinkSitTarget_Test.m4s
// GitID: 32c3ab5a6f17d96eafe4e110bb91a700e2fd589b

// This script uses the follow coordinate systems:
//
// rc     : Region Coordinates.
// oc     : Object, or root prim, Coordinates.
// link   : Coordinates of prim with link number 'link'.
// STc    : Coordinates of the touched face (see llDetectedTouchST)
//          (x and y in the ranger [0, 1] with (0, 0) in the bottom-left corner).
// center : Coordinates of the touched face with the origin in the center,
//          the orientation is the same as STc (z is the normal of the surface),
//          but scaled to be in meters.
// ac     : Avatar center Coordinates. This is really the same as stc, but
//          moved up an offset. Rotations passed to llLinkSitTarget are
//          relative to this coordinate system.

//----------------------------------------------------------------------------
//#line "/opt/secondlife/scripts/scripts/tests/llLinkSitTarget_Test.m4s:24"
// The sit position of the sit target.
// This value is *different* in SL and opensim, and even changed since opensim version 0.7.6.
// The currently measured values for the z component are:
// SL: 0.05, opensim 0.7.6: 0.0, opensim 0.8.x: -agentSize.z / 37.9075

// Use this in SL:
//vector sittarget_pos_ac = <0.0, 0.0, 0.05>;
//integer opensimver = 0; // SL

// Use this on opensim:
vector sittarget_pos_ac = <0.0, 0.0, 0.0>;
//integer opensimver = 7; // Opensim 0.7.6
integer opensimver = 8; // Opensim >= 0.8

//----------------------------------------------------------------------------
//#line "/opt/secondlife/scripts/scripts/tests/llLinkSitTarget_Test.m4s:40"
// face_link() returns a list of vectors in link coordinates.
//
// Returns: [bottomLeft_link, topRight_link].
//
// bottomLeft_link is the bottom left corner of the given face,
// that is - the corner where (S, T) = (0, 0).
// topRight_link is the top right corner of the given face,
// that is - the corner where (S, T) = (1, 1).
list face_link(integer link, integer face)
{
  list params = llGetLinkPrimitiveParams(link, [PRIM_TYPE, PRIM_SIZE]);
  vector bottomLeft_link;
  vector topRight_link;
  integer prim_type = llList2Integer(params, 0);
  if (prim_type == PRIM_TYPE_BOX)
  {
    topRight_link = llList2Vector(params, 7) * 0.5;
    if (face == 0)
    {
      bottomLeft_link = -topRight_link;
      bottomLeft_link.z = topRight_link.z;
    }
    else if (face == 1)
    {
      bottomLeft_link = -topRight_link;
      topRight_link.y = bottomLeft_link.y;
    }
    else if (face == 2)
    {
      bottomLeft_link = -topRight_link;
      bottomLeft_link.x = topRight_link.x;
    }
    else if (face == 3)
    {
      topRight_link.x = -topRight_link.x;
      bottomLeft_link = -topRight_link;
      bottomLeft_link.y = topRight_link.y;
    }
    else if (face == 4)
    {
      topRight_link.y = -topRight_link.y;
      bottomLeft_link = -topRight_link;
      topRight_link.x = bottomLeft_link.x;
    }
    else if (face == 5)
    {
      topRight_link.y = -topRight_link.y;
      topRight_link.z = -topRight_link.z;
      bottomLeft_link = -topRight_link;
      bottomLeft_link.z = topRight_link.z;
    }
    else
    {
      llWhisper(0, "This script only supports sitting on boxes without cuts or holes.");
    }
  }
  else
  {
    llWhisper(0, "This script does not work for non-box prims.");
  }
  return [bottomLeft_link, topRight_link];
}

//----------------------------------------------------------------------------
//#line "/opt/secondlife/scripts/scripts/tests/llLinkSitTarget_Test.m4s:104"
default
{
  touch_start(integer num)
  {
    integer n;
    // Run over everyone that might have touch a prim (it could be more than one!)
    for (n = 0; n < num; ++n)
    {
      // Get the link number of the touch prim.
      integer link = llDetectedLinkNumber(n);

      // Some conversion variables we need, copied directly from the table above.
      vector   link_pos_rc = llList2Vector(llGetLinkPrimitiveParams(link, [PRIM_POSITION]), 0);
      rotation link_rot_rc = llList2Rot(llGetLinkPrimitiveParams(link, [PRIM_ROTATION]), 0);

      // Get the face of the prim that was touched.
      integer face = llDetectedTouchFace(n);
      if (face == TOUCH_INVALID_FACE)
      {
        llInstantMessage(llDetectedKey(n), "Sorry, your viewer does not support touched faces.");
      }
      else
      {
        // llDetectedTouchPos returns the position where an object was touched in Region Coordinates.
        vector touch_pos_rc = llDetectedTouchPos(n);
        // Convert that into link coordinates.
        vector touch_pos_link = (touch_pos_rc - link_pos_rc) / link_rot_rc;

        // llDetectedTouchBinormal returns a normalized, directional vector in Region Coordinates.
        // So, it's not a position; it is a direction (dir), really representing a rotation.
        vector binormal_dir_rc = llDetectedTouchBinormal(n);

        // In general a direction vector can not be converted to a rotation (when the
        // starting reference vector is (almost) the same as the direction vector, which
        // might be the case here since we have no idea where the vector is pointing to).
        // Therefore we just leave it as a vector, but we treat it the same as a rotation:
        // leave the origin alone and just correct the rotation when converting to the
        // local coordinate system of the link prim.
        vector binormal_dir_link = binormal_dir_rc / link_rot_rc; // binormal_dir_rc * rc_rot_link

        // llDetectedTouchNormal returns a normalized, directional vector in Region Coordinates.
        // It is perpendicular to the surface (pointing outwards) and makes an angle of 90 degrees
        // with binormal_dir. Here we immediately convert it to _link coordinates.
        vector normal_dir_link = llDetectedTouchNormal(n) / link_rot_rc;

        // Calculate the tangent vector.
        vector tangent_dir_link = binormal_dir_link % normal_dir_link;

        // Get the [bottomLeft_link, topRight_link] positions of the surface 'face', in link coordinates.
        list surface_link = face_link(link, face);
        vector bottomLeft_pos_link = llList2Vector(surface_link, 0);
        vector topRight_pos_link = llList2Vector(surface_link, 1);

        // Calculate the center of the surface.
        vector center_pos_link = (bottomLeft_pos_link + topRight_pos_link) * 0.5;

        // The unit x vector in link coordinates is obviously <1, 0, 0>, center_rot_link is
        // therefore the rotation that is needed to rotate that to tangent_dir_link, which
        // represents <1, 0, 0> in center coordinates (the positive S direction).
        // The positive T direction, or <0, 1, 0> in center coordinates, is the binormal vector.
        // llAxes2Rot(fwd, left, up) returns the rotation needed to rotate <1, 0, 0> to fwd,
        // and <0, 1, 0> to left (and thus <0, 0, 1> to up).
        rotation center_rot_link = llAxes2Rot(tangent_dir_link, binormal_dir_link, normal_dir_link);

        // llDetectedTouchST returns a vector in the plane of the touched surface,
        // with x and y ranging from 0,0 in the "bottom left" corner till 1,1 for
        // opposite corner. The value is thus in ST Coordinates.
        vector touch_pos_STc = llDetectedTouchST(n);

        // Caluclate the size of the surface and rotate it in the plane of the center coordinates.
        vector center_link_scale = (topRight_pos_link - bottomLeft_pos_link) / center_rot_link;

        // Convert touch_pos to center coordinates.
        vector touch_pos_center = <(touch_pos_STc.x - 0.5) * center_link_scale.x, (touch_pos_STc.y - 0.5) * center_link_scale.y, 0.0>;

        // Calculate a unit vector from the touch position to the edge.
        vector edge_dir_center;
        // Find the nearest edge.
        vector edge_pos_center = touch_pos_center;
        if (center_link_scale.x * 0.5 - llFabs(edge_pos_center.x) > center_link_scale.y * 0.5 - llFabs(edge_pos_center.y))
        {
          integer positive = llRound(touch_pos_STc.y);
          edge_pos_center.y = (positive - 0.5) * center_link_scale.y;
          edge_dir_center = <0, positive * 2 - 1, 0>;
        }
        else
        {
          integer positive = llRound(touch_pos_STc.x);
          edge_pos_center.x = (positive - 0.5) * center_link_scale.x;
          edge_dir_center = <positive * 2 - 1, 0, 0>;
        }

        // Convert edge_dir_center to link coordinates (because llLinkSitTarget takes link coordinates).
        vector edge_dir_link = edge_dir_center * center_rot_link;

        // Calculate the rotation of the seated avatar from this direction, using the fact that we still have normal_dir_link.
        // This puts the fwd direction of the avatar in the direction edge_dir (from touch point to nearest edge)
        // and the up direction of the avatar perpendicular to the touched surface.
        rotation ac_rot_link = llAxes2Rot(edge_dir_link, normal_dir_link % edge_dir_link, normal_dir_link);

        // Where do we want to sit?
        // Currently lets just sit where the user clicked, without limits or corrections.
        vector sit_pos_center = touch_pos_center;

        // Convert the sit position to link coordinates (because llLinkSitTarget takes link coordinates).
        vector sit_pos_link = sit_pos_center * center_rot_link + center_pos_link;

        // This vector is really a function of avatar shape and what looks
        // aesthetically pleasing. However, it is only minimal a function
        // avatar height (the only thing we can roughly know from a script)
        // and much more a function of shape values that we can't know.
        // So for now, just leave this as a constant.
        //
        // In ascii art, a sitting avatar looks more or less like this:
        //                ^
        //                |z
        //            ac_ |
        //               \S           'S' is the sittarget_pos_ac, in SL it is 0.05m above(_ac) 'O'.
        //         <--x---O   ()      'O' is the Origin of ac (Avatar center (attach point) Coordinates).
        //                |   /|      '*' is the sittarget (as passed to llSitTarget) and is 0.4m below(_link) S.
        //    sittarget-->*__/_/
        //               _/   ^__ sit_pos_ac, 0.62 below(_ac) 'O', and 0.34 backwards(_ac).
        //                |
        vector sit_pos_ac = <-0.34, 0.0, -0.62>;

        // Correct grid dependent sittarget_pos_ac.
        if (opensimver == 7)
        {
          sittarget_pos_ac.z = 0.0;
        }
        else if (opensimver == 8)
        {
          vector agentSize = llGetAgentSize(llDetectedKey(n));
          sittarget_pos_ac.z = agentSize.z / -37.9075; // Constant determined on "OpenSim 0.8.0 Dev", April 2014.
          llSay(0, "sittarget_pos_ac = " + (string)sittarget_pos_ac);
        }

        // Find the translation of the ac coordinate system.
        vector ac_pos_link = sit_pos_link - sit_pos_ac * ac_rot_link;

        // Convert the sittarget from ac coordinates to link coordinates.
        vector sittarget_pos_link = sittarget_pos_ac * ac_rot_link + ac_pos_link;

        // Grid dependent vertical offset.
        //
        // This is a historical offset introduced by Linden Lab from the beginning.
        // It is actually a mistake because when the avatar is rotated this causes
        // the avatar to move sideways! Apparently also Philips Linden found rotations
        // hard to understand.
        if (opensimver == 7)
        {
          // opensim version 0.7.6
          // Obviously, this should have been the same as on SL.
          // Later it got "corrected", but in a wrong way (sittarget_pos_ac still differs with SL).
          sittarget_pos_link.z -= 0.418;
        }
        else
        {
          // SL and opensim >= 0.8.
          sittarget_pos_link.z -= 0.4;
        }

        // Set the sit target.
        llLinkSitTarget(link, sittarget_pos_link, ac_rot_link);
      }
    }
  }
}