How do you model the position and orientation of an airplane?
Position is easy, just represent it with a point in 3D space. But how do you specify its orientation — which direction it’s pointing?
At first glance, it seems a vector will do. After all, a vector points in some direction, right? If the plane is pointing east, represent its orientation by a unit vector pointing east.
Unfortunately, we quickly run into trouble when we try to roll. If we’re facing east, and we roll 90 degrees, we’re still facing east. Clearly we’re missing something.
Euler Angles
When real pilots talk about their orientation, they talk about roll, yaw, pitch. Pitch is going up or down, yaw is going left or right, roll is, well, roll.
Any change in orientation can be described by some combination of roll, yaw, pitch. This is the basis for Euler Angles. We use three angles to represent the airplane’s orientation.
This is all fine and dandy if we want to represent the orientation of a static object in space. But when we try to adjust our orientation, we start to run into problems.
You’re thinking, this should be simple! When we turn left or right, we just increment the yaw variable, right? Yes, it seems to work, at least initially. You can turn left and right, up and down, and roll around.
Implement it in Unity and play around a bit, however, and you begin to notice that things don’t quite behave the way you expect.
In this animation, I’m holding down the right button:
The plane does rotate to the right, but it’s not rotating relative to itself. Instead it’s rotating around some invisible y-axis. If it was rotating relative to itself, the green arrow shouldn’t be moving.
The problem becomes more and more severe when the pitch of the plane becomes higher and higher. The worst case is when the airplane is pointing straight up: then roll and yaw become the same thing! This is called gimbal lock: we have lost a degree of freedom and we can only rotate in 2 dimensions! Definitely not something desirable if we’re controlling a plane or spaceship.
It turns out that no matter what we do, we will suffer from some form of gimbal lock. As long as we use Euler Angles, there is one direction where if we turn too far, everything starts to screw up.
Practical Introduction to Quaternions
All is not lost, however. There is a way to represent orientation that represents all axes equally and does not suffer from gimbal lock. This mythical structure is called the quaternion. Unlike Euler Angles which describe your orientation relative to a fixed set of axes, quaternions do not rely on any fixed axis.
The drawback is that quaternions are unintuitive to understand for humans. There is no way to “look” at a quaternion and be able to visualize what rotation it represents. Fortunately for us, it’s not that difficult to make use of quaternions, even if we can’t visualize quaternions.
There is a lot of theory behind how quaternions work, but in this article, I will gloss over the theory and give a quick primer to quaternions, just the most common facts you need to use them. At the same time, I will implement the operations I describe in C#, so I can integrate them with Unity. If you don’t know C#, you can freely ignore the code.
Definition
A quaternion is an ordered pair of 4 real numbers (w,x,y,z). We write this as
The letters i,j,k are not variables. Rather, they are independent axes. If you like, you can think of the quaternions as a 4 dimensional vector space.
The defining property of quaternions is:
Play around with it a bit and you can derive 6 more identites:
If you’ve worked with complex numbers, this should seem familiar. Instead of 2 parts of a complex number (the real and imaginary parts), we have 4 parts for a quaternion.
The similarity doesn’t end here. Multiplying complex numbers represents a rotation in 2 dimensions. Similarly, multiplying by a quaternion represents a rotation in 3D.
One curious thing to note: we have 
Here’s a framework for a quaternion in C#:
public class Quat{
// Represents w + xi + yj + zk
public float w, x, y, z;
public Quat(float w, float x, float y, float z){
this.w = w;
this.x = x;
this.y = y;
this.z = z;
}
}
Normalizing Quaternions
The norm of a quaternion is
When we use quaternions to represent rotations, we typically want unit quaternions: quaternions with norm 1. This is straightforward: to normalize a quaternion, divide each component by the norm.
In C#:
public float Norm(){
return Mathf.Sqrt (w * w + x * x + y * y + z * z);
}
public Quat Normalize(){
float m = Norm ();
return new Quat (w / m, x / m, y / m, z / m);
}
Multiplying Quaternions
Multiplying is simple, just a little tedious. If we have two quaternions:
Then their product is this ugly mess:
In C#:
// Returns a*b
public static Quat Multiply(Quat a, Quat b){
float w = a.w * b.w - a.x * b.x - a.y * b.y - a.z * b.z;
float x = a.w * b.x + a.x * b.w + a.y * b.z - a.z * b.y;
float y = a.w * b.y + a.y * b.w - a.x * b.z + a.z * b.x;
float z = a.w * b.z + a.z * b.w + a.x * b.y - a.y * b.x;
return new Quat (w,x,y,z).Normalize();
}
Since multiplication is not commutative, I made this function static to avoid confusing left and right multiplication. Also, I normalize the product so that floating point errors don’t accumulate.
Constructing Rotation Quaternions
Every rotation operation can be written as a rotation of some angle, 
The following formula gives a quaternion that represents this rotation:
For our purposes,
That’s it: given any quaternion, multiplying on the left by our quaternion rotates it slightly around the x axis.
In C#, our code might look like this:
Quat qx = new Quat (Mathf.Cos (0.01 / 2), 0, 0, Mathf.Sin (0.01 / 2)); Quat qy = new Quat (Mathf.Cos (0.01 / 2), 0, Mathf.Sin (0.01 / 2), 0); Quat qz = new Quat (Mathf.Cos (0.01 / 2), Mathf.Sin (0.01 / 2), 0, 0);
Putting it together
That’s all we need to do interesting things with quaternions. Let’s combine everything we have. Here’s our quaternion class thus far:
public class Quat{
// Represents w + xi + yj + zk
public float w, x, y, z;
public Quat(float w, float x, float y, float z){
this.w = w;
this.x = x;
this.y = y;
this.z = z;
}
public float Norm(){
return Mathf.Sqrt (w * w + x * x + y * y + z * z);
}
public Quat Normalize(){
float m = Norm ();
return new Quat (w / m, x / m, y / m, z / m);
}
// Returns a*b
public static Quat Multiply(Quat a, Quat b){
float w = a.w * b.w - a.x * b.x - a.y * b.y - a.z * b.z;
float x = a.w * b.x + a.x * b.w + a.y * b.z - a.z * b.y;
float y = a.w * b.y + a.y * b.w - a.x * b.z + a.z * b.x;
float z = a.w * b.z + a.z * b.w + a.x * b.y - a.y * b.x;
return new Quat (w,x,y,z).Normalize();
}
public Quaternion ToUnityQuaternion(){
return new Quaternion (w, x, y, z);
}
}
Now we just need to read the input, perform our calculations, and output the rotation quaternion to Unity:
public class Airplane : MonoBehaviour {
public GameObject airplane;
public Quat quat = new Quat (0, 0, 0, -1);
public float speed = 0.01f;
void FixedUpdate(){
float inputX = Input.GetAxis("UpDown");
float inputY = Input.GetAxis("LeftRight");
float inputZ = Input.GetAxis("Roll");
Quat qx = new Quat (Mathf.Cos (speed * inputX / 2), 0, 0, Mathf.Sin (speed * inputX / 2));
Quat qy = new Quat (Mathf.Cos (speed * inputY / 2), 0, Mathf.Sin (speed * inputY / 2), 0);
Quat qz = new Quat (Mathf.Cos (speed * inputZ / 2), Mathf.Sin (speed * inputZ / 2), 0, 0);
quat = Quat.Multiply (qx, quat);
quat = Quat.Multiply (qy, quat);
quat = Quat.Multiply (qz, quat);
airplane.transform.rotation = quat.ToUnityQuaternion ();
}
}
In Unity, the input is not given to us as a single true/false value, but a float between -1 and 1. So holding right increases the LeftRight input gradually until it reaches 1, avoiding a sudden jump in movement.
What’s ToUnityQuaternion? Well, it turns out that Unity already has a Quaternion class that does everything here and much more, so all this could have literally been implemented in one line if we wanted.
Anyways, let’s see the result.
As you can see, holding right turns the plane relative to itself now, and the green arrow stays still. Hooray!
, it follows that when both ships travel to such a point P, they will meet at the same time.
, so it’s easy to (falsely) conclude that the angle between two faces must also be 
, 
, and 
; let 
and 
be the angles between the vectors (by symmetry we’re assuming that the two alpha’s are equal):
and 
be the perpendiculars to 
with 
and 
.
. Hence the cosine of the angle is:
and 
.
but 
:
, 
, and 
.
(since it’s part of a regular pentagon) and 
. Then 
and 
.
and 
. Let 
, 
, 
, 
, and 
. Write an algorithm to find triangles satisfying these conditions where a, b, c are all integers.
is completely unnecessary, but if we take it out, there doesn’t seem to be any way to relate the angle ratios to integer side lengths.
is similar to 
. Being an exterior angle, 
, and also 
. Both of the triangles have an angle of measure 
and another angle of measure 
, thus they are similar.
can be written as 
. Then we have the ratio
and 
is an isosceles triangle.
. Also with this construction, 
, and is also isosceles, hence 
.
is a perfect square, that is:
:
in the equation 
, we have the equation 
.
in the equation 
and simplifying, we get the original hyperbola, 
. Stretching the hyperbola horizontally by a factor of 2, this point gets transformed to 
. Even though each point has a different image point, however, the graph as a whole remains identical.
instead of 2, where 
. We call this combinations of transformations a hyperbolic rotation. Intuitively, the points on the hyperbola are ‘rotated’ downwards.
be the line passing through P tangent to the hyperbola. Let X be the intersection of 
. The theorem obviously holds here, because of symmetry. But by applying a hyperbolic rotation, we can transform the point 
into any other point on the hyperbola. Since ratios between line segments don’t change with a hyperbolic rotation, P is always the midpoint, and we are done.
evenly spaced points in the interval ![[0,1]](https://s0.wp.com/latex.php?latex=%5B0%2C1%5D&bg=ffffff&fg=393939&s=0&c=20201002)
for both the horizontal and vertical axis and connecting them.
on the x-axis, then obviously that point is connected to point 
on the y-axis. The slope of this line is 
and the equation of this line is 
. More conveniently, we can this expression as 
, intuitively the line formed when we choose 
as our x-axis point.
be the bounding curve for the blue lines. To find 
:
, so solving for 
. Then we have
, and the latitude (north-south location) as 
is the center of the earth, and 
a point on the equator.
on the surface of the globe will be projected in a straight line to point 
on the cylinder such that points 
; since 
, the length of 
is 
.
, the slope of 
becomes increasingly steep, and the 
south to the parallel. With a compass, it would be easy for you to maintain this constant heading for the entire duration of the trip.
is parallel to 
and 
is perpendicular to 
, the circumference of the Earth at the equator is 
. But the circumference at latitude 
.
in the diagram, it follows that angles 
and 
are equal. Next, 
, therefore 
. The result follows.
; in the projection any line at latitude 
or 
.
for some small number 
. Similarly, the space this piece of land occupies on the map is not just 
for some small number 
. Now as 
and integrating both sides gives:
be the denominator, 
, it turns out that the numerator is the derivative of 
. Thus we can substitute 
, 
. Substituting and solving for 
. Also, since 
for 
, we can remove the absolute value brackets. Therefore,
, the map is practically unusable.
, with all angles smaller than 
. Point 
is inside the triangle such that 
is at its minimum (
, 
, 
, 
, 
are all integers.
.
being the point such that 
is minimized, 
is equilateral, then construct the circumcircle of 
and the circumcircle of 
.
is cyclic. Since 
so we can divide it out to get
to both sides:
, so 
, or by substitution,
, and the sum 
, and 
:
, 
, and 
are all cyclic. As 
, 
, and similarly 
and 
.
):
is a perfect square.
, 
, 
are all pairs.
and 
in the list, search through the list to check if there exists some 
is a pair and 
is a pair. We index the list to drastically improve searching time.
triple is found and 
, then mark 
as found. This is easily implemented as a bit array of size 120000 with all bits initially 0 and flipped to 1 when a sum is encountered.
operations, which is in the billions. An implementation of this algorithm would take about 3 minutes.
can be generated by 
, 
, 
where 
and 
equation?
, and 
and dividing by 
, we get
. This is equivalent to finding all rational points on the ellipse 
:
to be the first point of the line. Then the equation of the line is 
, where 
is valid only if it is positive and greater than 1.
where 
.
, 
, and 
.
. We have 
which is divisible by 3. Next 
which is divisible by 3. Finally, 
, again divisible by 3! So if 
, then the resulting triple is not primitive!
, with 
.
, we have 
and also 
. Thus we can have 
where 
, and 
. Combining the two and substituting we get
operations. So we can loop up to 
, which cuts down the time down to about 2 seconds:
Lucky's Notes 