//大橢圓距離公式
#ifndef DATUOYUANJULI
#define DATUOYUANJULI
#include"jingweidu.h"
#include"jingweiduzhuanhua.h"
#include"tuoqiu.h"
#include<math.h>
#define pi 3.141592653589793
double datuoyuanjuli(dadijingweidu da1,dadijingweidu da2,tuoqiu tuo)
{
	double a=tuo.changzhou;
	double f=tuo.bianlv;
    dixinjingweidu di1,di2;
	di1=daditodixin(da1,beijing);//大地經緯度轉化到地心經緯度
	di2=daditodixin(da2,beijing);
	double s;//大橢圓距離
	double theta1,theta2,bi,m,r,deltaL;//中間變量
    if(da1.B==0&&da2.B==0)
	{
		s=a*fabs(da1.L-da2.L);
		return s;
	}
	else if(da1.B==0&&da2.B!=0)
		da1.B=SGN(da2.B)*pi*0.0001/3600/180;
	else if(da1.B!=0&&da2.B==0)
		da2.B=SGN(da1.B)*pi*0.0001/3600/180;
	deltaL=da2.L-da1.L;
	if(deltaL)
	{
		double alpha1,alpha2;
		double fenzi1=sin(deltaL);
		double fenmu1=cos(di1.psi)*(tan(di2.psi)-tan(di1.psi)*cos(deltaL));
		if(fenzi1>0&&fenmu1>0)
			alpha1=atan(fenzi1/fenmu1);
		else if(fenzi1<0&&fenmu1>0)
			alpha1=2*pi+atan(fenzi1/fenmu1);
		else 
			alpha1=pi+atan(fenzi1/fenmu1);

        double fenzi2=-sin(deltaL);
		double fenmu2=cos(di2.psi)*(tan(di1.psi)-tan(di2.psi)*cos(deltaL));
		if(fenzi2>0&&fenmu2>0)
			alpha2=atan(fenzi2/fenmu2);
		else if(fenzi2<0&&fenmu2>0)
			alpha2=2*pi+atan(fenzi2/fenmu2);
		else 
			alpha2=pi+atan(fenzi2/fenmu2);
		theta1=atan(cos(alpha1)/tan(di1.psi));
		theta2=atan(-cos(alpha2)/tan(di2.psi));
		bi=a*(1-f)/sqrt(1-(2*f-square(f))*square(cos(di1.psi))*square(sin(alpha1)));
	}
	else 
	{
		bi=a-a*f;
		theta1=pi/2-di1.psi;
		theta2=pi/2-di2.psi;
	}
	m=(square(a)-square(bi))/(square(a)+square(bi));
	r=sqrt((square(a)+square(bi))/2)*(1-square(m)/16);
	s=r*fabs(theta1-theta2-m/4*(sin(2*theta1)-sin(2*theta2))
		-square(m)/64*(sin(4*theta1)-sin(4*theta2)));
	if(SGN(da1.B)!=SGN(da2.B))
		s=pi*r-s;
	return s;
}
#endif