算法原理：
高斯平均引数正算公式推导的基本思想是：首先把勒让德级数在P1点展开改在大地线长度中点M展开，以使级数公式项数减少，收敛快，精度高；其次，考虑到求定点中M的复杂性，将M点用大地线两端点平均纬度及平均方位角相对应的m点来代替，并借助迭代运算，便可顺利实现大地主题正解。
算法流程：
首先输入大底线起点纬度B1，经度，大地方位角和大底线长度。其中经纬度与大地方位角角度弧度转换函数转换为弧度再参与计算。而后给定椭球参数（本次使用克氏椭球），给定Bm，Am初值，而后参与迭代计算。当两次迭代所得结果相差小于限值（此处考虑到计算机会省略小数点后几位，故程序中直接另两者相等）时，结束迭代。而后判断A21所在象限，最后输出结果。
#include <iostream>#include<cmath>using namespace std;double angle_to_radian(double degree, double min, double second){double flag = (degree < 0) ? -1.0 : 1.0;//判断正负double pi = 4 * atan(1);if (degree < 0){degree = degree * (-1.0);}double angle = degree + min / 60 + second / 3600;double result = flag * (angle * pi) / 180;return result;//cout<<result<<endl;eg}void radian_to_angle(double rad){double flag = (rad < 0) ? -1.0 : 1.0;double pi = 4 * atan(1);if (rad < 0){rad = rad * (-1.0);}double result = (rad * 180) / pi;double degree = int(result);double min = (result - degree) * 60;double second = (min - int(min)) * 60;cout << flag * degree << "°" << int(min) << "′" << second << "″";}int main(){double angle_to_radian(double degree, double min, double second);void radian_to_angle(double rad);double B1, L1, A12, B1d, B1f, B1m, L1d, L1f, L1m, A12d, A12f, A12m, S;cout << "请输入大底线起点纬度B1=" << "°" << "′" << "″" << endl;cin >> B1d >> B1f >> B1m;cout << "经度L1=" << "°" << "′" << "″" << endl;cin >> L1d >> L1f >> L1m;cout << "大地方位角A1=" << "°" << "′" << "″" << endl;cin >> A12d >> A12f >> A12m;cout << "大底线长度S=";cin >> S;cout << endl;B1 = angle_to_radian(B1d, B1f, B1m);L1 = angle_to_radian(L1d, L1f, L1m);A12 = angle_to_radian(A12d, A12f, A12m);//给定初值double p = 206265;double e2 = 0.006693421622966, e_2 = 0.006738525414683, c = 6399698.9017827110;//克拉索夫斯基椭球double B2 = B1;double A21 = A12;double L2;double deltaL = 0;double deltaB = 0;double deltaA = 0;double Lt, Bt, At;double i = 0;do{double Bm = (B1 + B2) / 2;double Am = (A12 + A21) / 2;double Vm = sqrt(1 + e_2 * cos(Bm) * cos(Bm));double Mm = c / (Vm * Vm * Vm), Nm = Mm * Vm * Vm;double deltaB0 = S * cos(Am) / Mm;double deltaL0 = S * sin(Am) * (1 / cos(Bm)) / Nm;double deltaAm = S * sin(Am) * tan(Bm) / Nm;Lt = deltaL, Bt = deltaB, At = deltaA;deltaL = S * sin(Am) * (1 / cos(Bm)) * (1 + deltaAm * deltaAm / (24) - deltaB0 * deltaB0 / (24)) / Nm;deltaB = S * cos(Am) * (1 + deltaL0 * deltaL0 / (12) + deltaAm * deltaAm / (24)) / Mm;deltaA = S * sin(Am) * tan(Bm) * (1 + deltaB0 * deltaB0 / (12) + deltaL0 * deltaL0 * cos(Bm) * cos(Bm) / (12) + deltaAm * deltaAm / (24)) / Nm;B2 = B1 + deltaB;L2 = L1 + deltaL;A21 = A12 + deltaA;i = i + 1;cout << "第" << i << "次迭代" << endl;cout << "B2=";radian_to_angle(B2);cout << endl << "L2=";radian_to_angle(L2);cout << endl << "A21=";radian_to_angle(A21);cout << endl;} while (Bt != deltaB);double pi = 4 * atan(1);if (A12 > pi){A21 = A12 + deltaA - pi;}else{A21 = A12 + deltaA + pi;}cout << "最终结果" << endl;cout << "B2=";radian_to_angle(B2);cout << endl << "L2=";radian_to_angle(L2);cout << endl << "A21=";radian_to_angle(A21);}
