# pendulum_waves_3D

Mar 26th, 2022
3,727
Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
1. function pendulum_waves_3D
2. % pendulum waves
3. clc;clear;close all
5. dt=0.01;
6. t=0:dt:T-dt;
7. d=0.1; % 各摆间距
8. opts=odeset('RelTol',1e-10,'AbsTol',1e-10, ...
9.     'MaxStep',1e-2,'InitialStep',1e-4);
10. for n=1:length(L)
11.     [~,u]=ode45(@(t,u)pw(t,u,g,L(n)),t,[theta 0],opts);
12.     z(:,n)=-L(n)*cos(u(:,1));
13.     y(:,n)=L(n)*sin(u(:,1));
14.     x(:,n)=zeros(size(y(:,n)))+n*d+d;
15. end
16. % 绘图
17. D=-3.2;%地面高度
18. W=1.6;
19. F=3.5;%臂长
20. [X,Y,Z]=sphere(20);
21. r=0.05;XS=r*X;YS=r*Y;ZS=r*Z;
22. for n=1:length(L)
23.     for ii=1:length(t)
24.        [XL(:,:,ii,n), ...
25.         YL(:,:,ii,n), ...
26.         ZL(:,:,ii,n) ] = line3d( ...
27.         [x(ii,n) 0 0],[x(ii,n) y(ii,n) z(ii,n)],0.01,10);
28.     end
29. end
30. %% 绘图
32. dpi=1.25;scr=[1920 1080];res=[1600 900];
33. fig=gcf;fig.Color='k';
34. fig.Position = [(scr-res)/2,res]/dpi;% 居中显示
35. %% 主视图
36. sbp1=subplot(1,2,1);
37. sbp1.Position=[0.05 0.05 0.5 0.9];
38. plot3([0 0 F],[0 0 0],[D 0 0],...
39.     '-','Color','k','LineWidth',2)
40. hold on;rotate3d on
41. [X,Y,Z]=line3d([0 0 D],[0 0 0],0.03,30);
42. zhu1=surf(X,Y,Z);zhu1.FaceColor='w';zhu1.EdgeColor='none';
43. [X,Y,Z]=line3d([F 0 0],[0 0 0],0.03,20);
44. zhu2=surf(X,Y,Z);zhu2.FaceColor='w';zhu2.EdgeColor='none';
45. qiu(0,0,0,0.03,20,'w');qiu(F,0,0,0.03,20,'w');
46. fill3([-.5 F+.5 F+.5 -.5],[-W -W W W],[D D D D],'w')
47. % 摆球与摆杆
48. for n=1:length(L)
49.     sp(n,1)=surf(XS+x(1,n),YS+y(1,n),ZS+z(1,n));
50.     sp(n,1).FaceColor='r';sp(n,1).EdgeColor='none';
51.     lines(n,1)=surf(XL(:,:,1,n),YL(:,:,1,n),ZL(:,:,1,n));
52.     lines(n,1).FaceColor='r';lines(n,1).EdgeColor='none';
53. end
54. hold off
55. g1 = light;g1.Visible = 'on';
56. g1.Color = 'w';
57. g1.Position = [5 0 3];
58. g1.Style = 'local';%点光源
59. material([0.3 0.7 0.8])% 环境、漫反射、镜面反射强度
60. axis equal;axis off
61. ax = gca;
62. ax.Projection = 'perspective';%透视
63. ax.CameraViewAngle = 6;%视野范围
64. ax.CameraPosition = [32 1 1];%相机位置
65. ax.CameraTarget = [F/2 0 D/2];%目标位置
66. %% 俯视图
67. sbp3=subplot(1,2,2);
68. sbp3.Position=[0.55 0.05 0.47 0.9];
69. plot3([0 0 F],[0 0 0],[D 0 0],...
70.     '-','Color','k','LineWidth',2)
71. hold on;rotate3d on
72. [X,Y,Z]=line3d([0 0 D],[0 0 0],0.03,20);
73. zhu1=surf(X,Y,Z);zhu1.FaceColor='w';zhu1.EdgeColor='none';
74. [X,Y,Z]=line3d([F 0 0],[0 0 0],0.03,20);
75. zhu2=surf(X,Y,Z);zhu2.FaceColor='w';zhu2.EdgeColor='none';
76. qiu(0,0,0,0.03,20,'w');qiu(F,0,0,0.03,20,'w');
77. fill3([-.5 F+.5 F+.5 -.5],[-W -W W W],[D D D D],'w')
78. % 摆球与摆杆
79. for n=1:length(L)
80.     sp(n,2)=surf(XS+x(1,n),YS+y(1,n),ZS+z(1,n));
81.     sp(n,2).FaceColor='r';sp(n,2).EdgeColor='none';
82.     lines(n,2)=surf(XL(:,:,1,n),YL(:,:,1,n),ZL(:,:,1,n));
83.     lines(n,2).FaceColor='r';lines(n,2).EdgeColor='none';
84. end
85. hold off
86. g1 = light;g1.Visible = 'on';
87. g1.Color = 'w';
88. g1.Position = [3 -2 3];
89. g1.Style = 'local';%点光源
90. material([0.3 0.7 0.8])% 环境、漫反射、镜面反射强度
91. axis equal;axis off
92. ax = gca;
93. ax.Projection = 'perspective';%透视
94. ax.CameraViewAngle = 6;%视野范围
95. ax.CameraPosition = [F/2 1 30];%相机位置
96. ax.CameraTarget = [F/2 0 D/2];%目标位置
97. ax.CameraUpVector = [-1 0 0];
98. %% 动画
99. mv = 1;% 是否录制视频，1录制，其它不录制
100. if mv == 1
101. v = VideoWriter('D:\Programs\MATLAB\蛇摆.mp4','MPEG-4');
102. v.Quality = 100;v.FrameRate = 60;
103. open(v);
104. end
105. tic
106. for ii=1:length(t)
107.     for k=1:2
108.         for n=1:length(L)
109.             sp(n,k).XData=XS+x(ii,n);
110.             sp(n,k).YData=YS+y(ii,n);
111.             sp(n,k).ZData=ZS+z(ii,n);
112.             lines(n,k).XData=XL(:,:,ii,n);
113.             lines(n,k).YData=YL(:,:,ii,n);
114.             lines(n,k).ZData=ZL(:,:,ii,n);
115.         end
116.     end
117.     drawnow
118.     if mv==1;F=getframe(fig);writeVideo(v,F);end
119.     clc;fprintf(['当前进度',num2str(100*ii/length(t),'%.1f'),'%%\n']);
120. end
121. if mv == 1;close(v);end
122. fprintf(['绘图用时',num2str(toc),'秒\n']);
123. %% 函数
124. function du=pw(~,u,g,L)
125. du=[u(2);-g/L*sin(u(1))];
126. function [X,Y,Z]=line3d(A,B,R,N)
127. n=(A-B)';C=[A;B];
128. a=cross([1 1 1]',n);
129. b=cross(a,n);
130. a=R*a/norm(a);b=R*b/norm(b);
131. t=linspace(0,2*pi,N);
132. r=a*sin(t)+b*cos(t);
133. X=C(:,1)+r(1,:);
134. Y=C(:,2)+r(2,:);
135. Z=C(:,3)+r(3,:);
136. function [X,Y,Z,s]=qiu(x,y,z,r,n,Color)
137. [X,Y,Z]=sphere(n);
138. X=r*X;Y=r*Y;Z=r*Z;
139. s=surf(X+x,Y+y,Z+z);
140. s.FaceColor=Color;
141. s.EdgeColor='none';