椭球面积计算公式

yonjay2

        double ParaE1 = m_de11;//WGS2000_E1;//0.00673949677548;  //WGS-2000椭球偏心率e的平方
        double E1 = ParaE1;
        // ParamA—ParamE
        double ParamA = 1.0 + 0.5 * E1 + 0.375 * E1 * E1 + 0.3125 * E1 * E1 * E1 + 0.2734375 * E1 * E1 * E1 * E1;
        double ParamB = 0.1666666666666666666666666667 * E1 + 0.1875 * E1 * E1 + 0.1875 * E1 * E1 * E1 + 0.1822916666666666666666666667 * E1 * E1 * E1 * E1;
        double ParamC = 0.0375 * E1 * E1 + 0.0625 * E1 * E1 * E1 + 0.078125 * E1 * E1 * E1 * E1;
        double ParamD = 0.0089285714285714285714285714 * E1 * E1 * E1 + 0.01953125 * E1 * E1 * E1 * E1;
        double ParamE = 0.0021701388888888888888888889 * E1 * E1 * E1 * E1;


        // 两点围成梯形面积
        double area=0;
        double lDiference = 0.0;    // 经差
        double bDiference = 0.0;    // 纬差
        double bSum = 0.0;          // 纬度和

        // 纬差||纬度和||经差
        bDiference = (B1 - B)/ 2.0;
        bSum = (B1 + B) / 2.0;
        lDiference = (L1 + L) / 2.0;


        double ItemValue[5];
        double bRadius=m_db;//WGS2000_eb;
        // 按照以下计算顺序：短半径*经差*短半径
        double num1 = 2.0 * bRadius * lDiference * bRadius;
        double num2 = cos(bSum);
        double num3 = sin(bDiference);
        double d1,d2,d3,d4,d5;
        d1=ItemValue[0] = num1 * ParamA * num2 * num3;
        d2=ItemValue[1] = num1 * ParamB * sin(3.0 * bDiference) * cos(3.0 * bSum);
        d3=ItemValue[2] = num1 * ParamC * sin(5.0 * bDiference) * cos(5.0 * bSum);
        d4=ItemValue[3] = num1 * ParamD * sin(7.0 * bDiference) * cos(7.0 * bSum);
        d5=ItemValue[4] = num1 * ParamE * sin(9.0 * bDiference) * cos(9.0 * bSum);
        area = ItemValue[0] - ItemValue[1] + ItemValue[2] - ItemValue[3] + ItemValue[4];