%
%   LTIDEMO2.M  MATLAB sample that demonstrates various analysis schemes for
%          Linear Time Invariant (LTI) systems.  DEMO2 in this series focuses
%          on frequency domain solutions and plotting techniques. 
%          See LTIDEMO1.M for the time domain simulations.
%  
%   Two areas of focus:
%     A. Step response using transfer function form (for comparison with DEMO1)
%     B. Generation and Interpretation of Bode plots
%
%   Simple RLC electrical circuit with applied voltage
%     LI'' + RI' + I/C = E'
%
%   This can be put into standard Laplace transform form as follows:
%      I(s) = G(s)E(s)    with      G(s) = s/(Ls^2 + Rs + 1/C)
%
%   File prepared by J. R. White, UMass-Lowell (last update: Feb. 2020)
%

      clear all,   close all,   nfig = 0;
%
%   define coefficients for specific problem
      L = 0.1;                % inductance (henry)
      Ca = 0.001;             % capacitance (farad)
      RR = [100 10 2];        % resistance (ohms)
      ir = menu('Choose R value in RLC circuit', ...
                'R = 100 ohms (over damped response)', ...
                'R = 10 ohms (under damped response)', ...
                'R = 2 ohms  (slightly damped response)');
      R = RR(ir);
%
%   create the transfer function polynomials and LTI object for this system
      us = 10;   num = [1 0];   den = [L R 1/Ca];   sys = tf(num,den);
%
%   Part A  Step Response of LTI System (to compare to DEMO1)
%
%   set integration time and time vector
      to = 0;   tf = 0.25;   nt = 251;   t = linspace(to,tf,nt); 
%
%   simulate and plot the time domain response for a step input
      nfig = nfig+1;   figure(nfig)
      y = step(sys,t);  y = us*y;   % normalize output by input strength (us volts)
      plot(t,1000*y,'b','LineWidth',2),grid
      title(['LTIDemo2:  Step Response for Series RLC Circuit (R = ', ...
           num2str(R),' ohms)'])
      ylabel('Current (ma)'), xlabel('Time (sec)')
%
%   Part B  Interpretation of Frequency Response Plots (Bode)
%
%   generate bode plot
      nfig = nfig+1;   figure(nfig),   bode(sys,'r'),grid
      h = findobj(gcf,'Color','red'); nh = length(h); 
      for i = 1:nh; set(h(i),'LineWidth',2); end
      title(['LTIDemo2:  Bode Plots for Series RLC Circuit (R = ', ...
           num2str(R),' ohms)'])
%
%   now simulate and plot system response for two different frequencies
      t1 = linspace(0,20,1001);   w1 = 1;
      u1 = sin(w1*t1)';           y1 = lsim(sys,u1,t1);
      n1 = polyval(num,w1*1i);     d1 = polyval(den,w1*1i);     f1 = abs(n1/d1);
      t2 = linspace(0,0.2,1001);  w2 = 100;
      u2 = sin(w2*t2)';           y2 = lsim(sys,u2,t2);
      n2 = polyval(num,w2*1i);     d2 = polyval(den,w2*1i);     f2 = abs(n2/d2);
%
      nfig = nfig+1;   figure(nfig)
      subplot(2,1,1),plot(t1,u1,'r:',t1,y1/f1,'b-','LineWidth',2),grid
      title(['LTIDemo2:  Input-Output Relationship for \omega = ',num2str(w1), ...
             ' rad/sec (R = ',num2str(R),' ohms)'])
      v1 = axis;   v1(3) = -1.1;   v1(4) = 1.1;   axis(v1);
      xlabel('time (sec)')
      gtext('E(t)'),gtext('I(t)'),gtext(['ratio = ',num2str(f1,'%5.4f')])
      subplot(2,1,2),plot(t2,u2,'r:',t2,y2/f2,'b-','LineWidth',2),grid
      title(['LTIDemo2:  Input-Output Relationship for \omega = ',num2str(w2), ...
             ' rad/sec (R = ',num2str(R),' ohms)'])
      v2 = axis;   v2(3) = -1.1;   v2(4) = 1.1;   axis(v2);
      xlabel('time (sec)')
      gtext('E(t)'),gtext('I(t)'),gtext(['ratio = ',num2str(f2,'%5.4f')])
%
%  end of demo
