算法原理
高斯投影是将椭球体上的元素投影到平面上的一种方式，是正形投影的一种。将地球椭球以一定经差范围分带（一般为6°带和3°带。本次使用6°带），然后分带投影。在投影面上，中央子午线和赤道的投影都是直线，并且以中央子午线和赤道的交点作为坐标原点，以中央子午线的投影为纵坐标轴，以赤道的投影为横坐标轴，便形成了高斯平面直角坐标系。
算法流程
首先输入待转换点纬度与经度，然后计算带号和待转换点与轴子午线的经差，再给定椭球参数（此处使用克氏椭球和1975国际椭球两种椭球参数），然后调用转换函数输出结果。转换函数中，先计算辅助值V，N，t,ղ，然后计算自赤道量起的子午线弧长X，最后代入公式计算x,y，输出结果。
代码
#include <iostream>#include<cmath>using namespace std;int main(){//输入已知值double angle_to_radian(double degree, double min, double second);void radian_to_angle(double rad);void Guass(double e2, double c, double B, double l);double B, L, Bd, Bf, Bm, Ld, Lf, Lm;cout << "请输入待转换点纬度B1=" << "°" << "′" << "″" << endl;cin >> Bd >> Bf >> Bm;cout << "经度L1=" << "°" << "′" << "″" << endl;cin >> Ld >> Lf >> Lm;B = angle_to_radian(Bd, Bf, Bm);L = angle_to_radian(Ld, Lf, Lm);//计算lint n = Ld / 6 + 1;int L0 = 6 * n - 3;double ld, lf, lm;/* double L0 = 111;double ld=3, lf=20, lm=0;*/if (Ld > L0){ld = Ld - L0;lf = Lf;lm = Lm;}else{if(Lm==0){lm = 0;lf = 60 - Lf;}else{lm = 60 - Lm;lf = 60 - 1.0 - Lf;}ld = Ld - L0 + 1;}double l = angle_to_radian(ld, lf, lm);//克氏椭球参数double long e2 = 0.006738525414683;double long c = 6399698.9017827110;//double X1 = 111134.861 * Bd - 16036.480 * sin(2 * B) + 16.828 * sin(4 * B) - 0.022 * sin(6 * B);//1975国际椭球参数double long e2_ = 0.006739501819473;double long c_ = 6399596.6519880105;cout << "按照克拉索夫斯基椭球计算所得结果为：" << endl;Guass(e2, c, B, l);cout << "按照1975年国际椭球计算所得结果为：" << endl;Guass(e2_, c_, B, l);}void Guass(double e2, double c, double B, double l){//辅助值double V = sqrt(1 + e2 * cos(B) * cos(B));double N = c / V;double t = tan(B);double ita = sqrt(e2 * cos(B) * cos(B));//计算Xdouble b0 = 1 - 3* e2 / 4 + 45* e2 * e2 / 64 - 175 * e2 * e2 * e2 / 256 + 11025 * e2 * e2 * e2 * e2 / 16384;double b2 = b0 - 1;double b4 = 15 * e2 * e2 / 32 - 175 * e2 * e2 * e2 / 384 + 3675 * e2 * e2 * e2 * e2 / 8192;double b6 = -35 * e2 * e2 * e2 / 96 + 735 * e2 * e2 * e2 * e2 / 2048;double b8 = 315 * e2 * e2 * e2 * e2 / 1024;double X = c * (b0 * B + (b2 * cos(B) + b4 * cos(B) * cos(B) * cos(B) + b6 * cos(B) * cos(B) * cos(B) * cos(B) * cos(B) + b8 * cos(B) * cos(B) * cos(B) * cos(B) * cos(B) * cos(B) * cos(B)) * sin(B));//计算x,ydouble x = X + N * sin(B) * cos(B) * l * l / 2 + N * sin(B) * cos(B) * cos(B) * cos(B) * (5 - t * t + 9 * ita * ita + 4 * ita * ita * ita * ita) * l * l * l * l / 24 +N * sin(B) * cos(B) * cos(B) * cos(B) * cos(B) * cos(B) * (6l - 58 * t * t + t * t * t * t) * l * l * l * l * l * l / 720;double y = N * cos(B) * l + N * cos(B) * cos(B) * cos(B) * (1 - t * t + ita * ita) * l * l * l / 6 + N * cos(B) * cos(B) * cos(B) * cos(B) * cos(B) * (5 - 18 * t * t + t * t * t * t + 14 * ita * ita - 58 * ita * ita * t * t) * l * l * l * l * l / 120;/* cout << "V=" << V << endl << "N=" << N << endl << "t=" << t << endl << "ita=" << ita << endl << "X=" << X << endl;*/cout << "x=" << x << endl << "y=" << y << endl;}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 << "″";}
