Причина несоответствия между методами расчета удельного момента количества движения?

Я пытаюсь численно решить стоимость дельта-V перехода между двумя близкими эллиптическими орбитами, которые наклонены друг к другу. Метод, который я использую, по существу вычисляет векторы скорости начальной орбиты в одном узле, конечной орбиты в противоположном узле, а затем вычисляет переходную орбиту по начальному углу траектории полета, начальному радиусу и конечному радиусу.

Одним из ключевых шагов является вычисление конкретного вектора углового момента и вектора эксцентриситета переходной орбиты, чтобы рассчитать косинусную матрицу перифокального направления к инерциальной для переходной орбиты. Однако, когда я вычисляю вектор углового момента h переходной орбиты в инерциальной системе отсчета из векторного произведения векторов положения и скорости в инерциальной системе отсчета, я обнаруживаю значительную ошибку (относительная ошибка составляет -3,9521). e-8) между величиной этого вектора и скалярным удельным угловым моментом, рассчитанным ранее в коде.

Это странно для меня, потому что скалярный угловой момент используется для вычисления вектора скорости. Я смущен тем, где происходит потеря точности.

Я попытался предоставить входные данные с большей точностью, особенно значение mu, которое я использовал, но это совсем не изменило относительную ошибку. Когда я использую тот же метод векторного произведения для вычисления удельного углового момента орбит 1 и 2, ошибка порядка машинной точности.

        mu = 3.98600437823e+14;

        thetaNT = -55.1582940061466; % deg
        eT = 0.022905923178296;
        aT = 7.243582592195826e+06; % m

        r1A = 7.146263097977215e+06; % m
        v1RA = -1.390985544431790e+02; % m/s
        v1ThetaA = 7.494958913236144e+03; % m/s

        eR1 = [0.355828643065080;-0.934551216774375;0];
        eTheta1 = [0.934551216774375;0.355828643065080;0];

        nCpf1 = [0.263190394679355,-0.840751409136755,0.473146789255815;
            0.880932410956014,0.00949753358184791,-0.473146789255815;
            0.393305102275257,0.541338032000730,0.743144825477394];
        nCpf2 = [0.107314578042381,-0.875080710676727,0.471929370924401;
            0.879361618777851,-0.137938482815824,-0.455736896003458;
            0.463903788257849,0.463903788257849,0.754709580222772];


        v1A = sqrt(v1RA^2 + v1ThetaA^2); % Total speed of orbit 1 at A

        hT = sqrt(aT*mu*(1-eT^2)); % Specific angular momentum of transfer orbit

        eRTB = [-cosd(thetaNT);sind(thetaNT+180);0];
        eThetaTB = [-sind(thetaNT+180);-cosd(thetaNT);0];

        % Calculation of radial speed and tangential speed
        vTRA = mu/hT*eT*sind(thetaNT);
        vTThetaA = mu/hT*(1+eT*cosd(thetaNT));

        vTA = sqrt(vTRA^2+vTThetaA^2);

        vTRB = mu/hT*eT*sind(thetaNT+180);
        vTThetaB = mu/hT*(1-eT*cosd(thetaNT));

        % Conversion of radius and speeds into radius and velocity vectors
        % in perifocal frames
        r1APF1 = r1A.*eR1;
        v1APF1 = v1RA.*eR1 + v1ThetaA.*eTheta1;

        vTBPFT = vTRB.*eRTB + vTThetaB.*eThetaTB;

        v2BPF2 = v2RB.*eR2 + v2ThetaB.*eTheta2;

        % Conversion to inertial reference frame
        r1AN = nCpf1*r1APF1;
        v1AN = nCpf1*v1APF1;

        v2BN = nCpf2*v2BPF2;

        rTAN = r1AN;
        vTAN = v1AN.*(vTA/v1A);

        % Calculation of angular momentum and eccentricity vectors in
        % inertial frame
        hTN = cross(rTAN, vTAN);
        eTN = cross(vTAN, hTN)./mu - rTAN./norm(rTAN);
        diffh = (norm(hTN)-hT)/hT
        diffe = (norm(eTN)-eT)/eT

Я бы ожидал, что diffh и diffe будут порядка машинной точности, примерно 2,2e-16, но они намного больше. В частности, diffh = -3,9689e-08 и diffe = 7,5474e-05.

ОБНОВЛЕНИЕ: ошибка появляется где-то в моем расчете радиального вектора и вектора скорости, если это поможет сконцентрироваться на вашем поиске.