高斯大地主题解算大地测量程序

// 高斯正大地主题解算.cpp : 定义控制台应用程序的入口点。
//

#include "stdafx.h"
#include "iostream"
#include "math.h"
#define PI 3.1415926
#define e2 0.00673949674227 //采用WGS-84椭球体数据
#define ee2 0.00669437999013
#define c 6399593.6258
#define a 6378137
#define p 206264.806247096355//弧度值转换位秒值
using namespace std;


//度分秒转化为度
double ZH(double du,double fen,double miao)
{
double HUDU;
if(du>0)
{
HUDU=(miao/60.0+fen)/60.0+du;
}
else
{
HUDU=-((miao/60.0+fen)/60.0+abs(du));
}//绝对值是防止输入的是负值角度
return HUDU;
}



int _tmain(int argc, _TCHAR* argv[])
{
//判断正反算
int k;
cout<<"请选择：1---正算 2---反算"<cin>>k;


double B1,L1,B2,L2,A1,A2,S;
double E[3],F[3],B[3],L[3],aa[3];//B,L存放B1和L1、E,F存放B2，L2、aa存放方位角A1



if(k==1)//正算
{


//用户输入已知值


cout<<"请输入大地线起点纬度B1："<cin>>B[0]>>B[1]>>B[2];
B1=ZH(B[0],B[1],B[2])*PI/180.0;//已经化为弧度值


cout<<"请输入大地线起点经度L1："<cin>>L[0]>>L[1]>>L[2];
L1=ZH(L[0],L[1],L[2])*PI/180.0;//已经化为弧度值

cout<<"请输入大地线大地方位角a："<cin>>aa[0]>>aa[1]>>aa[2];
A1=ZH(aa[0],aa[1],aa[2])*PI/180.0;//已经化为弧度值

cout<<"请输入大地线长度S："<cin>>S;


//计算起点的各项值
double W1=sqrt(1-(e2)*sin(B1)*sin(B1));

double V1=sqrt(1-(ee2)*cos(B1)*cos(B1));

double M1=a*(1-e2)/(W1*W1*W1);

double N1=a/W1;

//计算主项
double oB0=S*cos(A1)/M1; //弧度值

double oL0=S*sin(A1)*(1/cos(B1))/N1;//弧度制

double oA0=oL0*sin(B1);//弧度值

double Bm=B1+0.5*oB0;

double Am=A1+0.5*oA0;


double Vm,Nm,tm,nm;
int i=0;
double oB,oL,oA;


do{
//计算中点值



Vm=sqrt(1-(ee2)*cos(Bm)*cos(Bm));

Nm=a/( sqrt(1-(e2)*sin(Bm)*sin(Bm)) );

tm=tan(Bm);

nm=ee2*cos(Bm)*cos(Bm);




oB=(Vm*Vm)/Nm*cos(Am)*( 1 + S*S/(24*Nm*Nm)*( sin(Am)*sin(Am)*(2+3*tm*tm+3*nm*nm*tm*tm)+3*nm*nm*cos(Am)*cos(Am)*( -1+tm*tm-nm*nm-4*tm*tm*nm*nm ) ) );

oL=(1/Nm)*S*(1/cos(Bm))*sin(Am)*( 1+S*S/(24*Nm*Nm)*( sin(Am)*sin(Am)*tm*tm-cos(Am)*cos(Am)*( 1+nm*nm-9*tm*tm*nm*nm+nm*nm*nm*nm ) ) );

oA=(1/Nm)*S*sin(Am)*tm*( 1+S*S/(24*Nm*Nm)*( cos(Am)*cos(Am)*( 2+7*nm*nm+9*tm*tm*nm*nm+5*nm*nm*nm*nm ) +sin(Am)*sin(Am)*( 2+tm*tm+2*nm*nm ) ) );


Bm=B1+0.5*oB;
Am=A1+0.5*oA;

i++;

}while(i<4);//迭代计算3次


B2=(B1+oB)*180/PI;

L2=(L1+oL)*180/PI;

A2=(A1+oA+PI)*180/PI;
//输出结果
cout<<"终点大地经度B2："<<(int)B2<<"度"<cout<<"终点大地纬度L2："<<(int)L2<<"度"<

(L2-(int)L2)*60))<<"分"<cout<<"终点大地方位角A2："<<(int)A2<<"度"<
system("pause");

}

else//反算
{

//用户输入
cout<<"请输入大地线起点纬度B1："<cin>>B[0]>>B[1]>>B[2];
B1=ZH(B[0],B[1],B[2])*PI/180.0;

cout<<"请输入大地线起点经度L1："<cin>>L[0]>>L[1]>>L[2];
L1=ZH(L[0],L[1],L[2])*PI/180.0;

cout<<"请输入大地线终点纬度B2："<cin>>E[0]>>E[1]>>E[2];
B2=ZH(E[0],E[1],E[2])*PI/180.0;

cout<<"请输入大地线终点纬度L2："<cin>>F[0]>>F[1]>>F[2];
L2=ZH(F[0],F[1],F[2])*PI/180.0;

//辅助计算
double Bm=(B1+B2)/2;

double oB=B2-B1;

double oL=L2-L1;

double Wm=sqrt(1-(e2)*sin(Bm)*sin(Bm));

double Vm=sqrt(1-(ee2)*cos(Bm)*cos(Bm));

double Mm=a*(1-e2)/(Wm*Wm*Wm);

double Nm=a/Wm;

double tm=tan(Bm);

double nm=ee2*cos(Bm)*cos(Bm);

double U=oL*Nm*cos(Bm)+oL*oL*Nm*cos(Bm)/(24*Vm*Vm*Vm*Vm)*(1+nm*nm-9*nm*nm*tm*tm*nm*nm*nm*nm)+oL*oL*oL*(-Nm*cos(Bm)*cos(Bm)*cos(Bm)*tm*tm/24);

double V=oB*Nm/(Vm*Vm)+oB*oL*oL*Nm/(24*Vm*Vm)*cos(Bm)*cos(Bm)*(2+3*tm*tm+2*nm*nm)+oB*oB*oB*Nm/(8*Vm*Vm*Vm*Vm*Vm*Vm)*(nm*nm-tm*tm*nm*nm);

double oAm=tm*cos(Bm)*oL+oB*oB*oL*1/(24*Vm*Vm*Vm*Vm)*cos(Bm)*tm*(2+7*nm*nm+9*tm*tm*nm*nm+5*nm*nm*nm*nm)+oL*oL*oL*1/24*cos(Bm)*cos(Bm)*cos(Bm)*tm*(2+tm*tm+2*nm*nm);

double Am=atan(U/V);


double T;
if(oB>oL)
{
T=atan(abs(U/V));
}
else
{
double cc=abs(U/V);

T=PI/4+atan((1-cc)/(1+cc));
}


if(oB>0 && oL>=0)
{
Am=T;
}
else
{

if(oB<0 && oL>=0)
{
Am=PI-T;
}
else
{
if(oB<=0 && oL<0)
{
Am=PI+T;
}
else
{
if(oB>0 && oL<0 )
{
Am=2*PI-T;
}
else
{
Am=PI/2;
}
}
}

}

S=U/(sin(Am));

A1=(Am-oAm/2)*180/PI;

A2=(Am+oAm/2+PI)*180/PI;


//输出结果
cout<<"大地线长度："<cout<<"起点大地方位角："<<(int)A1<<"度"<cout<<"终点大地方位角："<<(int)A2<<"度"<system("pause");

}

return 0;


}