pendulum_waves_3D

function pendulum_waves_3D
% pendulum waves
clc;clear;close all
load('matlab.mat');
dt=0.01;
t=0:dt:T-dt;
d=0.1; % 各摆间距
opts=odeset('RelTol',1e-10,'AbsTol',1e-10, ...
    'MaxStep',1e-2,'InitialStep',1e-4);
for n=1:length(L)
    [~,u]=ode45(@(t,u)pw(t,u,g,L(n)),t,[theta 0],opts);
    z(:,n)=-L(n)*cos(u(:,1));
    y(:,n)=L(n)*sin(u(:,1));
    x(:,n)=zeros(size(y(:,n)))+n*d+d;
end
% 绘图
D=-3.2;%地面高度
W=1.6;
F=3.5;%臂长
[X,Y,Z]=sphere(20);
r=0.05;XS=r*X;YS=r*Y;ZS=r*Z;
for n=1:length(L)
    for ii=1:length(t)
       [XL(:,:,ii,n), ...
        YL(:,:,ii,n), ...
        ZL(:,:,ii,n) ] = line3d( ...
        [x(ii,n) 0 0],[x(ii,n) y(ii,n) z(ii,n)],0.01,10);
    end
end
%% 绘图
f=figure;f.MenuBar='none';%f.ToolBar='none';
dpi=1.25;scr=[1920 1080];res=[1600 900];
fig=gcf;fig.Color='k';
fig.Position = [(scr-res)/2,res]/dpi;% 居中显示
%% 主视图
sbp1=subplot(1,2,1);
sbp1.Position=[0.05 0.05 0.5 0.9];
plot3([0 0 F],[0 0 0],[D 0 0],...
    '-','Color','k','LineWidth',2)
hold on;rotate3d on
[X,Y,Z]=line3d([0 0 D],[0 0 0],0.03,30);
zhu1=surf(X,Y,Z);zhu1.FaceColor='w';zhu1.EdgeColor='none';
[X,Y,Z]=line3d([F 0 0],[0 0 0],0.03,20);
zhu2=surf(X,Y,Z);zhu2.FaceColor='w';zhu2.EdgeColor='none';
qiu(0,0,0,0.03,20,'w');qiu(F,0,0,0.03,20,'w');
fill3([-.5 F+.5 F+.5 -.5],[-W -W W W],[D D D D],'w')
% 摆球与摆杆
for n=1:length(L)
    sp(n,1)=surf(XS+x(1,n),YS+y(1,n),ZS+z(1,n));
    sp(n,1).FaceColor='r';sp(n,1).EdgeColor='none';
    lines(n,1)=surf(XL(:,:,1,n),YL(:,:,1,n),ZL(:,:,1,n));
    lines(n,1).FaceColor='r';lines(n,1).EdgeColor='none';
end
hold off
g1 = light;g1.Visible = 'on';
g1.Color = 'w';
g1.Position = [5 0 3];
g1.Style = 'local';%点光源
material([0.3 0.7 0.8])% 环境、漫反射、镜面反射强度
axis equal;axis off
ax = gca;
ax.Projection = 'perspective';%透视
ax.CameraViewAngle = 6;%视野范围
ax.CameraPosition = [32 1 1];%相机位置
ax.CameraTarget = [F/2 0 D/2];%目标位置
%% 俯视图
sbp3=subplot(1,2,2);
sbp3.Position=[0.55 0.05 0.47 0.9];
plot3([0 0 F],[0 0 0],[D 0 0],...
    '-','Color','k','LineWidth',2)
hold on;rotate3d on
[X,Y,Z]=line3d([0 0 D],[0 0 0],0.03,20);
zhu1=surf(X,Y,Z);zhu1.FaceColor='w';zhu1.EdgeColor='none';
[X,Y,Z]=line3d([F 0 0],[0 0 0],0.03,20);
zhu2=surf(X,Y,Z);zhu2.FaceColor='w';zhu2.EdgeColor='none';
qiu(0,0,0,0.03,20,'w');qiu(F,0,0,0.03,20,'w');
fill3([-.5 F+.5 F+.5 -.5],[-W -W W W],[D D D D],'w')
% 摆球与摆杆
for n=1:length(L)
    sp(n,2)=surf(XS+x(1,n),YS+y(1,n),ZS+z(1,n));
    sp(n,2).FaceColor='r';sp(n,2).EdgeColor='none';
    lines(n,2)=surf(XL(:,:,1,n),YL(:,:,1,n),ZL(:,:,1,n));
    lines(n,2).FaceColor='r';lines(n,2).EdgeColor='none';
end
hold off
g1 = light;g1.Visible = 'on';
g1.Color = 'w';
g1.Position = [3 -2 3];
g1.Style = 'local';%点光源
material([0.3 0.7 0.8])% 环境、漫反射、镜面反射强度
axis equal;axis off
ax = gca;
ax.Projection = 'perspective';%透视
ax.CameraViewAngle = 6;%视野范围
ax.CameraPosition = [F/2 1 30];%相机位置
ax.CameraTarget = [F/2 0 D/2];%目标位置
ax.CameraUpVector = [-1 0 0];
%% 动画
mv = 1;% 是否录制视频，1录制，其它不录制
if mv == 1
v = VideoWriter('D:\Programs\MATLAB\蛇摆.mp4','MPEG-4');
v.Quality = 100;v.FrameRate = 60;
open(v);
end
tic
for ii=1:length(t)
    for k=1:2
        for n=1:length(L)
            sp(n,k).XData=XS+x(ii,n);
            sp(n,k).YData=YS+y(ii,n);
            sp(n,k).ZData=ZS+z(ii,n);
            lines(n,k).XData=XL(:,:,ii,n);
            lines(n,k).YData=YL(:,:,ii,n);
            lines(n,k).ZData=ZL(:,:,ii,n);
        end
    end
    drawnow
    if mv==1;F=getframe(fig);writeVideo(v,F);end
    clc;fprintf(['当前进度',num2str(100*ii/length(t),'%.1f'),'%%\n']);
end
if mv == 1;close(v);end
fprintf(['绘图用时',num2str(toc),'秒\n']);
%% 函数
function du=pw(~,u,g,L)
du=[u(2);-g/L*sin(u(1))];
function [X,Y,Z]=line3d(A,B,R,N)
n=(A-B)';C=[A;B];
a=cross([1 1 1]',n);
b=cross(a,n);
a=R*a/norm(a);b=R*b/norm(b);
t=linspace(0,2*pi,N);
r=a*sin(t)+b*cos(t);
X=C(:,1)+r(1,:);
Y=C(:,2)+r(2,:);
Z=C(:,3)+r(3,:);
function [X,Y,Z,s]=qiu(x,y,z,r,n,Color)
[X,Y,Z]=sphere(n);
X=r*X;Y=r*Y;Z=r*Z;
s=surf(X+x,Y+y,Z+z);
s.FaceColor=Color;
s.EdgeColor='none';