I'm working with some friends on a browser based game where people can move on a 2D map. It's been almost 7 years and still people play this game so we are thinking of a way to give them something new. Since then the game map was a limited plane and people could move from (0, 0) to (MAX_X, MAX_Y) in quantized X and Y increments (just imagine it as a big chessboard).
We believe it's time to give it another dimension so, just a couple of weeks ago, we began to wonder how the game could look with other mappings:

• Unlimited plane with continous movement: this could be a step forward but still i'm not convinced.
• Toroidal World (continous or quantized movement): sincerely I worked with torus before but this time I want something more...
• Spherical world with continous movement: this would be great!

What we want Users browsers are given a list of coordinates like (latitude, longitude) for each object on the spherical surface map; browsers must then show this in user's screen rendering them inside a web element (canvas maybe? this is not a problem). When people click on the plane we convert the (mouseX, mouseY) to (lat, lng) and send it to the server which has to compute a route between current user's position to the clicked point.

What we have We began writing a Java library with many useful maths to work with Rotation Matrices, Quaternions, Euler Angles, Translations, etc. We put it all together and created a program that generates sphere points, renders them and show them to the user inside a JPanel. We managed to catch clicks and translate them to spherical coords and to provide some other useful features like view rotation, scale, translation etc. What we have now is like a little (very little indeed) engine that simulates client and server interaction. Client side shows points on the screen and catches other interactions, server side renders the view and does other calculus like interpolating the route between current position and clicked point.

Where is the problem? Obviously we want to have the shortest path to interpolate between the two route points. We use quaternions to interpolate between two points on the surface of the sphere and this seemed to work fine until i noticed that we weren't getting the shortest path on the sphere surface:

We though the problem was that the route is calculated as the sum of two rotations about X and Y axis. So we changed the way we calculate the destination quaternion: We get the third angle (the first is latitude, the second is longitude, the third is the rotation about the vector which points toward our current position) which we called orientation. Now that we have the "orientation" angle we rotate Z axis and then use the result vector as the rotation axis for the destination quaternion (you can see the rotation axis in grey):

What we got is the correct route (you can see it lays on a great circle), but we get to this ONLY if the starting route point is at latitude, longitude (0, 0) which means the starting vector is (sphereRadius, 0, 0). With the previous version (image 1) we don't get a good result even when startin point is 0, 0, so i think we're moving towards a solution, but the procedure we follow to get this route is a little "strange" maybe?

In the following image you get a view of the problem we get when starting point is not (0, 0), as you can see starting point is not the (sphereRadius, 0, 0) vector, and as you can see the destination point (which is correctly drawn!) is not on the route.

The magenta point (the one which lays on the route) is the route's ending point rotated about the center of the sphere of (-startLatitude, 0, -startLongitude). This means that if i calculate a rotation matrix and apply it to every point on the route maybe i'll get the real route, but I start to think that there's a better way to do this.

Maybe I should try to get the plane through the center of the sphere and the route points, intersect it with the sphere and get the geodesic? But how?

Sorry for being way too verbose and maybe for incorrect English but this thing is blowing my mind!

EDIT: This code version is related to the first image:

public void setRouteStart(double lat, double lng) {
    EulerAngles tmp = new EulerAngles (
            Math.toRadians(lat), 0, -Math.toRadians(lng));
    //set route start Quaternion
    qtStart.setInertialToObject(tmp);
    //do other stuff like drawing start point...
}

public void impostaDestinazione(double lat, double lng) {
    EulerAngles tmp = new AngoliEulero(
       Math.toRadians(lat), 0, -Math.toRadians(lng));
    qtEnd.setInertialToObject(tmp);
    //do other stuff like drawing dest point...
}

public V3D interpolate(double totalTime, double t) {
    double _t = t/totalTime;
    Quaternion q = Quaternion.Slerp(qtStart, qtEnd, _t);
    RotationMatrix.inertialQuatToIObject(q);
    V3D p = matInt.inertialToObject(V3D.Xaxis.scale(sphereRadius));
    //other stuff, like drawing point ...
    return p;
}

//mostly taken from a book!
public static Quaternion Slerp(Quaternion q0, Quaternion q1, double t) {
    double cosO = q0.dot(q1);
    double q1w = q1.w;
    double q1x = q1.x;
    double q1y = q1.y;
    double q1z = q1.z;
    if (cosO < 0.0f) {
            q1w = -q1w;
            q1x = -q1x;
            q1y = -q1y;
            q1z = -q1z;
            cosO = -cosO;
    }

    double sinO = Math.sqrt(1.0f - cosO*cosO);
      double O = Math.atan2(sinO, cosO);
      double oneOverSinO = 1.0f / senoOmega;
      k0 = Math.sin((1.0f - t) * O) * oneOverSinO;
      k1 = Math.sin(t * O) * oneOverSinO;

    // Interpolate
    return new Quaternion(
            k0*q0.w + k1*q1w,
            k0*q0.x + k1*q1x,
            k0*q0.y + k1*q1y,
            k0*q0.z + k1*q1z 
            );
}

A little dump of what i get (again check image 1):

Route info:
Sphere radius and center: 200,000, (0.0, 0.0, 0.0)

Route start:  lat 0,000 °, lng 0,000 °  @v: (200,000,    0,000,   0,000),   |v| = 200,000
Route end:    lat 30,000 °, lng 30,000 °  @v: (150,000,    86,603,   100,000),   |v| = 200,000
Qt dump: (w, x, y, z), rot. angle°, (x, y, z) rot. axis
Qt start: (1,000,    0,000,   -0,000,   0,000);  0,000 °; (1,000, 0,000, 0,000)
Qt end: (0,933,    0,067,   -0,250,   0,250);  42,181 °; (0,186, -0,695, 0,695)

Route start:  lat 30,000 °, lng 10,000 °  @v: (170,574,    30,077,   100,000),   |v| = 200,000
Route end:    lat 80,000 °, lng -50,000 °  @v: (22,324,    -26,604,   196,962),   |v| = 200,000
Qt dump: (w, x, y, z), rot. angle°, (x, y, z) rot. axis
Qt start: (0,962,    0,023,   -0,258,   0,084);  31,586 °; (0,083, -0,947, 0,309)
Qt end: (0,694,    -0,272,   -0,583,   -0,324);  92,062 °; (-0,377, -0,809, -0,450)