墨托投影正反算源码

// 所有经纬度参数均为弧度
// //Krasovsky
// //this->mSemiAxisA = 6378245.; //6378137.;
// //this->mSemiAxisB = 6356863.0188; //6356752.31424518;;
// //WGS84
// mSemiAxisA = 6378137.; //6378245.; 地球椭球体长轴半径
// mSemiAxisB = 6356752.31424518; //;6356863.0188;地球椭球体短轴半径
//
// mE2 = 1 - mSemiAxisB * mSemiAxisB / mSemiAxisA / mSemiAxisA; // 第一偏心率

/C
C 纬 圈 半 径
C 1990-12-29
C 输入参数: B 纬度数(弧度)
C 输出参数: FUN_LR 纬圈曲率半径 (米)
C 调 用: LR=FUN_LR(B)
C/

double CS57Projection::fun_lr(double lat_in)
{
double cosb, sinb;
cosb = cos(lat_in);
sinb = sin(lat_in);
return mSemiAxisA * cosb / sqrt(1. - mE2 * sinbsinb);
}
// 设置基准纬线，并计算基准纬线的纬圈半径
void CS57Projection::SetBaseLine(double base)
{
mCentralMeridian = 0.0;
mBaseLine = base;
lr0 = fun_lr(mBaseLine);
}
/**************************************************************
C 反 双 曲 正 切 函 数
C 1990-12-15
C***************************************************************/

double fun_atanh(double x_in)
{
if ((x_in > -1.) && (x_in < 1.))
return 0.5*log((1. + x_in) / (1. - x_in));
else
{
// QString str1 = QString::fromLocal8Bit(“fun_atanh反双曲正切函数错误报告：”);
// QString str2 = QString::fromLocal8Bit(“输入参数大于等于1或小于等于－1”);
// QMessageBox::critical(NULL, str1, str2,
// QMessageBox::Ok | QMessageBox::Default,
// QMessageBox::Cancel | QMessageBox::Escape, 0);
// //MessageBox(NULL, “\n error x_in of fun_atanh \n”, “”, MB_OK);
return 0;
}
}

// 等量纬度正解
/C
C 等 量 纬 度 函 数
C 1990-12-15
C 输入参数: LATITUDE 纬度数(弧度)
C 输出参数: FUN_Q 等量纬度 (弧度)
C 调 用: Q=FUN_Q(LATITUDE)
C/
double CS57Projection::fun_q(double lat)
{

double sinb;
sinb = sin(lat);
return fun_atanh(sinb) - mE * fun_atanh(sinb*mE);  // mE mE2的开方

}

/C
C 等 量 纬 度 反 解 函 数
C 1990-12-15
C 输入参数: FUN_Q 等量纬度 (弧度)
C 输出参数: LATITUDE 纬度数(弧度)
C 调 用: Q=FUN_invQ(LATITUDE)
C/

double CS57Projection::fun_invq(double q0)
{
double q1, q10, b0, dq;
q10 = q0;
b0 = asin(tanh(q0));
q1 = fun_q(b0);
dq = q0 - q1;
while (fabs(dq) > EPS10)
{
q10 = q10 + dq;
b0 = asin(tanh(q10));// inverse
q1 = fun_q(b0); // solution
dq = q0 - q1;
};
return(b0);

}

// GEO2MCT
void CS57Projection::getMctXY(double lat_in, double lon_in, double *x, double *y)
{
*x = fun_q(lat_in) * lr0;
*y = (lon_in - mCentralMeridian) * lr0; // mCentralMeridian 经线起算，一般为0
}

// MCT2GEO
void CS57Projection::getMctBL(double x, double y, double *lat_out, double *lon_out)
{
*lat_out = fun_invq(x / lr0);
*lon_out = y / lr0 + mCentralMeridian;
}