空间直角坐标直接求解大地纬度.docVIP

空间直角坐标直接求解大地纬度 PAGE 4空间直角坐标直接求解大地纬度1 计算公式依次使用下面的公式进行复数运算：2 VC++代码#include math.h#include complex/*******************************************************************\对一个复数 x 开 n 次方\*******************************************************************/std::complexdouble sqrtn(const std::complexdoublex,double n){ double r = _hypot(x.real(),x.imag()); //模 double a = atan2(x.imag(),x.real()); //辐角 n = 1.0 / n; r = pow(r,n); a *= n; return std::complexdouble(r * cos(a),r * sin(a));}/*******************************************************************\空间直角坐标变换为大地坐标a [in] 椭球长半轴f [in] 扁率倒数，不能为零X,Y,Z [in] 空间直角坐标B,L,H [in] 大地纬度，大地经度，大地高。纬度和经度的单位为弧度\*******************************************************************/void XYZtoBLH(double a,double f ,double X,double Y,double Z ,doubleB,doubleL,doubleH){ double b = a * (1.0 - 1.0 / f); double x = _hypot(X,Y); double y = Z; L = atan2(Y,X); f = 2.0 / ((a + b) * (a - b)); double ax = a * x * f; double by = b * y * f; std::complexdouble b4(-ax,by); std::complexdouble d4( ax,by); std::complexdouble P = (b4 * d4 + 4.0) / -3.0; std::complexdouble Q = (d4 + b4) * (d4 - b4) / 2.0; std::complexdouble D = sqrtn(Q * Q - P * P * P,2.0); std::complexdouble M(-0.5,-0.86602540378443864676372317075294); M = M * sqrtn(Q + D,3.0); std::complexdouble y3 = M + P / M; std::complexdouble m = sqrtn(b4 * b4 + 4.0 * y3,2.0); std::complexdouble n = (b4 * y3 - 2.0 * d4) / m; std::complexdouble z1; D = sqrtn((b4 + m) * (b4 + m) - 8.0 * (y3 + n),2.0); z1 = (-(b4 + m) + D) * 0.25; B = atan2(a * z1.imag(),b * z1.real()); ax = a * z1.imag(); by = b * z1.real(); H = (x * by + y * ax - a * b) / _hypot(ax,by);}3 特殊情况如下图所示：图 SEQ 图 \* ARABIC 1 子午平面