大地线正反算

// ConsoleApplication1.cpp : 定义控制台应用程序的入口点。
//
#include“stdafx.h“
#include
#include“string.h“
#include“math.h“
#include 
#include  //要用到格式控制符
using namespace std;
const   double   PI = 3.1415926535898;
double a = 6377397.155 f = 299.15281285;
int j=1;//判断当前参考椭球，默认为白塞尔椭球
void ChangeAF（）
{
	if （j == 1）
		cout << “Now : Ellipsoid Bessel 1841“ << endl;
	else if （j == 2）
		cout << “Now : BJ54“;
	cout << “Choose objective Coordinate System：（1：Bessel1841）（ 2：BJ54）“<		cin >> j;
		if （j == 1）
		{
			a = 6377397.155; f = 299.15281285;
		}
		else if （j == 2）
			a = 6378245; f = 298.3;
	
}

double DMS2RAD（double dmsAngle）
{
	int degAngle minAngle sign;
	double radAngle secAngle;
	sign = 1;
	if （dmsAngle < 0）
	{
		sign = -1;
		dmsAngle = fabs（dmsAngle）;
	}
	degAngle = （int）（dmsAngle + 0.0001）;
	minAngle = （int）（（dmsAngle - degAngle） * 100 + 0.0001）;
	secAngle = （dmsAngle - degAngle - minAngle / 100.0）*10000.0;
	radAngle = （degAngle + minAngle / 60.0 + secAngle / 3600.0）*PI/ 180.0;
	radAngle = radAngle * sign;
	return radAngle;
}
double RAD2DMS（double radAngle）
{
	int degAngle minAngle sign;
	double secAngle dmsAngle;
	sign = 1;
	if （radAngle < 0）
	{
		sign = -1;
		radAngle = fabs（radAngle）;
	}
	secAngle = radAngle*180.0 / PI*3600.0;
	degAngle = （int）（secAngle / 3600 + 0.0001）;
	minAngle = （int）（（secAngle - degAngle*3600.0） / 60.0 + 0.0001）;
	secAngle = secAngle - degAngle*3600.0 - minAngle*60.0;
	if （secAngle < 0）  secAngle = 0;
	dmsAngle = degAngle + minAngle / 100.0 + secAngle / 10000.0;
	dmsAngle = dmsAngle * sign;
	return dmsAngle;
}
void Direct（int ndouble Bfromdouble Lfromdouble Afromdouble S）
{
	
	double  Bi Li e1bV N M deltaS geodesicC dB dL Btemp Vtemp;
	b = a - a/f;	
	Bi = DMS2RAD（Bfrom）;
	Li= DMS2RAD（Lfrom）;
	Afrom = DMS2RAD（Afrom）;
	e1 = sqrt（（a*a - b*b） / b/b）;
	V = sqrt（1 + e1*e1*cos（Bi）*cos（Bi））;		
	N = a*a / b / V;
	M = N / V / V;
	deltaS = S / n;
	geodesicC = sin（Afrom）*N*cos（Bi）;

	for （int i = 1; i <=n; i++）
	{
		dB = deltaS*cos（Afrom） / M;
		Btemp = Bi + dB / 2.0;
		Vtemp = sqrt（1 + e1*e1*cos（Btemp）*cos（Btemp））;
		N = a*a / b / Vtemp;
		
		Afrom = asin（geodesicC / N / cos（Btemp））;
		M = a*a / b / Vtemp / Vtemp / Vtemp;

		dL = deltaS*sin（Afrom） / （N*cos（Btemp））;
		dB = deltaS*cos（Afrom） / M;
		Bi = Bi + dB;
		Li = Li + dL;
		Vtemp = sqrt（1 + e1*e1*cos（Bi）*cos（Bi））;
		M = a*a / b / Vtemp / Vtemp / Vtemp;
		N = a*a / b / Vtemp;
		Afrom = asin（geodesicC / N / cos（Bi））;
	}
	cout.setf（ios::fixed）;//不带指数