C++实现白贝塞尔法大地主题正反算

// Bessel.cpp : 此文件包含 “main“ 函数。程序执行将在此处开始并结束。
//

#include 
#include 
using namespace std;

#define e 0.081813334016931499   //克拉索夫斯基椭球体第一偏心率

int main（）
{
	int k B10 B11 L10 L11 A10 A11 B20 B21 L20 L21 A20 A21;
	double B12 L12 A12 B22 L22 A22;
	double B1 L1 A1 S B2 L2 A2 L pi;
	double A B C afa beta;
	double a1 a2 b1 b2 p q x y;
	double W1 W2 sinu1 sinu2 cosu1 cosu2 sinA0 cotsigma1 sin2sigma1 cos2sigma1 sigma0 sin2 cos2 sigma sins coss delta0 delta lamda;
	pi = 4 * atan（1）;
	cout << “白塞尔大地主题正算请输入1\n白塞尔大地主题反算请输入2“ << endl;
	cin >> k;
	if （k == 1）
	{
		cout << “请输入大地线起点纬度B经度L，大地方位角A（度，分，秒）及大地线长度S:“ << endl;
		cin>> B10>>B11>>B12>>L10>>L11>>L12>>A10>>A11>>A12>>S;
		B1 = （B10 + （float）B11 / 60 + B12 / 3600） * pi / 180;
		L1 = （L10 + （float）L11 / 60 + L12 / 3600） * pi / 180;
		A1 = （A10 + （float）A11 / 60 + A12 / 3600） * pi / 180;
		W1 = sqrt（1 - e * e * sin（B1） * sin（B1））;            //计算起点规划纬度
		sinu1 = sin（B1） * sqrt（1 - e * e） / W1;              //计算起点规划纬度
		cosu1 = cos（B1） / W1;                          //计算起点规划纬度
		sinA0 = cosu1 * sin（A1）;                                           //计算辅助函数值
		cotsigma1 = cosu1 * cos（A1） / sinu1;                                 //计算辅助函数值
		sin2sigma1 = 2 * cotsigma1 / （cotsigma1 * cotsigma1 + 1）;                //计算辅助函数值
		cos2sigma1 = （cotsigma1 * cotsigma1 - 1） / （cotsigma1 * cotsigma1 + 1）;    //计算辅助函数值
		A = 6356863.020 + （10708.949 - 13.474 * （1 - sinA0 * sinA0）） * （1 - sinA0 * sinA0）;
		B = （5354.469 - 8.798 * （1 - sinA0 * sinA0）） * （1 - sinA0 * sinA0）;
		C = （2.238 * （1 - sinA0 * sinA0）） * （1 - sinA0 * sinA0） + 0.006;
		afa = 691.46768 - （0.58143 - 0.00144 * （1 - sinA0 * sinA0）） * （1 - sinA0 * sinA0）;
		beta = （0.2907 - 1.0E-3 * （1 - sinA0 * sinA0）） * （1 - sinA0 * sinA0）;
		sigma0 = （S - （B + C * cos2sigma1） * sin2sigma1） / A;
		sin2 = sin2sigma1 * cos（2 * sigma0） + cos2sigma1 * sin（2 * sigma0）;
		cos2 = cos2sigma1 * cos（2 * sigma0） - sin2sigma1 * sin（2 * sigma0）;
		sigma = sigma0 + （B + 5 * C * cos2） * sin2 / A;
		delta = （afa * sigma + beta * （sin2 - sin2sigma1）） * sinA0;    //计算经度差改正数
		delta = delta / 3600 * pi / 180;
		sinu2 = sinu1 * cos（sigma） + cosu1 * cos（A1） * sin（sigma）;
		B2 = atan（sinu2 / （sqrt（1 - e * e） * sqrt（1 - sinu2 * sinu2）））;
		lamda = atan（sin（A1） * sin（sigma） / （cosu1 * cos（sigma） - sinu1 * sin（sigma） * cos（A1）））;
		if （sin（A1） > 0）
		{
			if （tan（lamda） > 0）
				lamda = fabs（lamda）;
			else
				lamda = pi - fabs（lamda）;
		}
		else
		{
			if （tan（lamda） > 0）
				lamda = fabs（lamda） - pi;
			else
				lamda = -1 * fabs（lamda）;
		}
		L2 = L1 + lamda - delta;
		A2 = atan（cosu1 * sin（A1） / （cosu1 * cos（sigma） * cos（A1） - sinu1 * sin（sigma）））;
		if （sin（A1） > 0）
		{
			if （tan（A2） > 0）
				A2 = pi + fabs（A2）;
			else
				A2 = 2 * pi - fabs（A2）;
		}
		else
		{
			if （tan（A2） > 0）
				A2 = fabs（A2）;
			else
				A2 = pi - fabs（A2）;
		}
		B2 = B2