代码几乎可以完成的工作：

• 计算差异向量
• 将矢量转换为spherical coordinates
• 将w_nx和wn_y添加到倾角和方位角（有关术语，请参见链接）
• 将修改的球坐标转换回笛卡尔坐标

有两个问题：

• 转换不正确，您所做的计算是针对两个倾斜向量（一个沿着x轴，另一个沿y轴）
• 即使计算是正确的，球坐标中的变换也不同于绕两个轴旋转

因此，在这种情况下，使用矩阵和向量数学将有所帮助：

b = p - m
b = RotationMatrixAroundX(wn_x) * b
b = RotationMatrixAroundY(wn_y) * b
n = m + b

基本旋转矩阵。

尝试使用矢量数学。决定旋转的顺序，首先沿着x，然后沿着y。
如果沿着z-axis旋转，则[z' = z]

x' = x*cos a - y*sin a;
y' = x*sin a + y*cos a;

对y-axis重复相同的步骤：[y'' = y']

x'' = x'*cos b - z' * sin b;
z'' = x'*sin b + z' * cos b;

再次沿着x-axis旋转：[x''' = x'']

y''' = y'' * cos c - z'' * sin c
z''' = y'' * sin c + z'' * cos c

最后是围绕某个特定“点”旋转的问题：
首先，从坐标中减去该点，然后应用旋转，最后将该点添加回结果。
这个问题，据我所见，是一个密切相关的“万向节锁定”。Angular w_ny不能相对于固定的xyz坐标系测量，但相对于坐标系，是通过应用Angular w_nx旋转。
作为kakTuZ observed，您的代码将点坐标转换为球坐标。这本身没有什么问题--有了经度和纬度，我们可以到达地球上的所有地方。如果我们不关心地球赤道平面相对于其绕太阳的轨道的倾斜，我也可以接受。
不沿着第一个w_ny旋转下一个参考轴的结果是，在赤道上彼此相距1 km的两个点在两极和90度纬度上彼此靠近，它们接触，尽管表面上的目的是使它们在旋转的任何地方都相距1 km。

如果你想转换坐标系而不是仅仅点，你需要3个Angular 。但是你是对的-对于转换点2个Angular 就足够了。详细信息请询问Wikipedia...
但是当你使用opengl时，你真的应该使用像glRotatef这样的opengl函数。这些函数将在GPU上计算--而不是像你的函数那样在CPU上计算。doc是here。

就像很多人说的那样，你应该使用glRotatef来旋转它以进行渲染。对于冲突处理，你可以通过将它的位置向量乘以渲染点堆栈顶部的OpenGL ModelView矩阵来获得它的世界空间位置。使用glGetFloatv获得该矩阵，然后将它乘以你自己的向量矩阵乘法函数，或者使用您可以在网上轻松获得的众多工具中的一种。
但是，这将是一个痛苦！相反，考虑使用GL反馈缓冲区。这个缓冲区将只存储图元将被绘制的点，而不是实际绘制图元，然后你可以从那里访问它们。
This是一个很好的起点。