DoOlemskoy.m

function [T, Y] = DoOlemskoy(y0,t,n,r0,gamma)
% DoOlemskoy(y0,t,n,r0,gamma) Ðåøåíèå ÑÎÄÓ ìåòîäîì Îëåìñêîãî.
%   DoOlemskoy(C) ðåøåíèå ÑÎÄÓ ìåòîäîì Îëåìñêîãî
%   Âåðñèÿ íà÷àëüíàÿ ïîýòîìó åùå íå âñå àñïåêòû ïðîãðàììû êàê ñëåäóåò
%   ïðîäóìàíû, â ñëåäñòâèè ìîãóò âîçíèêàòü ðàçíîãî ðîäà íåáîëüøèå
%   êàçóñû. Ïðîãðàììà ïðîèçâîäèò ðàñ÷åò óæå çàïðîãðàììèðîâàííîé ñõåìû
%   ñòðóêòóðíûì ìåòîäîì 5 ïîðÿäêà òî÷íîñòè. Âûâîäèò íà ýêðàí òðàåêòîðèþ
%   äâèæåíèÿ ýëåìåíòà øàðîâîãî ñêîïëåíèÿ, çàäàííîãî ïîòåíöèàëîì
%   Øóñòåðà-Ïëàììåðà.
%   Âõîäíûå ïàðàìåòðû: âåêòîð í.ó. y0 = (y1,y1',y2,y2'), âðåìÿ çà êîòîðîå
%   ïðîèâçîäèòñÿ îáñ÷åò (t), êîë-âî øàãîâ (n), âåëè÷èíà ïàðàìåòðà ìîäåëè ïîòåíöèàëà
%   ñêîïëåíèÿ r0 (r0), ïàðàìåòð "ãàììà" (gamma)

    global g_GAMMA g_R0
    g_GAMMA = gamma;
    g_R0 = r0;
    
    dt = t/n;
    %syms Y1 Y1_d Y2 Y2_d
    t_out(1) = 0;
    Y1(1)=y0(1);
    Y1_d(1)=y0(2);
    Y2(1)=y0(3);
    Y2_d(1)=y0(4);
    
%     hWB = waitbar(0, 'Ïîæàëóéñòà æäèòå...');
    
    for i = 1:n
        Y0 = [Y1(i), Y1_d(i), Y2(i), Y2_d(i)];
        [Y1(i+1),Y1_d(i+1),Y2(i+1),Y2_d(i+1)] = GetNextPoint(Y0, dt);
        T(i+1) = i*dt;
        
%         waitbar(i/n, hWB);
        %Y_next = GetNextPoint({Y1(n), Y1_d(n), Y2(n), Y2_d(n)}, t/n);
        %Y1(n+1) = Y_next(1);
        %Y1_d(n+1) = Y_next(2);
        %Y2(n+1) = Y_next(3);
        %Y2_d(n+1) = Y_next(4);
    end
    
%     close(hWB);
    
%     Y_out(1) = Y1;
%     Y_out(2) = Y2;
%     Y_out(3) = Y1_d;
%     Y_out(4) = Y2_d;
    
    Y = [Y1
        Y2
        Y1_d
        Y2_d];
%     set(0, 'CurrentFigure', hParentFigure);
%     plot(Y1,Y2,'-o');

    %----------------------------------------------------------------------
    function [y1, y1_d, y2, y2_d] = GetNextPoint (y0, h) %- âîâðàùàåò ñëåäóþùåå çíà÷åíèå ôóíêöèè
        W1 = getW1(h, y0(1), y0(2), y0(3));
        V1 = getV1(h, y0(1), y0(2), y0(3), y0(4), W1);
        W2 = getW2(h, y0(1), y0(2), y0(3), y0(4), W1, V1);
        V2 = getV2(h, y0(1), y0(2), y0(3), y0(4), W1, V1, W2);
        W3 = getW3(h, y0(1), y0(2), y0(3), y0(4), W1, V1, W2, V2);
        V3 = getV3(h, y0(1), y0(2), y0(3), y0(4), W1, V1, W2, V2, W3);

        y1 = y0(1) + h*y0(2) + (2*V1+V2)/6;
        y2_d = y0(4) + (W1+5*W2+4*W3)/10;
        y2 = y0(3) + h*y0(4) + h*(3*W1+10*W2+2*W3)/30;
        y1_d = y0(2) + (4*V1+5*V2+V3)/10;
    end

    %------------------------------------------------------------------------
    function result = getF1(y1, y2_d, y2)
        result = g_GAMMA*y2_d-(1/(1+((y1^2+y2^2)/(g_R0^2)))^(3/2))*(y1/(g_R0^2))+y1;
    end 

    %------------------------------------------------------------------------
    function result = getF2(y1, y2, y1_d)
        result = g_GAMMA*y1_d(1)-(1/(1+((y1(1)^2+y2(1)^2)/(g_R0^2)))^(3/2))*(y2(1)/(g_R0^2))+y1(1);
    end

    %------------------------------------------------------------------------
    function result = getW1(h, y1, y1_d, y2)
        f2_val = getF2(y1, y2, y1_d);
        result = h * f2_val;
    end  

    %------------------------------------------------------------------------
    function result = getV1(h, y1, y1_d, y2, y2_d, W1)
        f1_val = getF1(y1+(1/6)*h*y1_d, y2_d+(1/6)*W1, y2+(1/6)*h*y2_d);
        result = h * f1_val;
    end  

    %------------------------------------------------------------------------
    function result = getW2(h, y1, y1_d, y2, y2_d, W1, V1)
        f2_val = getF2(y1+(h*y1_d/3)+(h*V1/18), y2+(h*y2_d/3)+(h*W1/18), y1_d+V1/3);
        result = h * f2_val;
    end  

    %------------------------------------------------------------------------
    function result = getV2(h, y1, y1_d, y2, y2_d, W1, V1, W2)
        f1_val = getF1(y1+(2*h*y1/3)+h*V1/4, y2_d-(W1/12)+3*W2/4, y2+(2*h*y2_d/3)+(h*W1/15)+h*W2/8);
        result = h * f1_val;
    end   

    %----------------------------------------------------------------------
    function result = getW3(h, y1, y1_d, y2, y2_d, W1, V1, W2, V2)
        f2_val = getF2(y1+(5*h*y1/6)+(5*h*V1/18)+5*h*V2/72, y2+(5*h*y2_d/6)+(5*h*W1/144)+5*h*W2/16, y1_d+(5*V1/12)+5*V2/12);
        result = h * f2_val;
    end

    %------------------------------------------------------------------------
    function result = getV3(h, y1, y1_d, y2, y2_d, W1, V1, W2, V2, W3)
        f1_val = getF1(y1+h*y1_d+(5*h*V1/36)+5*h*V2/18, y2_d+(3*W1/4)-(5*W2/12)+2*W3/3, y2+h*y2_d-(h*W1/24)+5*h*W2/8);
        result = h * f1_val;
    end
   
end