%% 《理论力学(第3版)》（李俊峰，张雄编著）第13页：1.4节极坐标描述法
%% 根据输入的径向运动方程和极角运动方程绘制点的运动轨迹
%% 输入参数
fun_r = @(t)t;      % 径向运动方程（函数句柄）
fun_theta = @(t)cos(t);  % 极角运动方程（函数句柄）
t = 10;    % 运动时间
mode = 1;  % 是否生成视频文件，1，生成，0，不生成
name = 'Section_1_4';  % 视频文件名

%% 绘制

if mode ==  1
    aviobj = VideoWriter(name);  % 新建视频文件
    open(aviobj);  % 打开文件
end

t1 = 0:t/200:t;
t1(end) = t;

% 初始化笛卡尔坐标x、y，极径r，极角theta
x = zeros(1,length(t1));
y = zeros(1,length(t1));
r = zeros(1,length(t1));
theta = zeros(1,length(t1));

% 计算运动轨迹
for ii  =  1:1:length(t1)
    r(ii) = fun_r(t1(ii));
    theta(ii) = fun_theta(t1(ii));
    x(ii) = r(ii)*cos(theta(ii));  % 存储x坐标
    y(ii) = r(ii)*sin(theta(ii));  % 存储y坐标
end

% 计算绘图范围
r_max = 1.1*max(sqrt(x.^2+y.^2));
l1 = max(x);
l2 = max(y);
l = max(l1,l2);

% 绘图
for ii  =  1:1:length(t1)
    plot([-l,l],[0,0],'k','linewidth',2);  % 绘制x轴    
    axis ('equal');
    hold on;
    
    plot([0,0],[-l,l],'k','linewidth',2);  % 绘制y轴   
    set(gca,'xtick',[],'ytick',[],'xcolor','w','ycolor','w')
    
    arrow([l,0],[1,0],0.05*l,2,0,'k');   % 绘制箭头
    arrow([0,l],[0,1],0.05*l,2,0,'k');
    
    % 设定显示范围
    xlim([-r_max,r_max]);
    ylim([-r_max,r_max]);
    
    % 标注坐标轴   
    text(l,-0.1*l,'\itx','FontName','Times New Roman','FontSize',25)
    text(-0.15*l,l,'\ity','FontName','Times New Roman','FontSize',25)
    
    % 绘制标志杆（表示运动点的极径及延长线），M点及运动轨迹
    plot([0,r_max*cos(theta(ii))],[0,r_max*sin(theta(ii))],'k','linewidth',4);  % 绘制标志杆
    scatter(x(ii),y(ii),50,'o','MarkerEdgeColor','k','MarkerFaceColor','k','LineWidth',1);  %绘制M点
    plot(x(1:ii),y(1:ii),'k','linewidth',0.5);  % 绘制M点运动轨迹

    hold off;
    pause(0.01);  % 暂停0.01s，便于观察
    
    % 输出视频
    if mode == 1
        currFrame  =  getframe;
        writeVideo(aviobj,currFrame);
    end
end

if mode == 1
    close(aviobj);   % 关闭视频文件
end


%% 函数定义

% 功能：绘制箭头
% 输入参数
%    xy： 箭头坐标，1行2列向量
%    fx： 箭头方向，表示箭头末端到箭头尖端的方向，1行2列向量
%    sz： 箭头大小，即箭头三角形底边长度
%    sc： 箭头长宽比，（箭头为等腰三角形）箭头三角形腰长与底边长度之比
%    mode： 模式， 模式为0表示箭头底边位于给定的箭头坐标处，模式为1表示箭头尖端位于给定的箭头坐标处
%    color： 箭头颜色， 与MATLAB绘图颜色代号相同，如'k'表示黑色
function arrow(xy,fx,sz,sc,mode,color)
    sc = sqrt(sc^2-1/4);
    theta = atan2(fx(2),fx(1));
    
    xy1 = [0,sz/2];
    xy2 = [0,-sz/2];
    xy3 = [sc*sz,0];
    
    [temp1,temp2] = rotation(xy1(1),xy1(2),theta);
    xy1 = [temp1,temp2];
    
    [temp1,temp2] = rotation(xy2(1),xy2(2),theta);
    xy2 = [temp1,temp2];
    
    [temp1,temp2] = rotation(xy3(1),xy3(2),theta);
    xy3 = [temp1,temp2];
    
    if mode==0   
        xy1=[xy1(1)+xy(1),xy1(2)+xy(2)];
        xy2=[xy2(1)+xy(1),xy2(2)+xy(2)];
        xy3=[xy3(1)+xy(1),xy3(2)+xy(2)];
    else
        xy1=[xy1(1)+xy(1)-xy3(1),xy1(2)+xy(2)-xy3(2)];
        xy2=[xy2(1)+xy(1)-xy3(1),xy2(2)+xy(2)-xy3(2)];
        xy3=[xy(1),xy(2)];    
    end
    
    if nargin == 5
        fill([xy1(1),xy2(1),xy3(1)],[xy1(2),xy2(2),xy3(2)],'k');
    else
        fill([xy1(1),xy2(1),xy3(1)],[xy1(2),xy2(2),xy3(2)],color);
    end
end

% 功能： 将矢量[x,y]逆时针旋转theta角
% 输入参数
%    x,y： 旋转前坐标，1行n列向量
%    theta： 旋转角度
% 输出参数
%    x_new, y_new： 旋转后坐标，1行n列向量
function [x_new,y_new] = rotation(x,y,theta)
    R = [cos(theta),-sin(theta)
         sin(theta), cos(theta)];
    temp = R*[x;y];
    x_new = temp(1,:);
    y_new = temp(2,:);
end