与 GMAT（通用任务分析工具）相比，重力扰动状态传播的 C++ 代码显示错误

Question

使用 Runge kutta 7 (8) 方法开发了引力摄动模型和传播模型。该模型使用 EGM96 模型的非标准化系数。与 GMAT 相比，它显示出大约几米的误差。可能是什么原因？

结果有几米的差异。 C++ 是否有精度限制？它会导致错误吗？

   vector<double> GravityPerturbation::gravity(int t, vector<double> y)
{
r2 = y[0] * y[0] + y[1] * y[1] + y[2] * y[2];
    r1 = sqrt(r2);
    r3 = r2 * r1;

    if (lgrav == 0 && mgrav == 0)
    {
        // keplerian motion
    }
    else if (lgrav == 2 && mgrav == 0)
    {
         j2 only

         r5 = r2 * r3;

         d1 = -1.5 * Constants::j2 * Constants::req * Constants::req * Constants::mu / r5;

        d2 = 1 - 5 * y[2] * y[2] / r2;

         agrav[0] = agrav[0] + y[0] * d1 * d2;

         agrav[1] = agrav[1] + y[1] * d1 * d2;

         agrav[2] = agrav[2] + y[2] * d1 * (d2 + 2);

       
      

        for (int i = 0; i < agrav.size(); ++i)
        {
            agrav[i] = -Constants::mu * y[i] / pow(r1, 3);
        }

    }
    else

    {
 
       
// user defined degree and order GravityPerturbation model
        sr2 = y[0] * y[0] + y[1] * y[1];
        sr1 = sqrt(sr2);
        sphi = y[2] / r1;

        

    
    // date and time input

        int dd, mm, yyyy;
        double uthr, utmin, utsec;
        dd = 1, mm = 1, yyyy = 2000;
        uthr = 0.0, utmin = 0.0, utsec = 0.0;

        double jdate0 = julian(dd, mm, yyyy) + uthr / 24.0 + utmin / 1440.0 + utsec / 86400.0;

        double gst0 = gast1(jdate0);

        // Right ascension of greenwich
        double pmt = fmod(gst0 + Constants::omega * t, 2.0 * Constants::PI);
        

        double fun_tan = atan3(y[1], y[0]);

        // east longitude of the spacecraft
        double lamda = realmod(fun_tan - pmt, 2.0 * Constants::PI);

        //! Degree of GravityPerturbation

        lgrav = n = 2;

        //! Order of GravityPerturbation

        mgrav = m = 2;

        double b = asin(sphi);

        im = mgrav;

        if (mgrav < lgrav)
        {
            im = mgrav + 1;
        }

        P_leg = legendre(n, im, sphi);

        Angle = angle(m, lamda, b);

        e1 = e2 = e3 = 0.0;

        d1 = Constants::req / r1;
        d2 = 1.0;

        int il, il1, im1, im2;

        for (il = 0; il < lgrav; il++)
        {
            f1 = f2 = f3 = 0.0;
            il1 = il + 1;

            for (im = 0; im <= il1; im++)
            {
                if (im > mgrav)
                {
                    break;
                }

                im1 = im;
                im2 = im + 1;

                int pos = findCoefIndex(il1, im1);
              
                d3 = ccoeff[pos] * Angle[0][im1] + scoeff[pos] * Angle[1][im1];

              
                f1 = f1 + d3 * P_leg[il1][im1];
                if (im2 <= il1)
                {
                    f2 = f2 + d3 * (P_leg[il1][im2] - Angle[2][im1] * P_leg[il1][im1]);

                    
                }
                if (im2 > il1)
                {
                    f2 = f2 - d3 * Angle[2][im1] * P_leg[il1][im1];
                }
                if (im != 0)
                {
                    f3 = f3 + im * (scoeff[pos] * Angle[0][im1] - ccoeff[pos] * Angle[1][im1]) * P_leg[il1][im1];

                    
                }

                
            }
            

            d2 = d2 * d1;
            e1 = e1 + d2 * (il1 + 1.0) * f1;
            e2 = e2 + d2 * f2;
            e3 = e3 + d2 * f3;
        }

        d3 = Constants::mu / r1;

        e1 = -e1 * d3 / r1;
        e2 = e2 * d3;
        e3 = e3 * d3;

        d1 = e1 / r1;
        d2 = y[2] / (r2 * sr1) * e2;
        d3 = e3 / sr2;

        agrav[0] = agrav[0] + (d1 - d2) * y[0] - d3 * y[1];
        agrav[1] = agrav[1] + (d1 - d2) * y[1] + d3 * y[0];
        agrav[2] = agrav[2] + d1 * y[2] + sr1 * e2 / r2;

       

    }

    for (int i = 0; i < 3; i++)
    {
        agrav[i] = agrav[i] - Constants::mu * y[i] / r3;
    }
   

    return agrav;