% Моделирование стрельбы ядром из пушки
function yadro
    close all
    %hold on
    alpha= pi/4; % угол стрельбы
    v0= 100;    % начальная скорость
    h= 0.00001;   % шаг по времени точек в массиве xy
                % должен быть постоянным, иначе
                % анимация будет неправильной
            % Будем отслежживать момент падения
    options= odeset('Events',@eventf,...
        'OutputFcn',@outf,...
        'MaxStep',0.01);
    [t,xy]= ode45(@yrp,0:h:100,...
        [0 v0*cos(alpha) 0 v0*sin(alpha)],options);
    xy(end,:)
end
function dxy=yrp(t,z)
    r= 0.3; % Радиус ядра
    mu= 0.3*1.2*pi*r^2/2;   % коэффициент трения для шара
    m= 30;  % масса ядра
    x= z(1); vx= z(2); y= z(3); vy= z(4);
    dxy= [vx;-mu*sqrt(vx^2+vy^2)*vx/m;...
        vy;-mu*sqrt(vx^2+vy^2)*vy/m-9.81];
end
function [value,isterminal,...
    direction]= eventf(t,xy)
        % Если ядро упало, прекращаем интегрирование
    value= xy(3);
    isterminal= 1;
    direction= -1;
end
function r=outf(t,y,flag)
    if strcmp(flag,'init')
        a=axes('XLim',[0 900],'YLim',[0 400]);
        a.NextPlot= 'add';
    end
    if isempty(flag)
        plot(y(1),y(3),'b*')
        hold on
        drawnow;
    end
    r= 0;
end