大地坐标与空间直角坐标系的相互转化

//大地坐标转化为空间直角坐标
#include <iostream>
#include <fstream>
#include <string>
#include<iomanip>
#include<cmath>

using namespace std;
#define PI 3.14159265358979


double angle_to_radian(double degree, double min, double second)
{
	double angle = degree + min / 60 + second / 3600;
	double result = (angle * PI) / 180;
	return result;
}

double radian_to_angle(double L,double B)
{
	double result = (L * 180) / PI;
	double degree = int(result);
	double min = (result - degree) * 60;
	double second = (min - int(min)) * 60;

	double result1 = (B * 180) / PI;
	double degree1 = int(result1);
	double min1 = (result1 - degree1) * 60;
	double second1 = (min1 - int(min1)) * 60;
	
	cout <<"L="<< int(degree) << "o" << int(min) << "'" << fixed << setprecision(6) << second << "''" << "\n";
	cout <<"B="<< int(degree1) << "o" << int(min1) << "'" << fixed << setprecision(6) << second1 << "''" << "\n";
	return result;
}

double  LBHtoXYZ(double L, double B, double H, double& X, double& Y, double& Z)

{
	int modle=-1;
	while (modle != 0)
	{
		cout << "请选择椭球参数：[1]GRS80椭球 [2]GRS75椭球 [3]克拉索夫斯基椭球" << endl;
		cin >> modle;
		
		switch (modle)
		{

		case 1://使用GRS80椭球
		{
			double a = 6378137.0;
			double e2 = 0.00669438002290;
			//转换为弧度   
			

			double WW = 1 - e2 * sin(B) * sin(B);
			double N = a / sqrt(WW); //卯酉圈曲率半径
			double Daita_h = 0;      //高程异常,默认为0
			H = Daita_h + H;

			X = (N + H) * cos(B) * cos(L);
			Y = (N + H) * cos(B) * sin(L);
			Z = (N * (1 - e2) + H) * sin(B);
			cout << fixed << setprecision(6) << X << "\n" << fixed << setprecision(6) << Y << "\n" << fixed << setprecision(6) << Z << endl;
			
			break;
		}

		case 2: //使用GRS75椭球
		{
			double a = 6378140.0;
			double e2 = 0.00669438499959;

			double fac1 = 1 - e2 * sin(B) * sin(B);
			double N = a / sqrt(fac1); //卯酉圈曲率半径
			double Daita_h = 0;      //高程异常,默认为0
			H = Daita_h + H;

			X = (N + H) * cos(B) * cos(L);
			Y = (N + H) * cos(B) * sin(L);
			Z = (N * (1 - e2) + H) * sin(B);
			cout << fixed << setprecision(6) << X << "\n" << fixed << setprecision(6) << Y << "\n" << fixed << setprecision(6) << Z << endl;
			
			break;
		}
		case 3: //使用克拉索夫斯基椭球
		{
			double a = 6378245.0;
			double e2 = 0.00669342162297;

			double fac1 = 1 - e2 * sin(B) * sin(B);
			double N = a / sqrt(fac1); //卯酉圈曲率半径
			double Daita_h = 0;      //高程异常,默认为0
			H = Daita_h + H;

			X = (N + H) * cos(B) * cos(L);
			Y = (N + H) * cos(B) * sin(L);
			Z = (N * (1 - e2) + H) * sin(B);
			cout << fixed << setprecision(6) << X << "\n" << fixed << setprecision(6) << Y << "\n" << fixed << setprecision(6) << Z << endl;
			
			break;
		}
		}
		
	}      return 0;
}

