大地问题解算.cpp

由大地经纬度转为空间直角坐标

#include<iostream.h>
#include<math.h>
void main()
	{
	//正解
	double s=40000.000,epi2=0.0067395018195,B1=32*3.14159265359/180,ita12,L1=118+10/60+5/3600;
	ita12=epi2*cos(B1)*cos(B1);
	cout.precision(11);
	cout<<"ita1^2\t"<<ita12<<endl;
	double c=6399596.652,N1;
	N1=c/sqrt(1+ita12);
	cout<<"N1\t"<<N1<<endl;
	double A1=45*3.14159265359/180,V1;
	V1=sqrt(1+ita12);
	cout<<"V1\t"<<V1<<endl;
	cout<<"A1\t"<<A1<<endl<<endl;
	double l[3],b[3],a[3],deil,deib,deic,t,m;
	b[0]=s*cos(A1)*V1*V1/N1;
	l[0]=s*sin(A1)/(N1*cos(B1));
	m=l[0]*cos(B1);
	cout<<"第一次赋值: \nb,l,m\t"<<0<<'\t'<<b[0]<<'\t'<<l[0]<<'\t'<<m<<endl;
	for(int i=1;i<3;i++)
		{ 
		t=tan(B1);//B1=0.5(B1+B2)
		cout<<"t\t"<<i-1<<'\t'<<t<<endl<<endl;
		deil=m*m*t*t/24-b[i-1]*b[i-1]*(1+ita12-9*ita12*t*t)/(24*V1*V1*V1*V1);
		deib=m*m*(2+3*t*t+2*ita12)/24+b[i-1]*b[i-1]*ita12*(t*t-1-ita12-4*ita12*t*t)/(8*V1*V1*V1*V1);
		deic=m*m*V1*V1/12+b[i-1]*b[i-1]*(3+8*ita12+5*ita12*ita12)/(24*V1*V1*V1*V1);
		cout<<"deil,deib,deic\t"<<i<<'\t'<<deil<<'\t'<<deib<<'\t'<<deic<<endl;
		l[i]=s*sin(A1)*(1+deil)/(N1*cos(B1));
		b[i]=s*cos(A1)*V1*V1*(1+deib)/N1;
		a[i]=s*sin(A1)*t*(1+deil+deic)/N1;
		cout<<"l,b,a\t"<<l[i]<<'\t'<<b[i]<<'\t'<<a[i]<<endl;
		L1+=(l[i]*180/3.14159265359);
		/*cout<<"L1,B1\t"<<L1<<'\t'<<B1*180/3.14159265359<<endl<<endl<<endl;*/

		//重新赋值
		A1+=(a[i]/2+(90*3.14159265359/180));
		B1+=(b[i]/2);
		m=l[i]*cos(B1);
		cout<<"m\t"<<i<<'\t'<<m<<endl;
		ita12=epi2*cos(B1)*cos(B1);
		cout<<"ita\t"<<i<<'\t'<<ita12<<endl;
		V1=sqrt(1+ita12);
		cout<<"V\t"<<i<<'\t'<<V1<<endl;
		N1=c/sqrt(1+ita12);
		cout<<"N1\t"<<i<<'\t'<<N1<<endl;
		cout<<"A1\t"<<A1<<endl;
		cout<<"B1\t"<<B1<<endl;
		}
	}