%《MATLAB及其在理工課程中的應用指南》第七章例7－3－1程序
% 阻尼系數對單自由度振動的影響
% 西安電子科技大學出版社出版  陳懷琛編著 1999年10月初版，2004年11月第二版
%
clear,wn=10;tf=2;
for i=1:2
    if i==1 x0=1;v0=0; else x0=0;v0=1;end
    for j=1:10
        zeta(j)=0.1*j;					                    % 對不同的ζ
        wd(j)=wn*sqrt(1-zeta(j)^2); 		                % 求wd
        a=sqrt((wn*x0*zeta(j)+v0)^2+(x0*wd(j))^2)/wd(j); 	% 求振幅A 
        phi=atan2(wd(j)*x0,v0+zeta(j)*wn*x0); 	            % 用atan2是為了求四象限相角
        t=0:tf/1000:tf;				                        % 設定自變量數組
        x(j,:)=a*exp(-zeta(j)*wn*t).*sin(wd(j)*t+phi); 	    % 求過渡過程
    end
    subplot(1,2,i),plot(t,x(1,:), t,x(2,:), t,x(3,:), t,x(4,:), t,x(5,:),...	% 繪圖
        t,x(6,:), t,x(7,:), t,x(8,:), t,x(9,:), t,x(10,:))
    set(gcf,'color','w')                                    % 置圖形背景色為白色
    grid on,pause
end
figure(2),mesh(x)                                           % 畫三維圖
set(gcf,'color','w'),pause,                                 % 置圖形背景色為白色