Advertisement
Y-BH

spherical_pendulum_3D

Apr 10th, 2022
2,009
Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
  1. % 球面摆运动模型
  2. clc;clear;close all
  3. % 模型参数、初值、数值计算
  4. g=10;L=1.0;
  5. thetad=15;%释放初始角度DEG
  6. theta=thetad*pi/180;%释放初始角度
  7. phid=90;
  8. phi=phid*pi/180;
  9. % 初值
  10. dphi=0.8386*sqrt(g/L)/sin(theta);
  11. dtheta=0.9*sqrt(g/L);
  12. T=120;%最大计算时间
  13. odeopts=odeset('RelTol',1e-10,'AbsTol',1e-10, ...
  14.     'MaxStep',1e-2,'InitialStep',1e-4,...
  15.     'Events',@(t,u) dyn_evns(t,u,L,theta));
  16. [t,u] = ode45(@(t,u) pendum(t,u,g,L),...
  17.     0:0.01:T,[theta phi dtheta dphi],odeopts);
  18. % 质点位置
  19. x=L.*sin(u(:,1)).*cos(u(:,2));
  20. y=L.*sin(u(:,1)).*sin(u(:,2));
  21. z=-L.*cos(u(:,1));tic
  22. xc=0.05*cos(0:0.1:2*pi);yc=0.05*sin(0:0.1:2*pi);
  23. %% 模型绘图
  24. D=-1.2*L;
  25. % 支架
  26. f=figure;f.MenuBar='none';
  27. plot3([0 0 0 0],[L L -L -L],[D 0 0 D],...
  28.     '-','Color','k','LineWidth',2)
  29. hold on
  30. [X1,Y1,Z1]=line3d([0 L D],[0 L 0],0.02,30);
  31. [X2,Y2,Z2]=line3d([0 -L D],[0 -L 0],0.02,30);
  32. [X3,Y3,Z3]=line3d([0 L 0],[0 -L 0],0.02,30);
  33. zhu1=surf(X1,Y1,Z1);zhu1.FaceColor='w';zhu1.EdgeColor='none';
  34. zhu2=surf(X2,Y2,Z2);zhu2.FaceColor='w';zhu2.EdgeColor='none';
  35. heng=surf(X3,Y3,Z3);heng.FaceColor='w';heng.EdgeColor='none';
  36. qiu(0,L,0,0.02,10,'w');
  37. qiu(0,-L,0,0.02,10,'w');
  38. qiu(0,0,0,0.03,30,'w');
  39. % 3D小球
  40. [X,Y,Z]=sphere(50);
  41. r=0.05;X=r*X;Y=r*Y;Z=r*Z;
  42. sp=surf(X,Y,Z);
  43. sp.FaceColor='w';
  44. sp.EdgeColor='none';
  45. % 3D绳子
  46. n1=[x y z];n0=zeros(size(n1));
  47. for ii=1:length(x)
  48.   [ Xl(:,:,ii), ...
  49.     Yl(:,:,ii), ...
  50.     Zl(:,:,ii) ] = line3d(n1(ii,:),n0(ii,:),0.01,20);
  51. end
  52. sl=surf(Xl(:,:,1),Yl(:,:,1),Zl(:,:,1));
  53. sl.FaceColor='w';
  54. sl.EdgeColor='none';
  55. sl.AmbientStrength=0.3;
  56. % 固定点
  57. plot3(0,0,0,'.','Color','#202020','markersize',20)
  58. a=animatedline('LineStyle','-','LineWidth',2,'Color','#AAAAAA');
  59. tracks=animatedline('Color','#D95319','LineWidth',2);
  60. fill3(1.5*[L -L -L L],1.5*[L L -L -L],[D D D D],'w')
  61. p=patch(xc,yc,zeros(size(xc))+D+0.001,10*16/255*[1 1 1]);
  62. p.EdgeColor='none';
  63. h = light;
  64. h.Style='local';
  65. axis equal;axis off;axis([-L L -L L D 0])
  66. ax=gca;ax.Clipping='off';
  67. ax.Projection = 'perspective';%透视
  68. ax.CameraViewAngle = 10;%视野范围
  69. ax.CameraTarget = [0 0 -0.7];%目标位置
  70. dpi=1.25;scr=[1920 1080];res=[1280 960];
  71. fig=gcf;fig.Color='k';
  72. fig.Position = [(scr-res)/2,res]/dpi;% 居中显示
  73. rotate3d on
  74. %% 动画
  75. mv = 1;% 是否录制视频,1录制,其它不录制
  76. if mv == 1
  77. v = VideoWriter('D:\Programs\MATLAB\球面摆3D.mp4','MPEG-4');
  78. v.Quality = 100;v.FrameRate = 60;
  79. open(v);
  80. end
  81. for n=1:3000
  82.     sp.XData=X+x(n);
  83.     sp.YData=Y+y(n);
  84.     sp.ZData=Z+z(n);
  85.     sl.XData=Xl(:,:,n);
  86.     sl.YData=Yl(:,:,n);
  87.     sl.ZData=Zl(:,:,n);
  88.     p.XData=xc+x(n);p.YData=yc+y(n);
  89.     addpoints(tracks,x(n),y(n),z(n));
  90.     addpoints(a,x(n),y(n),D);
  91.     ax.CameraPosition = [10*sin(t(n)/20) 10*cos(t(n)/20) 6];%相机位置
  92.     drawnow% limitrate
  93.     if mv==1;F=getframe(fig);writeVideo(v,F);end
  94. end
  95. for n=3001:3:length(t)
  96.     sp.XData=X+x(n);
  97.     sp.YData=Y+y(n);
  98.     sp.ZData=Z+z(n);
  99.     sl.XData=Xl(:,:,n);
  100.     sl.YData=Yl(:,:,n);
  101.     sl.ZData=Zl(:,:,n);
  102.     p.XData=xc+x(n);p.YData=yc+y(n);
  103.     addpoints(tracks,x(n),y(n),z(n));
  104.     addpoints(a,x(n),y(n),D);
  105.     ax.CameraPosition = [10*sin(t(n)/20) 10*cos(t(n)/20) 6];%相机位置
  106.     drawnow% limitrate
  107.     if mv==1;F=getframe(fig);writeVideo(v,F);end
  108. end
  109. if mv == 1;close(v);end
  110. fprintf(strcat('绘图保存用时',num2str(toc,4),'秒\n'));
  111. %% 函数
  112. function du = pendum(~,u,g,L)
  113. theta=u(1);phi=u(2);dtheta=u(3);dphi=u(4);
  114. ddtheta=-g/L*sin(theta)+dphi^2*sin(theta)*cos(theta);
  115. ddphi=-2*dtheta*dphi*cot(theta);
  116. du=[dtheta dphi ddtheta ddphi]';
  117. end
  118. function [position,isterminal,direction] = dyn_evns(~,y,L,theta)
  119. position = y(1)/L-sin(theta)+theta*1e-5;
  120. isterminal = 0;
  121. direction = 1;
  122. end
  123. function [X,Y,Z]=line3d(A,B,R,N)
  124. n=(A-B)';C=[A;B];
  125. a=cross([1 1 1]',n);
  126. b=cross(a,n);
  127. a=R*a/norm(a);b=R*b/norm(b);
  128. t=linspace(0,2*pi,N);
  129. r=a*sin(t)+b*cos(t);
  130. X=C(:,1)+r(1,:);
  131. Y=C(:,2)+r(2,:);
  132. Z=C(:,3)+r(3,:);
  133. end
  134. function [X,Y,Z,s]=qiu(x,y,z,r,n,Color)
  135. [X,Y,Z]=sphere(n);
  136. X=r*X;Y=r*Y;Z=r*Z;
  137. s=surf(X+x,Y+y,Z+z);
  138. s.FaceColor=Color;
  139. s.EdgeColor='none';
  140. end
Advertisement
RAW Paste Data Copied
Advertisement