Hello,

Rather than the normal way of having latitude/longitude, with a north/south pole radiating lines for longitude, and then latitude lines being parallel circles of various radii, I instead have *two* poles, so a north/south pole and an east/west pole both radiate lines. So both sets of lines being like segments of an orange, coming together at the poles.

As in this rough diagram, the two poles shown as red and blue:



blue pole marked x is also the x axis (+ve to the right)
red pole marked y is also the y axis (+ve going upwards)
+z axis is coming out of the screen towards you

Now, I am managing to get a 3D point from a given lat/long, using these equations:

real_x = worldRadius * cos(latitude) * sin(longitude)

real_y = worldRadius * sin(latitude) * cos(longitude)

real_z = worldRadius * cos(latitude) * cos(longitude)

However, I also need to be able to do the reverse, from a given 3D point on the sphere, calculate the lat/long. This is giving me great trouble.

I thought this might work, and it does in *some* situations, however in other cases one or the other (or even both) lat and long are out by 180 degrees (and when put into the previous equations come out to a point directly opposite to what they should be):

lat = atan2(real_y, real_z)

lng = atan2(real_x, real_z)

Any thoughts? Math is really not my strong point.

You are talking about Spherical Polar Coordinates

http://en.wikipedia.org/wiki/Spherical_coordinate_system

Should have all the info you need...

What's an example of a point where using the two atan2s doesn't do what you want? I think that ought to work...

What's an example of a point where using the two atan2s doesn't do what you want? I think that ought to work...

Here's some test code I'm using:

		worldRadius = 10.0
		
		for i in range(0,10):
			# choose a random lat/long point on sphere
			latitude = random.uniform(-2*pi, 2*pi) # up/down, y
			longitude = random.uniform(-2*pi, 2*pi) # left/right, x
			
			# Find the real world 3D position on the sphere surface
			real_x = worldRadius * cos(latitude) * sin(longitude)
			real_y = worldRadius * sin(latitude) * cos(longitude)
			real_z = worldRadius * cos(latitude) * cos(longitude)
			
			# Now attempt to do that backwards and get the lat/long again from this real world position.
			latitudeAgain = atan2(real_y, real_z)
			longitudeAgain = atan2(real_x, real_z)
			
			# Put all lat/longs in the range 0 - 2*PI for easy comparison
			latitude += 2*pi
			longitude += 2*pi
			latitudeAgain += 2*pi
			longitudeAgain += 2*pi
			latitude = fmod(latitude, 2*pi)
			longitude = fmod(longitude, 2*pi)
			latitudeAgain = fmod(latitudeAgain, 2*pi)
			longitudeAgain = fmod(longitudeAgain, 2*pi)
			
			# Also try to get the position once again from the newly got lat/long
			other_x = worldRadius * cos(latitudeAgain) * sin(longitudeAgain)
			other_y = worldRadius * sin(latitudeAgain) * cos(longitudeAgain)
			other_z = worldRadius * cos(latitudeAgain) * cos(longitudeAgain)
			
			print "Lats: ",latitude," - ",latitudeAgain," Diff: ",abs(latitude-latitudeAgain)
			print "Lons: ",longitude," - ",longitudeAgain," Diff: ",abs(longitude-longitudeAgain)
			print "Pos1: ",real_x,real_y,real_z
			print "Pos2: ",other_x,other_y,other_z

And the test output:

Lats:  1.57903473033  -  1.57903473033  Diff:  0.0
Lons:  5.10226596478  -  1.96067331119  Diff:  3.14159265359
Pos1:  0.0762007311343 3.80061735995 -0.0313117278778
Pos2:  -0.0762007311343 -3.80061735995 0.0313117278778

Lats:  5.93480811564  -  2.79321546205  Diff:  3.14159265359
Lons:  3.39572137905  -  3.39572137905  Diff:  0.0
Pos1:  -2.36299962085 3.30408958826 -9.0973998891
Pos2:  2.36299962085 -3.30408958826 9.0973998891

