rkf78qs.m

调用RKF78的子过程

function [T,X,DT] = RKF78QS(T0,X0,DT0)
% RKF78QS            调用RKF78的子过程，步长控制采用《人造卫星精密轨道确定》
% T                  积分时间,s
% X                  积分变量:位置和速度_列向量,m m/s
% DT0                积分步长,s
% 
% 编写于2005年06月09日
% 修改于2005年06月09日(No.1):

global NREJ NREJT NSTP MAXREJ DTFAIL DTFIX B BLO BUP REMIN DTINC DTDEC RELERR ABSERR DTOLD
% 调用RKF7(8)------------------------------------------------------------------------------------------
% 判断是否选择定步长积分
DT = DT0;
if DTFIX == 'T' 
    [X,TE] = RKF78(T0,X0,DT);
    T      = T0 + DT;   
else
    % 每推进1个DT,步长尝试的次数低于最大尝试次数
    while NREJ < MAXREJ        
        [X,TE] = RKF78(T0,X0,DT);
        AE   = ABSERR/RELERR;
        XMAG = abs(X) + abs(X0) + AE;
        RTE  = max(TE ./ XMAG);
        RTE  = RTE /(7.5*RELERR);
        % 判断误差是否满足精度要求
        if RTE >= 1
            % 此步长计算失败
            DTFAIL = 'T';
            NREJ   = NREJ + 1;
            NREJT  = NREJT + 1;
            % 此步长过大,减小步长
            PCT   = DTDEC ;
            if RTE < BUP
                PCT   = B/(10.*RTE)^.125;                
            end
            DT    = PCT * DT;
            DTOLD = DT;
        else
            T      = T0 + DT;
            break;
        end  
    end
    if NREJ >= MAXREJ
        X  = X0;
        T  = T0;
        error('本次积分尝试步数过多,积分失败！IFLAG = 7');
    else
        % 下一次积分步长
        PCT = DTINC;
        if RTE > BLO  
            PCT = B/RTE^.125;
        end
        if DTFAIL == 'T'
            PCT = min(PCT,1.0);
        end
        DT    = DT * PCT;
        DTOLD = DT;
    end
end
NSTP = NSTP + 1;