%Fig_3_2.m %response of a damped harmonic oscillator with resonance at w0 % clear all; clf; FS = 18; %label fontsize 18 FSN = 16; %number fontsize 16 LW = 2; %linewidth % Change default axes fonts. set(0,'DefaultAxesFontName', 'Times New Roman'); set(0,'DefaultAxesFontSize', FSN); % Change default text fonts. set(0,'DefaultTextFontname', 'Times New Roman'); set(0,'DefaultTextFontSize', FSN); %plot curves in different colors ss=char('black','green','red','blue','magenta','lightBlue'); %plot characters in different colors s=char('k','g','r','b','m','c'); % red, green, yellow, magenta, blue, black, white % gray, darkGreen, orange, and lightBlue %*********************************************************************** eye=complex(0,1); %square-root of minus one delta=pi*1e-6; %small addition to avoid divide by zero npoints=300; %number of points in plot w0=1; %center frequency wmax=3*w0; %maximum frequency gamma=[2 1 0.5 0.25]; %broadening term in frequency nplots=length(gamma); %number of plots peak=1/(w0*gamma(nplots)) %peak w02=w0^2; dx=wmax/npoints; %frequency step x=[0:(dx+delta):wmax]; %frequency [start:increment:stop] w02x2=(w02-x.^2); gammax=zeros(npoints,nplots); gamma2=zeros(npoints,nplots); y=zeros(npoints,nplots); for n=1:1:npoints; for k=1:1:nplots; gammax(n,k)=eye*x(n)*gamma(k); y(n,k)=1./((w02x2(n))-gammax(n,k)); end end % y=1./((w02x2)+gammax); %complex amplitude normalized to unit force and unit mass % y1=abs(y); %absolute value of amplitude y3=real(y); %real amplitude y4=imag(y); %imaginary amplitude y2=atan2(y4,y3); %phase calculated using four quadrants %*********************************************************************** % plot figures %*********************************************************************** for i=1:1:nplots figure(1); plot(x,y1(:,i),s(i),'LineWidth',LW); hold on cindex=mod(i,6) ; %color index for curves ttlp=['\color{',ss(cindex,:),'}{\gamma}\color{black} = ',... num2str(gamma(1,i),'%5.2f')]; %color plot lables to match curves text(2*w0,(0.9*peak*(1.0-(nplots-i)/(2*nplots))),ttlp); axis([0,wmax,0,1.1*peak]); %grid on xlabel('Frequency, \omega'); ylabel('Amplitude (| \it{x}\rm_{p0}(\omega)|)'); figure(2); plot(x,y2(:,i)/pi,s(i),'LineWidth',LW); hold on axis([0,wmax,0,1]); %grid on xlabel('Frequency, \omega'); ylabel('Phase, \phi(\omega/\pi) (radians)'); end hold off