Lats:  2.36983896676  -  5.51143162035  Diff:  3.14159265359
Lons:  3.43495760648  -  0.293364952892  Diff:  3.14159265359
Pos1:  2.07248505761 -6.67597897677 6.86069121889
Pos2:  2.07248505761 -6.67597897677 6.86069121889

Lats:  2.71182886148  -  2.71182886148  Diff:  4.4408920985e-16
Lons:  5.72284029962  -  2.58124764603  Diff:  3.14159265359
Pos1:  4.83148051356 3.52937619571 -7.70042639429
Pos2:  -4.83148051356 -3.52937619571 7.70042639429

Lats:  4.41217625779  -  1.2705836042  Diff:  3.14159265359
Lons:  3.85396423175  -  0.712371578164  Diff:  3.14159265359
Pos1:  1.9329383337 7.22964312956 2.23807587188
Pos2:  1.9329383337 7.22964312956 2.23807587188

Lats:  0.135838374462  -  3.27743102805  Diff:  3.14159265359
Lons:  1.64115450181  -  1.64115450181  Diff:  0.0
Pos1:  9.88336821577 -0.0952011605369 -0.696525460481
Pos2:  -9.88336821577 0.0952011605369 0.696525460481

Lats:  5.03750570722  -  5.03750570722  Diff:  0.0
Lons:  1.46803830938  -  1.46803830938  Diff:  0.0
Pos1:  3.17734475072 -0.972036026841 0.327651707624
Pos2:  3.17734475072 -0.972036026841 0.327651707624

Lats:  0.350045505601  -  0.350045505601  Diff:  0.0
Lons:  5.436873812  -  5.436873812  Diff:  0.0
Pos1:  -7.0342906832 2.27283766805 6.22558689914
Pos2:  -7.0342906832 2.27283766805 6.22558689914

Lats:  0.0171541809941  -  0.0171541809941  Diff:  0.0
Lons:  1.13445981395  -  1.13445981395  Diff:  0.0
Pos1:  9.06172667879 0.0724937989497 4.22559887977
Pos2:  9.06172667879 0.0724937989497 4.22559887977

Lats:  1.69526847861  -  4.8368611322  Diff:  3.14159265359
Lons:  2.92736309538  -  6.06895574897  Diff:  3.14159265359
Pos1:  -0.263938369219 -9.69580670745 1.21312953854
Pos2:  -0.263938369219 -9.69580670745 1.21312953854

As you can see it works in a lot of cases, shown by the "diff" between the original and the newly gotton lat/lon being 0.0 (and Pos1 and Pos2 being equal).

Funnily, when both the lat and long are out by 180 degrees the correct position is also returned.

In the cases of failure, only one of the returned lat or long is out by 180 degrees, and Pos2 is directly opposite (negative) to Pos1.

So it's really just the sign of the represented point that wrong in some cases. The actual absolute values are correct.

Your parametric formulas to obtain x, y and z from latitude and longitude don't describe a sphere (the squares of the numbers obtained don't sum to a constant).

This is some test code showing how I think it should be done:

#include <iostream>
#include <cstdlib>
#include <cmath>

double const tau = 8*std::atan(1.0);

double U(double a, double b) {
  return a + (b-a) * std::rand() / RAND_MAX;
}

int main() {
  for (int k=0; k<10; ++k) {
    double z = U(-1,1);
    double alpha = U(0,tau);
    double x = std::cos(alpha)*std::sqrt(1.0-z*z);
    double y = std::sin(alpha)*std::sqrt(1.0-z*z);

    double a1 = std::atan2(y, z);
    double a2 = std::atan2(x, z);

    std::cout << '(' << x << ',' << y << ',' << z << ") "
              << a1 << ' ' << a2 << '\n';

    double new_z = (std::cos(a1) > 0 ? 1 : -1)
      * std::sqrt(1.0/(std::pow(std::tan(a1),2)+
                       std::pow(std::tan(a2),2)+1));
    double new_x = std::tan(a2) * new_z;
    double new_y = std::tan(a1) * new_z;

    std::cout << '(' << new_x << ',' << new_y << ',' << new_z << ")\n";
  }
}

THANKS! That works perfectly. My original stuff was perhaps inspired by the math school of trial and error.