正常重力计算

计算方法

代码

弧度转换角度

double fun1(double x0[3])
//度分秒转化为弧度函数
{
double a = x0[0] + x0[1] / 60 + x0[2] / 3600;
a = a / 180*3.14159265358979;
return a;
}

计算椭球面上高 h 处的正常重力函数

double fun2(double x0[3], double h)//x0[3]为纬度度分秒，h为椭球面上高h(在椭球面上h=0)
//计算椭球面上高h处的正常重力函数
{
//弧度转化
double x = fun1(x0);
//WGS84参考椭球参数
double a = 6378137, GM = 3.986004418e14, w = 7.292115e-5, f = 1.0 / 298.257223563;
//中间量的定义与计算
double b = a - a*f;
double m = w*w*a*a*b / GM, E = sqrt(a*a - b*b);
double e2 = E / b;
double q0 = 1.0 / 2 * ((1 + 3 * b*b / E / E)*atan2(E,b)- 3 * b / E);
double dq0 = 3 * (1 + b*b / E / E)*(1 - b / E *atan2(E , b)) - 1;
//计算赤道上正常重力
double ra = GM / a / b*(1 - m - m*e2*dq0 / 6.0/ q0);
//计算两极上正常重力
double rb = GM / a / a*(1 + m*e2*dq0 / 3.0 / q0);
//计算椭球面上的正常重力
double r = 100*(a*ra*cos(x)*cos(x) + b*rb*sin(x)*sin(x)) / sqrt(a*a*cos(x)*cos(x) + 
b*b*sin(x)*sin(x));
//计算椭球以上的正常重力值
double rh = r*(1 - 2 * (1 + f + m - 2 * f*sin(x)*sin(x))*h / a + 3 * h*h / a / a);
return rh;
}