Advertisement
Y-BH

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
  4. load('matlab.mat');
  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. %% 绘图
  31. f=figure;f.MenuBar='none';%f.ToolBar='none';
  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';
Advertisement
RAW Paste Data Copied
Advertisement