%Fig_3_3.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

% 
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
    cindex=mod(i,6) ;           %color index for curves
figure(1); 
plot(x,y3(:,i),s(i),'LineWidth',LW);
axis([0,wmax,1.1*min(y3(:,4)),1.1*max(y4(:,4))]);
hold on
ttlp=['\color{',ss(cindex,:),'}{\gamma}\color{black} = ',...
    num2str(gamma(1,i),'%5.2f')];
%color plot lables to match curves
text(2*w0,(max(y3(:,4))*(1-(nplots-i)/(nplots))),ttlp);   
%grid on
plot(x,y4(:,i),strcat(s(i),':'),'LineWidth',LW);%axis([0,wmax,0,1.1*peak]);
xlabel('Frequency, \omega');
ylabel('Re (solid), Im (dash) (\it{x}\rm_{p0}(\omega))');
end

hold off
