# 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;
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. % 支架
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);
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);
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