%Lab 3 part 1 a-i
clear all;
b=1;
a=[1,1];
t=[0:.01:15]';
a_all=[.25, .5, 1, 1.5, 2];
x=ones(size(t));

for i=1:5
    clf
    a(1)=a_all(i);
    hold on  
    plot(t,x);
    xlabel('Time in seconds');
    ylabel('Amplitude');
    title(['Amplitude of Step and Response for a_{0}=', num2str(a_all(i))]);
    y=lsim(b,a,x,t); % output
    plot(t,y,'r'), grid on 
    axis([0 15 0 1.5])
    legend('Step', 'Response')
    line(t,0.99);
    ind=find(y<=0.99);              % ind is composed of points from y that are <= .99 
    t_trans(i)=t(max(ind));          % the maximimum value in ind is where y(t)<= .99    
    disp('Please hit a key')
    pause,
end

    t_trans
    clf
    hold off
    axis([0 10 0 2])
    plot (a_all, t_trans), grid on
    disp('Please hit a key')
    xlabel('a_{0}');
    ylabel('t_{trans}'); 
    title('a_{0} vs. t_{trans}')
    pause,
    
 %Lab 3 part 1 b-i,j
 clear all;
 M=1.5;
 a0=1/(20*pi);
 fvector=([2,4,6,8,10,12]);
 t=[0:.001:15]';
 a=[a0,1];
 b=1;

for i= 1:6
     period=1/fvector(:,i);
     f=fvector(:,i);
     ind=find(t>=5 & t<=(5+period));
     u=ones(size(t));           %We do not need since we are graphing positive time.  
     x=[M*cos(2*pi*f*t)];
     y=lsim(b,a,x,t)
     t_ss=t(ind);
     x_ss=x(ind);
     y_ss=y(ind);
     hold off
     clf
     legend('off');
     hold on
     axis([5 (5+period) -2 2])               % set axis to display 1 oscillation at steady state
     plot(t_ss,[y_ss,x_ss]);
     grid on
     max_x=max(x(ind)) 
     max_y=max(y(ind)) 
     xlabel('Time in seconds');
     ylabel('Amplitude');
     legend ('y (output)', 'x (input)')    
     title('Click on rising edges of waveforms (input first).');
     [xp,yp]=ginput(2);                      % get mouse input on graph to determine phase shift by inspection
     phase = xp(1) - xp(2);
     phase = (phase/period) * 2 * pi
     title(['Steady State Response for f=',num2str(f),' Hz, phase shift=',num2str(phase), 'r, X Max=',num2str(max_x),', Y Max=',num2str(max_y),])
     disp('Please hit a key')
     freq_resp(f/2) = max_y/max_x;            % store B/M in freq_resp vector for future graphing (y value)
     f_storage(f/2)=f ;                       % store f in f_storage vector for future graphing (x value)
     theta2(f/2)=phase;                       % store phase shift in theta vector for future graphing (x value)
     pause, 
 end

  %Lab 3 part 1 b-k
  hold off
  clf
  legend('off');
  hold on
  xlabel('Frequency (Hz)');
  ylabel('B/M');  
  legend ('B/M', 'Frequency')
  title('Magnitude Response (Low-pass)');
  plot (f_storage, freq_resp), grid on
  pause,
  hold off
  clf
  legend('off');
  hold on
  xlabel('Frequency (Hz)');
  ylabel('\theta (Radians)');  
  legend ('\theta', 'frequency');
  title('Phase Response')
  plot (f_storage, theta2), grid on
  pause
  
   
  %Lab 3 part 2-ai
  clear all;
  b=[0,0,2];
  a=[1,sqrt(2),1];
  t=[0:.01:30]';
  x=ones(size(t)); % DC input 
  clf
  plot(t,x);
  hold on
  xlabel('Time in seconds');
  ylabel('Amplitude');
  y=lsim(b,a,x,t); % output
  plot(t,y,'r'), grid on
  legend('Step (input)', 'Response (output)');
  max_y=max(y) % find peak value of y
  line(t,2.01)
  line(t,1.99)
  title('Zoom in on intersection where y enters boundary and press a key.');
  pause,
  title('Click on intersection.')
  [xp,yp]=ginput(1);
  title(['[2.a.i] Response of System Has a PeakPeak value of y:', num2str(max_y),'   y enters steady state at t=', num2str(xp(1))])
  pause
  
  %Lab 3 part 2-aii
  clear all;
  b=[0,0,2];
  a=[1,.5,1];
  t=[0:.01:30]';
  x=ones(size(t)); % DC input 
  clf
  plot(t,x);
  hold on
  xlabel('Time in seconds');
  ylabel('Amplitude');
  y=lsim(b,a,x,t); % output
  plot(t,y,'r'), grid on
  legend('Step (input)', 'Response (output)');
  max_y=max(y) % find peak value of y
  line(t,2.01)
  line(t,1.99)
  title('Zoom in on intersection where y enters boundary and press a key.');
  pause,
  title('Click on intersection.')
  [xp,yp]=ginput(1);
  title(['[2.a.ii] Response of System Has a PeakPeak value of y:', num2str(max_y),'   y enters steady state at t=', num2str(xp(1))])
  pause
  
  %Lab 3 part 2-bi
  clear all;
  b=[0,0,(20*pi)^2];
  t=[0:.01:30]';
  a=[1,(20*pi)/.707,(20*pi)^2];
  f=[0;.01;50]';
  w=2*pi*f;
  H=freqs(b,a,w);
  x=ones(size(t));
  y=lsim(b,a,x,t);
  t_ss=t(ind);
  x_ss=x(ind);
  y_ss=y(ind);
  hold off
  clf
  legend('off');
  hold on
  axis([5 (5+period) -2 2])               % set axis to display 1 oscillation at steady state
  plot(t_ss,[y_ss,x_ss]);
  grid on
  max_x=max(x(ind)) 
  max_y=max(y(ind)) 
  xlabel('Time in seconds');
  ylabel('Amplitude');
  legend ('y (output)', 'x (input)')    
  title('Click on rising edges of waveforms (input first).');
  [xp,yp]=ginput(2);                      % get mouse input on graph to determine phase shift by inspection
  phase = xp(1) - xp(2);
  phase = (phase/period) * 2 * pi
  title(['Steady State Response for f=',num2str(f),' Hz, phase shift=',num2str(phase), 'r, X Max=',num2str(max_x),', Y Max=',num2str(max_y),])
  disp('Please hit a key')
  freq_resp(f/2) = max_y/max_x;            % store B/M in freq_resp vector for future graphing (y value)
  f_storage(f/2)=f ;                       % store f in f_storage vector for future graphing (x value)
  theta2(f/2)=phase;                       % store phase shift in theta vector for future graphing (x value)
  pause,  
  
  

  