half_cycle

function half_cycle
% 质点离开竖直光滑半圆轨道的轨迹
% 模型参数
R=1;L=4*R;g=10; % 轨道半径、水平轨长、重力加速度
vu=sqrt(g*R);tu=sqrt(R/g); % 单位速度、单位时间
T=100*tu;%模拟时间上限
%半圆轨道
tht=-0.5*pi:0.01:0.5*pi;
xC=R*cos(tht);yC=R+R*sin(tht);
%圆轨道
tht=-pi:0.01:pi;
xO=R*cos(tht);yO=R+R*sin(tht);
% ode设置
opts=odeset('RelTol',1e-6,'AbsTol',1e-6, ...
    'MaxStep',tu*1e-2,'InitialStep',tu*1e-2, ...
    'Events',@(t,u) luodi(t,u,R));
v1=sqrt(2)*vu;%恰好到达圆轨道半高的速度
v2=sqrt(5)*vu;%恰好经过圆轨道顶点的速度
v3=sqrt(7/2)*vu;%恰好落在圆轨道底部的速度
% 恰好通过最高点轨迹
[~,yd] = ode45(@(t,u) half_c(t,u,g,R), ...
    [0 T],[0 0 -0.5*pi v2 v2 0]',opts);
% 恰好落至圆轨道底部轨迹
[~,yb] = ode45(@(t,u) half_c(t,u,g,R), ...
    [0 T],[0 0 -0.5*pi v3 v3 0]',opts);
NN=16;%小球个数
% 微分方程数值计算
for ii=1:NN
    v0=sqrt(ii/2)*vu;%初速度
    [tN{ii},uN{ii}]=ode45(@(t,u)half_c(t,u,g,R), ...
        0:0.005:T,[0 0 -0.5*pi v0 v0 0]',opts);
    LT(ii)=length(tN{ii});
end
% 数组填充对齐
for ii=1:NN
    dL=max(LT)-LT(ii);
    if dL>0
        u=uN{ii};
        uN{ii}=[u;repmat(u(end,:),dL,1)];
    end
    u=uN{ii};
    X(:,ii)=u(:,1);Y(:,ii)=u(:,2);
end
%% 模型绘图
pos1 = [0 0 1 1];pos2 = [0.05 0.05 0.9 0.9];
subplot('Position',pos2)
dpi=1.25;scr=[1920 1080];res=[1280 720];
fig=gcf;fig.Color='w';
fig.Position = [(scr-res)/2,res]/dpi;% 居中显示
%半圆轨道
plot(xC,yC,'Color','k','LineWidth',4);hold on
% 圆轨道
plot(xO,yO,'-.','Color','k','LineWidth',0.5)
%水平轨道
line([-L 0],[0 0],'linestyle','-','Color','k','LineWidth',4)
% 恰好通过最高点轨迹
pd=plot(yd(:,1),yd(:,2),'--','Color','#444444','LineWidth',1.5);
% 恰好落至圆轨道底部轨迹
pb=plot(yb(:,1),yb(:,2),'-','Color','#444444','LineWidth',1);
% 动画线条
for k=NN:-1:1
    p(k)=animatedline('Color','#D95319','LineWidth',1.5);
end
% 小球
pp=plot(0,0,'.','markersize',30,'Color','#D95319');
% 图例
legend([pd pb],{'恰好通过最高点轨迹','恰好落至圆轨道底部轨迹'},...
    'FontSize',16,'FontName','黑体','Location','northwest');
legend('boxoff')
axis([-4 1 0 2]);axis equal;axis tight% axis off
%% 动画部分
mv = 0;% 是否录制视频，1录制，其它不录制
if mv == 1
v = VideoWriter('D:\Programs\MATLAB\小球脱轨.mp4','MPEG-4');
v.Quality = 10;v.FrameRate = 30;
open(v);
end
for ii=1:max(LT)
    for jj=NN:-1:1
        u=uN{jj};
        addpoints(p(jj),u(ii,1),u(ii,2));
    end
    pp.XData=X(ii,:);pp.YData=Y(ii,:);
    drawnow% limitrate
    if mv == 1;F = getframe(fig);writeVideo(v,F);end
end
if mv == 1;close(v);end
%% 微分方程
function du = half_c(~,u,g,R)
x=u(1);y=u(2);t=u(3);v=u(4);vx=u(5);vy=u(6);
a=v^2/R;
if y==0 && x<0
    dx=v;dy=0;dtheta=0;
    dv=0;dvx=0;dvy=0;
elseif x>=0 && g*sin(t)<a
    dx=-v*sin(t);dy=v*cos(t);
    dtheta=v/R;dv=-g*cos(t);
    dvx=-a*cos(t)-dv*sin(t);
    dvy=-a*sin(t)+dv*cos(t);
else
    dx=vx;dy=vy;dtheta=0;
    dv=0;dvx=0;dvy=-g;
end
du=[dx dy dtheta dv dvx dvy]';% 应返回列向量
%% 事件函数
function [position,isterminal,direction] = luodi(~,u,R)
if u(1)>0 && u(2)<R && u(3)>0
    position = u(2)+sqrt(R^2-u(1)^2)-R;% 落在圆轨道
else
    position = u(2);% 落回水平面
end
isterminal = 1;
direction = -1;