%Fig_3_4.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');
s=char('k','g','r','b','m','c');%plot characters in different colors
% 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);
R=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));
    gamma2(n,k)=(x(n)*gamma(k))^2;
    R(n,k)=gamma2(n,k)/(((w02x2(n))^2)+gamma2(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,R(:,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')];
text(2.25*w0,(1-(i-1)/(2*nplots)),ttlp); %color plot lables to match curves
%grid on
axis([0,wmax,0,1.1*max(R(:,4))]);
xlabel('Frequency, \omega');
ylabel('Response, \it{R}\rm (\omega)');
end

hold off
