这篇文章描述2件事情：

1、已知球面上任意点，求过该点、地心、与北极点的平面方程（即过该点的经线平面方程）；

2、绕过球心的任意轴旋转平面得到新平面的方程

一、已知球面上任意点，求过该点、地心、与北极点的平面方程（即过该点的经线平面方程）

输入经纬度，输出过该点，穿过地心、与北极点的平面方程，输出参数是平面方程的参数。

平面方程基本形式为：A*X+B*Y+C*Z+D=0 

%% 已知球面上任意点，求过该点、地心、与北极点的平面方程（即过该点的经线平面方程）
%% 输入参数：任意点的经纬度(LNG,LAT)
%% 输出参数：平面方程的参数，平面方程的表达式为A*X+B*Y+C*Z+D=0
function [A,B,C,D]=GPS2EQPlane(LNG,LAT)
%% 地理常数
R=6371;%地球半径，单位km

%% 经纬度转地理坐标系坐标
xt=R*cosd(LAT)*cosd(LNG);
yt=R*cosd(LAT)*sind(LNG);
zt=R*sind(LAT);

%% 平面计算
% 计算单位法向量
nx = -yt/(sqrt(xt^2+yt^2));
ny = xt/(sqrt(xt^2+yt^2));
nz = 0;

% 法向量
n = [nx; ny; nz];

%% 求解平面方程的参数形式
A = n(1);
B = n(2);
C = n(3);
D = -(A * xt + B * yt + C * zt);
end

二、绕过球心的任意轴旋转平面得到新平面

%% 绕过球心的任意轴旋转平面得到新平面
% 输入参数
% P:球面上的一点坐标，在直角坐标系下（x,y,z)
% N：原平面的法向量
% theta：绕旋转轴旋转的角度，单位度
% 输出参数：新平面的参数A*X+B*Y+C*Z+D=0
function [newA, newB, newC, newD] = rotatePlane(P, N, theta)
%% 地理常数
R=6371;%地球半径，单位km

%% 计算点 P到球心的向量
OP = P; % 假设球心在原点，P 是点 P 的笛卡尔坐标
% 旋转轴的单位化
A = OP / norm(OP);
    
%% 计算右乘旋转矩阵
C=cosd(theta);
S=sind(theta);
R=[C+A(1)^2*(1-C) A(1)*A(2)*(1-C)-A(3)*S A(1)*A(3)*(1-C)+A(2)*S;
   A(1)*A(2)*(1-C)+A(3)*S C+A(2)^2*(1-C) A(2)*A(3)*(1-C)-A(1)*S;
   A(1)*A(3)*(1-C)-A(2)*S A(2)*A(3)*(1-C)+A(1)*S C+A(3)^2*(1-C)];

%% 旋转原平面的法向量
rotated_N = N*R;

%% 计算新的平面方程的系数
newA = rotated_N(1);
newB = rotated_N(2);
newC = rotated_N(3);
newD = -(newA * P(1) + newB * P(2) + newC * P(3));
end

三、示例

3.1 代码1示例

以海南凤凰机场为例：

三亚凤凰国际机场位于经度:109.414871、纬度:18.303421。

3.1.1 测试代码

clc;clear;close all
%% 经纬度和地理坐标系的转换仿真
%% 输入参数
R=6371;%地球半径，单位km
lng=109.414871;%经度
lat=18.303421;%纬度


%% 算法计算
xt=R*cosd(lat)*cosd(lng);
yt=R*cosd(lat)*sind(lng);
zt=R*sind(lat);

%% 结果展示（绘制三维球体地图）
plot_Globe

%% 结果展示（绘制三维点）
% 绘制三维点图
plot3(xt, yt, zt, 'o', 'MarkerSize', 10, 'MarkerEdgeColor', 'r', 'MarkerFaceColor', 'g');

%% 平面计算
[A,B,C,D]=GPS2EQPlane(lng,lat);
%% 结果展示（绘制三维面）
plot_plane(A,B,C,D, 7000);

% %% 旋转平面
% theta =90; % 旋转角度
% [newA, newB, newC, newD] = rotatePlane([xt yt zt],[A B C], theta);

% %% 结果展示（绘制三维面）
% plot_plane(newA, newB, newC, newD, 7000);

3.1.2 结果展示

3.2 代码二示例

以经度:0、纬度:0为例：

3.1.1 测试代码

clc;clear;close all
%% 经纬度和地理坐标系的转换仿真
%% 输入参数
R=6371;%地球半径，单位km
lng=109.414871;%经度
lat=18.303421;%纬度

%% 算法计算
xt=R*cosd(lat)*cosd(lng);
yt=R*cosd(lat)*sind(lng);
zt=R*sind(lat);

%% 结果展示（绘制三维球体地图）
plot_Globe

%% 结果展示（绘制三维点）
% 绘制三维点图
plot3(xt, yt, zt, 'o', 'MarkerSize', 10, 'MarkerEdgeColor', 'r', 'MarkerFaceColor', 'g');

%% 平面计算
[A,B,C,D]=GPS2EQPlane(lng,lat);
% %% 结果展示（绘制三维面）
% plot_plane(A,B,C,D, 7000);

%% 旋转平面
theta =45; % 旋转角度
[newA, newB, newC, newD] = rotatePlane([xt yt zt],[A B C], theta);

%% 结果展示（绘制三维面）
plot_plane(newA, newB, newC, newD, 7000);

3.1.2 结果展示

旋转45度后

旋转90度后

旋转180度后