double  XYZtoLBH(double X, double Y, double Z, double& L, double& B, double& H)
{
	int modle = -1;
	while (modle != 0) {
		cout << "请选择椭球参数：[1]GRS80椭球 [2]GRS75椭球 [3]克拉索夫斯基椭球" << endl;
		cin >> modle;
		switch (modle)
		{
		case 1:
		{double a = 6378137.0;
		double e2 = 0.00669438002290;


		double tanB1 = Z / (sqrt(X * X + Y * Y));
		double tanB = 1 / (sqrt(X * X + Y * Y)) * (Z + a * e2 * tanB1 / (sqrt(1 + tanB1 * tanB1 - e2 * tanB1 * tanB1)));
		double Delta = abs(tanB - tanB1);

		while (a != 0)
		{
			if (Delta >= (5e-10))
			{
				tanB1 = tanB;
				tanB = (1 / (sqrt(X * X + Y * Y))) * (Z + (a * e2 * tanB1 / sqrt(1 + tanB1 * tanB1 - e2 * tanB1 * tanB1)));
				Delta = abs(tanB - tanB1);
			}
			else
			{
				break;
			}
		}
		B = atan(tanB);
		double N = a / sqrt(1 - e2 * sin(B) * sin(B)); //卯酉圈曲率半径
		H = (sqrt(X * X + Y * Y) / cos(B)) - N;
		L = atan(Y / X);
		radian_to_angle(L, B);
		cout << "H=" << fixed << setprecision(6) << H;
		break;
		}
		case 2:
		{double a = 6378140.0;
		double e2 = 0.00669438499959;


		double tanB1 = Z / (sqrt(X * X + Y * Y));
		double tanB = 1 / (sqrt(X * X + Y * Y)) * (Z + a * e2 * tanB1 / (sqrt(1 + tanB1 * tanB1 - e2 * tanB1 * tanB1)));
		double Delta = abs(tanB - tanB1);

		while (a != 0)
		{
			if (Delta >= (5e-10))
			{
				tanB1 = tanB;
				tanB = (1 / (sqrt(X * X + Y * Y))) * (Z + (a * e2 * tanB1 / sqrt(1 + tanB1 * tanB1 - e2 * tanB1 * tanB1)));
				Delta = abs(tanB - tanB1);
			}
			else
			{
				break;
			}
		}
		B = atan(tanB);
		double N = a / sqrt(1 - e2 * sin(B) * sin(B)); //卯酉圈曲率半径
		H = (sqrt(X * X + Y * Y) / cos(B)) - N;
		L = atan(Y / X); 
		radian_to_angle(L, B);
		cout << "H=" << fixed << setprecision(6) << H; 
		break;
		}
		case 3:
		{double a = 6378245.0;
		double e2 = 0.00669342162297;

		double tanB1 = Z / (sqrt(X * X + Y * Y));
		double tanB = 1 / (sqrt(X * X + Y * Y)) * (Z + a * e2 * tanB1 / (sqrt(1 + tanB1 * tanB1 - e2 * tanB1 * tanB1)));
		double Delta = abs(tanB - tanB1);

		while (a != 0)
		{
			if (Delta >= (5e-10))
			{
				tanB1 = tanB;
				tanB = (1 / (sqrt(X * X + Y * Y))) * (Z + (a * e2 * tanB1 / sqrt(1 + tanB1 * tanB1 - e2 * tanB1 * tanB1)));
				Delta = abs(tanB - tanB1);
			}
			else
			{
				break;
			}
		}
		B = atan(tanB);
		double N = a / sqrt(1 - e2 * sin(B) * sin(B)); //卯酉圈曲率半径
		H = (sqrt(X * X + Y * Y) / cos(B)) - N;
		L = atan(Y / X); }
		radian_to_angle(L, B);
		cout << "H=" << fixed << setprecision(6) << H;
		break;
		}

		return 0;
		
	}
}



int main()
{
	double L, B, H, X, Y, Z;
	double Ld, Lm, Ls, Bd, Bm, Bs;
	int transform = 0;
	cout << "请选择   [1]大地坐标转大地空间直角坐标  [2]大地空间直角坐标转大地坐标" << endl;
	cin >> transform;
	switch (transform)
	{
	case 1:
	{cout << "请输入\n L:\n B:\n H：\n(上标用空格表示)" << endl;
	cin >> Ld >> Lm >> Ls;
	L = angle_to_radian(Ld, Lm, Ls);
	cin >> Bd >> Bm >> Bs;
	B = angle_to_radian(Bd, Bm, Bs);
	cin >> H;
	LBHtoXYZ(L, B, H, X, Y, Z);
	}

	case 2:
	{cout << "请输入X: Y: Z：" << endl;
	cin >> X >> Y >> Z;
	XYZtoLBH(X, Y, Z, L, B, H);
	}
	}
}