球面坐标系比较常见的场景就是地理上的经纬度了，以南北极为地轴，地轴通过球心的垂直面为赤道面。
如果把地轴转动，北极从北极点转到任意一点，下面计算球面上的点在这个新的球面坐标系的坐标值（经纬度）。
计算球面上任意两点相对于球心的夹角
原理：设a,b是两个不为0的向量，它们的夹角为θ，则有向量公式：cosθ=a*b/|a||b|
假设北极从北极点[0°,90°]移动到点1[longitude1,latitude1]，对于任意的点2[longitude2,latitude2]，其新的纬度与点1到球心与点2到球心的夹角相关。
将点1和点2换算成直角坐标系坐标：点1(x1,y1,z1)、点2(x2,y2,z2)
球心到点1的向量就是a=(x1,y1,z1)；球心到点2的向量就是b=(x2,y2,z2)；
则有：
1
| cosθ=(x1,y1,z1)*(x2,y2,z2)/(|(x1,y1,z1)|*|(x2,y2,z2)|)
|
其中：
1
2
3
| (x1,y1,z1)*(x2,y2,z2)=(x1x2+y1y2+z1z2)
|a|=√(x1^2+y1^2+z1^2)
|b|=√(x2^2+y2^2+z2^2)
|
用坐标表示的向量公式：
1
| cosθ=(x1x2+y1y2+z1z2)/(√(x1^2+y1^2+z1^2)*√(x2^2+y2^2+z2^2))
|
通过这个公式可以算出点1和点2的夹角(取值范围是：[0,π]；夹角为锐角时，cosθ>0；夹角为钝角时,cosθ<0)
计算球面上任意点2相对于点1的纬度
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
| private static final double HD = (Math.PI / 180); //Math类三角函数计算以弧度为单位
public static double getLatitudeAngle(double longitude1, double latitude1, double longitude2, double latitude2) {
//计算球心到点1的向量
double z1 = Math.sin(latitude1 * HD);
double xy1 = Math.cos(latitude1 * HD);
double x1 = xy1 * Math.cos(longitude1 * HD);
double y1 = xy1 * Math.sin(longitude1 * HD);
//计算球心到点2的向量
double z2 = Math.sin(latitude2 * HD);
double xy2 = Math.cos(latitude2 * HD);
double x2 = xy2 * Math.cos(longitude2 * HD);
double y2 = xy2 * Math.sin(longitude2 * HD);
//根据向量公式计算出cosθ的值
double cos = (x1 * x2 + y1 * y2 + z1 * z2)
/ (Math.sqrt(Math.pow(x1, 2) + Math.pow(y1, 2) + Math.pow(z1, 2))
* Math.sqrt(Math.pow(x2, 2) + Math.pow(y2, 2) + Math.pow(z2, 2)));
//浮点数异常范围处理：由于浮点数运算的精度问题，计算出cosθ的值可能在[-1,1]之外，这时Math.acos方法会得出NaN
if (cos > 1) {
cos = 1;
} else if (cos < -1) {
cos = -1;
}
//夹角换算成纬度值
return (90 - (Math.acos(cos) / HD));
}
|
计算三维直角坐标系上两个平面的夹角
1.已知不在同一直线上的三个点的坐标，求三个点所在平面的平面方程：
假设三个点的坐标分别为P1(x1,y1,z1)，P2(x2,y2,z2)，P3(x3,y3,z3)
假设通过P1，P2，P3三点的平面方程为A(x - x1) + B(y - y1) + C(z - z1) = 0
化简为一般式：Ax + By + Cz + D = 0
1
2
3
| Ax1 + By1 + Cz1 + D = 0
Ax2 + By2 + Cz2 + D = 0
Ax3 + By3 + Cz3 + D = 0
|
求A、B、C、D的值：
1
2
3
4
| A = (y2 - y1)*(z3 - z1) - (z2 -z1)*(y3 - y1);
B = (x3 - x1)*(z2 - z1) - (x2 - x1)*(z3 - z1);
C = (x2 - x1)*(y3 - y1) - (x3 - x1)*(y2 - y1);
D = -(A * x1 + B * y1 + C * z1)
|
2.求平面的法向量：
平面”Ax + By + Cz + D = 0”的法线向量为：(A,B,C)
3.两个平面法向量的夹角和这两个平面的夹角相等，两个向量的夹角求法见前面的公式，即：
1
| cosθ=(x1x2+y1y2+z1z2)/(√(x1^2+y1^2+z1^2)*√(x2^2+y2^2+z2^2))
|
计算球面上任意一点(Point2)对另外一点(Point1)的经度角(新坐标系以南极点为0°经度)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
| private static final double HD = (Math.PI / 180); //Math类三角函数计算以弧度为单位
public static double getLongitudeAngle(double longitude1, double latitude1, double longitude2, double latitude2) {
//计算点1的坐标
double z1 = Math.sin(latitude1 * HD);
double xy1 = Math.cos(latitude1 * HD);
double x1 = xy1 * Math.cos(longitude1 * HD);
double y1 = xy1 * Math.sin(longitude1 * HD);
//计算点2的坐标
double z2 = Math.sin(latitude2 * HD);
double xy2 = Math.cos(latitude2 * HD);
double x2 = xy2 * Math.cos(longitude2 * HD);
double y2 = xy2 * Math.sin(longitude2 * HD);
//球心和北极点的坐标
double x0 = 0, y0 = 0, z0 = 0;
double xn = 0, yn = 0, zn = 1;
//求 球心，点1，北极点的平面法向量F1(A1,B1,C1)
double A1 = (y1 - y0) * (zn - z0) - (z1 - z0) * (yn - y0);
double B1 = (x1 - x0) * (zn - z0) - (z1 - z0) * (xn - x0);
double C1 = (y1 - y0) * (xn - x0) - (x1 - x0) * (yn - y0);
//求 球心，点1，点2的平面方程法向量F2(A2,B2,C2)
double A2 = (y1 - y0) * (z2 - z0) - (z1 - z0) * (y2 - y0);
double B2 = (x1 - x0) * (z2 - z0) - (z1 - z0) * (x2 - x0);
double C2 = (y1 - y0) * (x2 - x0) - (x1 - x0) * (y2 - y0);
//求两个平面法向量的夹角(F1、F2)
double labc = (Math.sqrt(Math.pow(A1, 2) + Math.pow(B1, 2) + Math.pow(C1, 2))
* Math.sqrt(Math.pow(A2, 2) + Math.pow(B2, 2) + Math.pow(C2, 2)));
if (labc == 0) {
//平面的三点在同一条直线，经度值显示0°
return 0; //浮点数计算时0f/0f会得出NaN
}
double cos = (A1 * A2 + B1 * B2 + C1 * C2) / labc;
//浮点数异常范围处理：由于浮点数运算的精度问题，计算出cosθ的值可能在[-1,1]之外，这时Math.acos方法会得出NaN
if (cos > 1) {
cos = 1;
} else if (cos < -1) {
cos = -1;
}
//求出两个平面(法向量)的夹角[0°,180°]
double theta = Math.acos(cos) / HD;
//计算经度的正负号和校准基线(新的坐标系：南极点为0°经度,北极点为180°经度)
int sign;
if (longitude1 > longitude2) {
sign = (360 + longitude2 - longitude1) >= 180 ? 1 : -1;
} else {
sign = (longitude2 - longitude1) >= 180 ? 1 : -1;
}
return ((sign > 0 ? theta : (360 - theta)) + 180) % 360;
}
|
结论
球坐标极轴北极点移到点1(longitude1,latitude1)时，球面上任意一点(longitude2,latitude2)在新的球面坐标系的坐标用上述工具方法表示为：
1
| (Plongitude,Platitude) = (getLongitudeAngle(longitude1,latitude1,longitude2,latitude2), getLatitudeAngle(longitude1,latitude1,longitude2,latitude2))
|
注意事项
1.经度计算时可能出现球心、点1、点2在同一条直线的情况，即求夹角时(cosθ=a*b/|a||b|)可能出现分母为0的情况，需特殊处理；
2.java中的浮点数计算精度问题，可能会出现计算出的cosθ值在[-1,1]之外，需特殊处理；
源码请参考：
https://github.com/guangGG/AndroidStarrySky