%Matlab script for ECE100 Lab 1
%here the resonant frequency w_0 = 1 for convenience. It is just a scaling factor
w=0.01:.01:20;                  %defines a vector of frequency samples, starting at w=0.01
%				 and ending at w=20, with steps of 0.01
%
%note: to operate on a vector term by term you need the . operators
%for example to square a scalar a use a^2 but to square each term in
%a vector v use v.^2
%to multiply a vector and a scalar use a*v but to multiply two vectors
%term by term use v.*y
%to invert a vector term by term use 1./v
%the ; after a statement suppresses printing of the result
%when in doubt about use of a function (like plot) do "help plot" at the command line
%
figure				%create a new figure window and give it a number
for z=[0.1 0.3 0.7 1 3 10];     %z is the damping factor
    h=1./(1 + j*2*z*w -w.^2);   %calculates the h(w) complex vector for each z in the list
    h2=abs(h).^2;               %need squared magnitude of complex h(w) vector
    plot(w,h2)                  %plot on linear axes
    xlabel('Frequency \omega/\omega_0')	%characters preceded by \ are in TeX convention
    ylabel('Magnitude |H(\omega)|^2')
    title(['ECE100 Lab 1, ' datestr(now,'dd-mmm-yyyy HH:MM:SS') ', W. A. Coles'])
    hold on
end
grid on
%redo with axes in db vs log w
figure
for z=[0.1 0.3 0.707 1 3 10];
    h=1./(1 + j*2*z*w -w.^2);
    h2=abs(h).^2;
    hdb=10*log10(h2);           %convert squared magnitude to dB    
    semilogx(w,hdb)		%plot it
    k=find(w==1);		%now find the place to annotate the plot at w=1
    text(w(k),hdb(k)-1,['\zeta = ' num2str(z)])
    xlabel('Frequency \omega/\omega_0')
    ylabel('Magnitude |H(\omega)| dB')
    title(['ECE100 Lab 1, ' datestr(now,'dd-mmm-yyyy HH:MM:SS') ' , W. A. Coles'])
    axis tight
    hold on
end
grid on
%redo in full Bode format with db over phase vs log w
w=0.01:.01:20;
figure
for z=[0.1 0.3 0.707 1 3 10];
    h=1./(1 + j*2*z*w -w.^2);
    h2=abs(h).^2;
    hdb=10*log10(h2);
    subplot(2,1,1)
    semilogx(w,hdb)
    k=find(w==1);
    text(w(k),hdb(k)-3,num2str(z))
    hold on
    grid on
    axis tight
    ylabel('Magnitude |H(\omega)| dB')
    title(['ECE100 Lab 1, ' datestr(now,'dd-mmm-yyyy HH:MM:SS') ' , W. A. Coles'])
    ph=phase(h)*180/pi;         %get phase of h(w) and convert to degrees
    subplot(2,1,2)
    semilogx(w,ph)
    k=find(ph < -30,1);
    text(w(k),ph(k),num2str(z))
    hold on
    grid on
    axis tight
    xlabel('Frequency \omega/\omega_0')
    ylabel('Phase (deg)')
end