#include "stdafx.h"
#include "ZuoBiao.h"
#include <math.h>
#include "stdlib.h" 
//WGS-84椭球体参数 

const double a = 6378137.0;//长半轴 

const double flattening = 1 / 298.257223563;//扁率 

const double delta = 0.0000001;

// 两个数的差的平方，坐标差的平方。
inline double diffsqr(double x1, double x2){
	return (x1 - x2)* (x1 - x2);
}

//直接坐标表示，求距离 
inline double distance(double x1, double y1, double z1
	, double x2, double y2, double z2)
{
	return sqrt(diffsqr(x1, x2) + diffsqr(y1, y2) + diffsqr(z1, z2));
}
//点表示，求距离 
double distance(const CRDCARTESIAN &p1, const CRDCARTESIAN &p2)
{
	return sqrt(diffsqr(p1.x, p2.x) + diffsqr(p1.y, p2.y) + diffsqr(p1.z, p2.z));
}

//pcc：指向所转换出的笛卡尔坐标的指针； 
//pcg：指向待转换的大地坐标的指针； 
//dSemiMajorAxis：参考椭球的长半轴； 
//dFlattening：参考椭球的扁率。 
//由大地坐标转换为笛卡尔坐标 
void GeodeticToCartesian(PCRDCARTESIAN pcc, PCRDGEODETIC pcg, double dSemiMajorAxis, double dFlattening)
{
	double e2;//第一偏心率的平方 
	double N;//卯酉圈半径 
	e2 = 2 * dFlattening - dFlattening*dFlattening;
	N = dSemiMajorAxis / sqrt(1 - e2*sin(pcg->latitude)*sin(pcg->latitude));
	pcc->x = (N + pcg->height)*cos(pcg->latitude)*cos(pcg->longitude);
	pcc->y = (N + pcg->height)*cos(pcg->latitude)*sin(pcg->longitude);
	pcc->z = (N*(1 - e2) + pcg->height)*sin(pcg->latitude);
}

//获取两个经纬度之间距离
//经度
//维度
double GetJuLi(double longitude1, double latitude1, double longitude2, double latitude2)
{
	double dRet = 0;
	tagCRDGEODETIC aa;
	aa.height = 0;
	aa.longitude = longitude1;
	aa.latitude = latitude1;

	CRDCARTESIAN bb;
	bb.x = 0;
	bb.y = 0;
	bb.z = 0;

	tagCRDGEODETIC cc;
	cc.height = 0;
	cc.longitude = longitude2;
	cc.latitude = latitude2;

	CRDCARTESIAN dd;
	dd.x = 0;
	dd.y = 0;
	bb.z = 0;

	try
	{
		GeodeticToCartesian(&bb, &aa, a, flattening);

		GeodeticToCartesian(&dd, &cc, a, flattening);

		dRet = distance(bb, dd);
	}
	catch (...)
	{
		dRet = 100;
	}

	return dRet;
}