benhuan.txt

绘制苯环的matlab源文件

18 苯环

主程序文件：bh.m 

function bh
m=1; k=50;
S=m/k*diag(ones(1,6));
P=2*diag(ones(1,6))-diag(ones(1,5),-1)-diag(ones(1,5),+1);
P(1,6)=-1; P(6,1)=-1;
[JM,JBB]=eig(S\P)
JB=abs(sqrt(JBB));
a1=[0.2, 0, 0, 0, 0, 0, 0.1];
a2=[ 0, 0.2, 0, 0, 0, 0, 0.1];
a3=[ 0, 0, 0.2, 0, 0, 0, 0.1]; 
a4=[ 0, 0, 0, 0.2 0, 0, 0.1]; 
a5=[ 0, 0, 0, 0, 0.2 0, 0.1];
a6=[ 0, 0, 0, 0, 0, 0.2 0.1];
phi1=0; phi2=0; phi3=0; phi4=0; phi5=0; phi6=0 ; 
t=0:0.01:4; r=1;

figure
for kk=1:7
axis([-1.5*r 1.5*r -1.5*r 1.5*r]);
axis equal
hold on
plot(r.*cos(0:0.1:2*pi),r.*sin(0:0.1:2*pi),'y')
if kk==7
title('一般模式')
else
ti1='简正模 \cdot\cdot\cdot'; ti2=int2str(kk);
ti=[ti1,ti2];
title(ti);
end
for i=1:6
ss=a1(kk)*JM(i,1)*cos(JB(1,1)*t+phi1)+a2(kk)*JM(i,2)*cos(JB(2,2)*t+phi2)...
+a3(kk)*JM(i,3)*cos(JB(3,3)*t+phi3)+a4(kk)*JM(i,4)*cos(JB(4,4)*t+phi4)...
+a5(kk)*JM(i,5)*cos(JB(5,5)*t+phi5)+a6(kk)*JM(i,6)*cos(JB(6,6)*t+phi6);
x{i}= r.*cos((i)*pi/3-ss./r);
y{i}= r.*sin((i)*pi/3-ss./r);
end

for i=1:5
[xp,yp]=plotstring(x{i}(1),y{i}(1),x{i+1}(1),y{i+1}(1),r);
h{i}=plot(xp,yp,'erasemode','xor','linewidth',1);
hh{i}=plot(x{i}(1),y{i}(1),'erasemode','xor','marker','o',...
'markersize',25,'linewidth',2.5,'color','r');
ii=7-i;
hhh{i}=text(x{i}(1)-0.08,y{i}(1)+0.01,int2str(ii),'fontsize',14,...
'erasemode','xor');
end
[xp6,yp6]=plotstring(x{6}(1),y{6}(1),x{1}(1),y{1}(1),r);
h6=plot(xp6,yp6,'erasemode','xor','linewidth',1);
hh6=plot(x{6}(1),y{6}(1),'erasemode','xor','marker','o',...
'markersize',25,'linewidth',2.5,'color','r');
hhh6=text(x{6}(1)-0.08,y{6}(1)+0.01,'1','fontsize',14,'erasemode','xor');

for j=2:2:401
for i=1:5
[xp,yp]=plotstring(x{i}(j),y{i}(j),x{i+1}(j),y{i+1}(j),r);
set(h{i},'xdata',xp,'ydata',yp);
set(hh{i},'xdata',x{i}(j),'ydata',y{i}(j));
set(hhh{i},'position',[x{i}(j),y{i}(j)]);
end
[xp6,yp6]=plotstring(x{6}(j),y{6}(j),x{1}(j),y{1}(j),r); 
set(h6,'xdata',xp6,'ydata',yp6);
set(hh6,'xdata',x{6}(j),'ydata',y{6}(j));
set(hhh6,'position',[x{6}(j),y{6}(j)]);
drawnow;
end 
cla
end
close(gcf);


function [xp,yp]=plotstring(xa,ya,xb,yb,r)
[xa1,ya1]=cart2pol(xa,ya);
xa1=xa1+0.13;
[xa,ya]=pol2cart(xa1,ya1);

[xb1,yb1]=cart2pol(xb,yb);
xb1=xb1-0.13;
[xb,yb]=pol2cart(xb1,yb1);

a=0.13;n=5;
q2=[];
d=sqrt((xa-xb)^2+(ya-yb)^2);
w=2*pi*n/d;
x=xa:0.02:(xa+d);
y=a.*sin(w.*(x-xa));

if xa>xb
q1=pi+atan((ya-yb)/(xa-xb));
else
q1=atan((ya-yb)/(xa-xb));
end
xd=xa+(x-xa).*cos(q1);
yd=ya+(x-xa).*sin(q1);

for i=1:length(y)
if xd(i)<0
q2a=pi+atan((yd(i))/(xd(i)));
q2=[q2 q2a];
else
q2a=atan((yd(i))/(xd(i)));
q2=[q2 q2a];
end
end
xp=(r+y).*cos(q2);
yp=(r+y).*sin(q2);