子午线弧长正反解程序－地图投影基础

子午线弧长正反解是地图数学基础解算重要内容之一，本程序以地图投影类为基础，以类函数的形式提供子午线弧长正反解的实现。该代码在VS2010上测试，计算精度为1E10-8米，1E10-9度。

file *.h

class CMapProjection
{
public:
CMapProjection(void);        // 构造函数
~CMapProjection(void);      //析构函数  
protected:
double mSemiAxisA,mSemiAxisB; // 地球椭球体长轴、短轴半径
double mE2; // 第一偏心率的平方

// 子午线弧长计算系数

double m0 ; // = mSemiAxisA * (1-mE2);
double m2 ; // = 3*mE2*m0/2;
double m4 ; // = 5*mE2*m2/4;
double m6 ; // = 7*mE2*m4/6;
double m8 ; // = 9*mE2*m6/8;
double a0 ; // = m0+m2/2+3*m4/8+5*m6/16+35*m8/128;
double a2 ; // = m2/2+m4/2+15*m6/32+7*m8/16;
double a4 ; // = m4/8+3*m6/16+7*m8/32;
double a6 ; // = m6/32+m8/16;
double a8 ; // = m8/128;

public:
// 子午线弧长正解
double latitude2merdian(double latitude_in);
// 子午线弧长反解
double merdian2latitude(double x_in);
// 设置椭球长短轴
void SetParameterEarth(double AxisA, double AxisB);

};

file *.cpp

CMapProjection::CMapProjection(void)
{
this->mSemiAxisA = 6378137.; //默认为1984椭球,6378245为克拉索夫斯基椭球.;
this->mSemiAxisB = 6356752.31424518; //默认为1984椭球,6356863.0188为克拉索夫斯基椭球;
this->mE2 = (mSemiAxisA * mSemiAxisA - mSemiAxisB * mSemiAxisB)/(mSemiAxisA * mSemiAxisA); // 第一偏心率的平方

// 字午线弧长计算系数

this->m0 = mSemiAxisA * (1-mE2);
this->m2 = 3*mE2*m0/2;
this->m4 = 5*mE2*m2/4;
this->m6 = 7*mE2*m4/6;
this->m8 = 9*mE2*m6/8;
this->a0 = m0+m2/2+3*m4/8+5*m6/16+35*m8/128;
this->a2 = m2/2+m4/2+15*m6/32+7*m8/16;
this->a4 = m4/8+3*m6/16+7*m8/32;
this->a6 = m6/32+m8/16;
this->a8 = m8/128;
}

CMapProjection::~CMapProjection(void)
{
}

// 子午线弧长正解
double CMapProjection::latitude2merdian(double latitude_in)
{
return a0*latitude_in-a2*sin(2*latitude_in)/2+a4*sin(4*latitude_in)/4-a6*sin(6*latitude_in)/6+a8*sin(8*latitude_in)/8;
}

// 子午线弧长反解
double CMapProjection::merdian2latitude(double x_in)
{
double b0 = x_in/a0;
while(1)
{
double F = -a2*sin(2*b0)/2 + a4*sin(4*b0)/4 -a6*sin(6*b0)/6 +a8*sin(8*b0)/8;
double B = (x_in - F)/a0;
if(abs(b0-B)<1.0E-6) return B;
b0 = B;
}
}

// 设置椭球长短轴
void CMapProjection::SetParameterEarth(double AxisA, double AxisB)
{
mSemiAxisA = AxisA;
mSemiAxisB = AxisB;
mE2 = (mSemiAxisA * mSemiAxisA - mSemiAxisB*mSemiAxisB)/(mSemiAxisA * mSemiAxisA);

// 字午线弧长计算系数
this->m0 = mSemiAxisA * (1-mE2);
this->m2 = 3*mE2*m0/2;
this->m4 = 5*mE2*m2/4;
this->m6 = 7*mE2*m4/6;
this->m8 = 9*mE2*m6/8;
this->a0 = m0+m2/2+3*m4/8+5*m6/16+35*m8/128;
this->a2 = m2/2+m4/2+15*m6/32+7*m8/16;
this->a4 = m4/8+3*m6/16+7*m8/32;
this->a6 = m6/32+m8/16;
this->a8 = m8/128;
}