高斯平面坐标与大地坐标转换

• 包含两部分：
  • 高斯平面坐标xy转为大地坐标BL。
  • 大地坐标BL转为高斯平面坐标xy。
• 用C++语言编写，能够读取txt坐标文件，并将结果导出到新的txt文件。
• 椭球模型为WGS84，可以根据实际需要更改参数。
• 不足之处请指正！

张华海.应用大地测量学(第4版).中国矿业大学出版社.2016

1. 高斯平面坐标xy转为大地坐标BL

convert.h

#pragma once
#include<math.h>
#include<iomanip>
#include <cmath>
using namespace std;

#define PI  3.14159265358979323846;

typedef struct
{
    double pos1;
    double pos2;
}position;

//WGS84椭球参数
double a = 6378137.0;
double f = 1 / 298.257223563;

double e2 = 2 * f - f * f;   //第一偏心率平方
double e_2 = e2 * (1 + e2);  //第二偏心率平方
double c = a / sqrt(1 - e2);
double beita0 = 1 - 3 * e_2 / 4 + 45 * e_2 * e_2 / 64 - 175 * e_2 * e_2 * e_2 / 256 + 11025 * e_2 * e_2 * e_2 * e_2 / 16384;
double beita2 = beita0 - 1;
double beita4 = 15 * e_2 * e_2 / 32 - 175 * e_2 * e_2 * e_2 / 384 + 3675 * e_2 * e_2 * e_2 * e_2 / 8192;
double beita6 = -35 * e_2 * e_2 * e_2 / 96 + 735 * e_2 * e_2 * e_2 * e_2 / 2048;
double beita8 = 315 * e_2 * e_2 * e_2 * e_2 / 1024;

double FxB(double B)
{
    double temp = (c * beita4 + (c * beita6 + c * beita8 * cos(B) * cos(B)) * cos(B) * cos(B)) * cos(B) * cos(B);
    return c * beita2 + temp * sin(B) * cos(B);
}

double FxBl(double B, double l)
{
    double W = sqrt(1 - e2 * sin(B) * sin(B));
    double N = a / W;  //卯酉圈曲率半径
    double t2 = tan(B) * tan(B);
    double yita2 = e_2 * cos(B) * cos(B);

    double a2 = N * sin(B) * cos(B) / 2;
    double a4 = N * sin(B) * cos(B) * cos(B) * cos(B) * (5 - t2 + 9 * yita2 + 4 * yita2 * yita2) / 24;
    double a6 = N * sin(B) * cos(B) * cos(B) * cos(B) * cos(B) * cos(B) * (61 - 58 * t2 + t2 * t2) / 720;

    return a2 * l * l + a4 * l * l * l * l + a6 * l * l * l * l * l *