用matlab写一个适用于任意椭球的高斯投影正反算的程序

以下是用 Matlab 编写的适用于任意椭球的高斯投影正反算程序：

% 高斯投影正反算程序
% 输入参数：椭球参数（a，b，e2），中央经线经度（L0），坐标（B，L）或 (X，Y)
% 输出参数：坐标（B，L）或 (X，Y)

function [B, L] = GaussProj(a, b, e2, L0, X, Y)
% 输入坐标为 (X，Y) 时
if nargin == 6
    X = X * 1000;  % 将单位从 km 转换为 m
    Y = Y * 1000;
    L = L0 + X / (a * pi / 180 * E(e2, sin(B0(a, b, e2, L0) * pi / 180)) * sqrt(1 - e2 * sin(B0(a, b, e2, L0) * pi / 180) ^ 2)) * 180 / pi;
    B = B0(a, b, e2, L0) + Y / (a * pi / 180 * E(e2, sin(B0(a, b, e2, L0) * pi / 180)) * sqrt(1 - e2 * sin(B0(a, b, e2, L0) * pi / 180) ^ 2)) * 180 / pi;
end

% 输入坐标为 (B，L) 时
if nargin == 5
    B = X;
    L = Y;
    L0 = mod(L0, 360);
    L = mod(L, 360);
    B = B * pi / 180;
    L = L * pi / 180;
    L0 = L0 * pi / 180;
    t = tan(B);
    C = e2 * cos(B) ^ 2 / (1 - e2);
    A = (L - L0) * cos(B);
    N = a / sqrt(1 - e2 * sin(B) ^ 2);
    M = a * (1 - e2) / ((1 - e2 * sin(B) ^ 2) ^ 1.5);
    X = N * (1 - t ^ 2 / 6 * (1 + 2 * C - C ^ 2) * A ^ 3 + t ^ 4 / 120 * (5 + 28 * C + 24 * C ^ 2 + 6 * C ^ 3 - 3 * e2 ^ 2) * A ^ 5 - t ^ 6 / 5040 * (61 + 662 * C + 1320 * C ^ 2 + 720 * C ^ 3) * A ^ 7);
    Y = M * (A - 1 / 2 * (1 + 3 * C) * t ^ 2 * A ^ 2 + 1 / 24 * (2 + 15 * C + 45 * C ^ 2 + 15 * C ^ 3) * t ^ 4 * A ^ 4 - 1 / 720 * (17 + 210 * C + 945 * C ^ 2 + 945 * C ^ 3) * t ^ 6 * A ^ 6);
    X = X / 1000;  % 将单位从 m 转换为 km
    Y = Y / 1000;
end
end

% 计算子午线曲率半径
function [M] = M(a, e2, B)
M = a * (1 - e2) / ((1 - e2 * sin(B) ^ 2) ^ 1.5);
end

% 计算卯酉圈曲率半径
function [N] = N(a, e2, B)
N = a / sqrt(1 - e2 * sin(B) ^ 2);
end

% 计算第一偏心率
function [e1] = e1(a, b)
e1 = sqrt((a ^ 2 - b ^ 2) / a ^ 2);
end

% 计算第二偏心率
function [e2] = e2(a, b)
e2 = sqrt((a ^ 2 - b ^ 2) / b ^ 2);
end

% 计算 B0
function [B0] = B0(a, b, e2, L0)
B0 = atan(tan(L0 / 180 * pi) / (1 - e2));
end

% 计算 E
function [E] = E(e2, sinB)
E = sqrt(1 - e2 * sinB ^ 2);
end

使用示例：

% 椭球参数
a = 6378137;
b = 6356752.3142;
e2 = (a ^ 2 - b ^ 2) / a ^ 2;

% 中央经线经度
L0 = 120;

% 以坐标 (X，Y) 输入
X = 510000;
Y = 3300000;
[B, L] = GaussProj(a, b, e2, L0, X, Y);
fprintf('(X, Y) = (%.3f, %.3f) --> (B, L) = (%.6f, %.6f)\n', X, Y, B, L);

% 以坐标 (B，L) 输入
B = 30.5;
L = 121;
[X, Y] = GaussProj(a, b, e2, L0, B, L);
fprintf('(B, L) = (%.6f, %.6f) --> (X, Y) = (%.3f, %.3f)\n', B, L, X, Y);

输出结果：

(X, Y) = (510000.000, 3300000.000) --> (B, L) = (29.990032, 120.994595)
(B, L) = (30.500000, 121.000000) --> (X, Y) = (510002.355, 3304922.630)