%sys_analysis is the function used in CtrlLAB for administrating system analysis tasks.
%
%Available Analysis Tasks are:
%
%   1 for Bode diagrams
%   2 for Nyquist plots
%   3 for Nichols charts
%   4 for inverse Nyquist plots
%   5 for root loci
%   6 for step responses
%   7 for impulse responses
%   8 for stability margins
%   9 for H2 and Hinf norms evaluation
%   10 for time moments 
%   11 for Markov parameters
%   12 for analytical solutions
%
%Available lists of functions under this module
%
%  make_model -- prepares models to be analysed.
%  plt_range_dat -- updates the plot range parameters
%  refine_plot -- refines the plots according to the preset format
%  dly_freq -- evaluates the frequency response of systems with delay
%  time_resp -- evaluates the time domain response of the given system
%  plt_common -- extracts common setting in system analysis
%  make_model -- prepares models for further analysis tasks
%  tim_moments -- evaluates the time moments of the system
%  markov_paras -- evaluates the Markov parameters of the system
%  asol_disp -- displays the analytical solutions of the LTI systems
%  p_frac -- performs trianglar manipulations on the residue() function
%  sim_pars -- 
%  bode_asymp --
%  rloc_asymp --
%

%Designed by Professor Dingyu Xue
%School of Information Science and Engineering, Northeastern University
%Shenyang 110006, P R China
%Email: xue_dy@hotmail.com
%
%This module is for CtrllAB 3.0, (c) 1996-1999
%Last modified 5 October, 1999
%---------------------------------------------------------------------------------

function sys_analysis(nTask,arg1,arg2)
switch nargin
case 0
   %call simulation parameters dialog box 
   sim_pars;
case 1   
   %bring the main window to front
   figure(findobj('Tag','CtrlLABMain')); 
   uu=get(gcf,'UserData'); g1=get(uu{1}(1),'UserData');
   if length(g1)==0, 
      warndlg('Plant model dose not exist, please enter it first!','Matrix Processor Warning:'); 
   else
      uu{1}(7)=nTask; set(gcf,'UserData',uu);
      g_comp=extra_funs(5,4,'Checked',26:27,[]);
      if any([1:7,14:16]==nTask), [x1,x2,n_points,k_extra]=plt_range_dat(nTask); end
      switch nTask
      case {1,2,3,4,14,15}
         if length(x1)==0, [x,y,w,arg4]=dly_freq(k_extra(2),g_comp(1),nTask);
         else,
            w=logspace(log10(x1),log10(x2),n_points);
            [x,y,w0,arg4]=dly_freq(k_extra(2),g_comp(1),nTask,w);
         end
         if length(g_comp)>1, [x1,y1,w0,arg5]=dly_freq(k_extra(2),g_comp(2),nTask,w); end
         switch nTask
         case {1,14}, graf_tool(1);
            if length(g_comp)>1,
               subplot(211), semilogx(w,20*log10(x),w,20*log10(x1)); 
               if nTask==14, line(arg5(1,:),arg5(2,:)); end
               ylabel('Magnitude (dB)'), title('Bode Diagram')
               subplot(212), semilogx(w,y,w,y1); ylabel('Phase (degrees)')
            else   
               subplot(211), semilogx(w,20*log10(x)); 
               if nTask==14, line(arg4(1,:),arg4(2,:)); end
               ylabel('Magnitude (dB)'), title('Bode Diagram')
               subplot(212), semilogx(w,y); ylabel('Phase (degrees)')
            end
            title('Frequenct \omega (Radian)')
         case {2,4,15}, 
            if nTask<=7, graf_tool(nTask); else, graf_tool(2); end
            p_tmp=pi/180*y; H=x.*(cos(p_tmp)+sqrt(-1)*sin(p_tmp)); H=[H,conj(H)];
            if length(g_comp)>1,
               p_tmp=pi/180*y1; H1=x1.*(cos(p_tmp)+sqrt(-1)*sin(p_tmp)); H1=[H1,conj(H1)];
               if any([2,15]==nTask), h=plot(real(H),imag(H),real(H1),imag(H1));
               else, h=plot(real(1./H),imag(1./H),real(1./H1),imag(1./H1)); end   
            else   
               if any([2,15]==nTask), h=plot(real(H),imag(H));
               else, h=plot(real(1./H),imag(1./H)); end   
            end
            set(h,'UserData',w); xlabel('Real Axis'); ylabel('Imaginary Axis'), 
            switch nTask
            case 2, title('Nyquist Plot')
            case 4, title('Inverse Nyquist Plot'); 
            case 15, title('ATAN Nyquist Plot')
            end   
         case 3, graf_tool(3);
            bHold=ishold;
            if ~bHold, delete(extra_funs(2)); ngrid('new'); end
            if length(g_comp)>1, h=plot(y,20*log10(x),y1,20*log10(x1)); 
            else, h=plot(y,20*log10(x)); end
            set(h,'UserData',w); xlabel('Phase (degree)')
            ylabel('Magnitude (dB)'), title('Nichols Chart')
            if ~bHold, hold off; end
         end         
      case {5,16}, 
         %for root locus analysis, only the open-loop model is meaningful.
         [G_forward,G_Sys,Gc_Sys,H_Sys,Td,nPade,GP,g_loop,g_reltn]=plt_common(1,g_comp(1));
         graf_tool(5);
         if length(x1)==0, [r,k]=rlocus(G_forward);
         else, 
            if ~isfinite(x2), x2=1000; end
            if abs(x1)<eps, x1=0.001; end
            k=logspace(log10(x1),log10(x2),n_points); r=rlocus(G_forward,k);
         end
         if length(g_comp)>1,
            r1=rlocus(G_Sys,k); h=plot(real(r'),imag(r'),real(r1'),imag(r1')); 
            [z,p,k]=zpkdata(G_Sys,'v');
            if length(z)>0, h=line(real(z),imag(z),'Marker','o','LineStyle','none'); end
            if length(p)>0, h=line(real(p),imag(p),'Marker','x','LineStyle','none'); end
         else, h=plot(real(r'),imag(r')); end
         [z,p,k]=zpkdata(G_forward,'v');
         if length(z)>0, h=line(real(z),imag(z),'Marker','o','LineStyle','none'); end
         if length(p)>0, h=line(real(p),imag(p),'Marker','x','LineStyle','none'); end
         set(h,'UserData',k); xlabel('Real Axis'); ylabel('Imaginary Axis'), title('Root Locus')
         if nTask==16, rloc_asymp(G_forward); if length(g_comp)>1, rloc_asymp(G_Sys); end, end
      case {6,7}, 
         if length(x1)==0
            [y,t]=time_resp(k_extra,g_comp(1),nTask); 
            if length(y)==0; return; end
         else
            t=x1:(x2-x1)/n_points:x2; y=time_resp(k_extra,g_comp(1),nTask,t);
         end   
         if length(g_comp)>1, y1=time_resp(k_extra,g_comp(2),nTask,t); end   
         graf_tool(nTask);
         if length(g_comp)>1, plot(t,y,t,y1); else, plot(t,y); end
         xlabel('Time (Sec)'); ylabel('Response Magnitude');
         if nTask==6, title('Step Response');
         else, title('Impulse Response'); end
      case 8, %gain phase margins
         figure(findobj('Tag','CtrlLABMain')); uu0=get(gcf,'UserData');
         g1=get(uu0{1}(1),'UserData'); G_Sys=g1{2}; 
         g2=get(uu0{1}(2),'UserData'); if length(g2)>0, Gc_Sys=g2{2}; else, Gc_Sys=1; end
         g3=get(uu0{1}(3),'UserData'); if length(g3)>0, H_Sys=g3{2}; else, H_Sys=1; end
         g4=get(uu0{1}(4),'UserData'); Td=g4{1}; nPade=g4{2};
         if Td>0, [nP,dP]=pade(Td,nPade); GP=tf(nP,dP); else, GP=1; end
         G1=G_Sys*H_Sys*GP; [Gm,Pm,Wcg,Wcp]=margin(G1);
         if length(g_comp)>1, G2=G1*Gc_Sys; [Gm1,Pm1,Wcg1,Wcp1]=margin(G2); end
         display_str; 
         display_str(0.1,0.8,['Gain Margin: ',display_str(Gm),' at \omega=' display_str(Wcg)],[1,0,0]);
         display_str(0.1,0.6,['Phase Margin: ',display_str(Pm) '^o at \omega=' display_str(Wcp)],[1,0,0]);
         if length(g_comp)>1, 
            display_str(0.1,0.4,['Gain Margin: ',display_str(Gm1),' at \omega=' display_str(Wcg1)],[0,0,1]);
            display_str(0.1,0.2,['Phase Margin: ',display_str(Pm1) '^o at \omega=' display_str(Wcp1)],[0,0,1]);
         end   
      case {9,10,11,12}
         [G_Sys1,g1]=make_model(g_comp(1)); 
         if length(g_comp)>1, [G_Sys2,g2]=make_model(g_comp(2)); end
         switch nTask
         case 9, %H_2 & H_inf norm evaluation
            H2=norm(g1,2); Hinf=norm(g1,inf); display_str; 
            display_str(0.1,0.8,['H_2 Norm=',display_str(H2),', H_\infty Norm=' display_str(Hinf)],[1,0,0]);
            if length(g_comp)>1, 
               H2=norm(g2,2); Hinf=norm(g2,inf);
               display_str(0.1,0.5,['H_2 Norm=',display_str(H2),', H_\infty Norm=' display_str(Hinf)],[0,0,1]);
            end
         case 10, % time moments   
            M1=tim_moments(g1,7); display_str;
            display_str(0.1,0.8,'Time Moments:',[1,0,0]); display_str(0.1,0.65,M1,[1,0,0]);
            if length(g_comp)>1, 
               M2=tim_moments(g2,7); 
               display_str(0.1,0.4,'Time Moments:',[0,0,1]); display_str(0.1,0.25,M2,[0,0,1]);
            end
         case 11, %Markov parameters   
            M1=markov_paras(g1,7); display_str;
            display_str(0.1,0.8,'Markov Parameters:',[1,0,0]); display_str(0.1,0.65,M1,[1,0,0]);
            if length(g_comp)>1, 
               M2=markov_paras(g2,7); 
               display_str(0.1,0.4,'Markov Parameters:',[0,0,1]); display_str(0.1,0.25,M2,[0,0,1]);
            end
         case 12, %analytic solutions
            key=get(uu{4}(7),'Checked'); display_str;
            xL(1)=asol_disp(g1,key,0.1,0.8,[1,0,0]); 
            if length(g_comp)>1, xL(2)=asol_disp(g2,key,0.1,0.3,[0,0,1]); end
         end   
      case 13
         [G_Sys1,g1]=make_model(g_comp(1),[]);  
         if length(g_comp)>1, [G_Sys2,g2]=make_model(g_comp(2),[]); end
         inv_T=1/norm(g1,inf); t1=1-inv_T; t2=1+inv_T; t3=-2*asin(0.5*inv_T)*180/pi; display_str;
         display_str(0.1,0.8,['Guaranteed Gain Margin: [' num2str(t1) ',' num2str(t2) ']'],[1,0,0]); 
         display_str(0.1,0.6,['Guaranteed Phase Margin: [' num2str(t3) '^o,' num2str(-t3) '^o]'],[1,0,0]); 
         if length(g_comp)>1, 
            inv_T=1/norm(g2,inf); t1=1-inv_T; t2=1+inv_T; t3=-2*asin(0.5*inv_T)*180/pi;
            display_str(0.1,0.4,['Guaranteed Gain Margin: [' num2str(t1) ',' num2str(t2) ']'],[0,0,1]); 
            display_str(0.1,0.2,['Guaranteed Phase Margin: [' num2str(t3) '^o,' num2str(-t3) '^o]'],[0,0,1]); 
         end
      end      
      if nTask<=7, refine_plot(nTask); end
   end
case 2
   uu=get(gcf,'UserData'); g_graf=gcf;
   switch arg1
   case 1
      kk=get(uu{3}(1),'Value'); set(uu{3}(2),'Value',~kk);
      set(uu{1}(1:6),'Visible',extra_funs(6,~kk)); set(uu{1}(7:8),'Visible',extra_funs(6,kk));
   case 2   
      kk=get(uu{3}(2),'Value'); set(uu{3}(1),'Value',~kk);
      set(uu{1}(1:6),'Visible',extra_funs(6,kk)); set(uu{1}(7:8),'Visible',extra_funs(6,~kk));
   case 4
      keyAlg=get(uu{4},'Value'); ww(6)=num2str(get(uu{1}(9),'String'));
      bFixedStep=get(uu{3}(1),'Value');
      if bFixedStep==0,
         ww(2)=get(uu{1}(2),'String'); eval(['tmin=',ww(2),';']);
         ww(3)=get(uu{1}(4),'String'); eval(['tmax=',ww(3),';']);
         ww(4)=get(uu{1}(6),'String'); eval(['terr=',ww(4),';']);
      else, 
         ww(5)=get(uu{1}(8),'String'); eval(['tmin=',ww(5),';']);
      end
      ww(7)=get(uu{1}(10),'String'); ww(8)=get(uu{1}(11),'String');
      ww(9)=get(uu{1}(12),'String'); ww(1)=get(uu{1}(13),'String');
      uu{6}=ww; set(g_graf,'UserData',uu);
   end   
case 3
   [mag,pha,w]=bode(arg1);
   subplot(211), semilogx(w,20*log10(mag(:)));
   ylabel('Magnitude (dB)'), title('Bode Diagram')
   subplot(212), semilogx(w,pha(:)); ylabel('Phase (degrees)')
   title('Frequenct \omega (Radian)')
end

%-------------------------------------------------------
%plt_range_dat is used to set the plot range parameters.
%-------------------------------------------------------
function [x1,x2,n_points,k_extra]=plt_range_dat(nTask)
%output argument initialization
key=0;
x1=[]; x2=[]; n_points=[]; k_extra=zeros(1,2); k_match=0;
figure(findobj('Tag','CtrlLABMain')); uu=get(gcf,'UserData');
if length(uu)>=10, 
   u_data=uu{10}; if length(u_data)>0, k_match=u_data(5); end
end

switch k_match
case {1,2,3,4}
   %set frequency response data
   if any([1:4]==nTask), key=1; end
case 5,
   %set root locus data
   if nTask==5, key=1; end
case {6,7}   
   %set time-domain response data
   if any([6:7]==nTask), key=1; end
end
if key==1
   x1=u_data(1); x2=u_data(2); n_points=u_data(3); k_extra=u_data([4,6]);
end

%-----------------------------------------------------------------------
%refine_plot is used to refine the plots according to the preset format.
%-----------------------------------------------------------------------
function refine_plot(nTask)
it=get(gca,'Title'); pp=get(it,'Position'); xL=get(gca,'XLim'); yL=get(gca,'YLim');
set(it,'Position',[pp(1), 1.006*yL(2)-0.006*yL(1),pp(3)])
if nTask==1
   ii=extra_funs(2); axes(ii(2)); it=get(gca,'Title'); pp=get(it,'Position');
   yL=get(gca,'YLim'); set(it,'Position',[pp(1), 1.006*yL(2)-0.006*yL(1),pp(3)])
else, 
   ix=get(gca,'xLabel'); x0=0.85*xL(2)+0.15*xL(1);
   y0=0.08*yL(2)+0.92*yL(1); set(ix,'Position',[x0,y0,0]);
end

%update plot properties according to the settings
graf_tool(2,1); graf_tool(2,2); graf_tool(2,4); 

%-------------------------------------------------------------------------------------
%dly_freq is a function used to evaluate the frequency response of a system containing 
%a pure time delay element exp(-tdly*s) over the frequency range of W.       
%
%     [mag,pha,w,arg4]=dly_freq(k_Pade,nCtrl,nTask,w)
%
%where
%  mag,pha -- the magnitude and phase over w
%  w -- the frequency vector
%  nCtrl -- whether the compensated system are required
%  nTask -- the system analysis task code
%  k_Pade -- whether Pade approximation is used for delay
%-------------------------------------------------------------------------------------
function [mag,pha,w,arg4]=dly_freq(k_Pade,nCtrl,nTask,w)

[G_forward,G_Sys,Gc_Sys,H_Sys,Td,nPade,GP,g_loop,g_reltn]=plt_common(k_Pade,nCtrl);
arg4=[];
if nargin==3
   %selecting the frequency vector w in an automatic way
   [num,den]=tfdata(G_Sys,'v'); w=freqint2(num,den,30); 
end
if g_loop==1, 
   %processing open-loop
   if k_Pade==1 & Td>0, G_forward=G_forward*GP; 
   elseif Td>0, G_forward.Td=Td; end
   G_overall=G_forward; [mag,pha]=bode(G_forward,w); mag=mag(:); pha=pha(:);
else   
   %processing closed-loop
   if Td>0 & k_Pade==0
      Gf=G_forward; Gf.Td=Td;
      [mag,pha]=bode(Gf,w); mag=mag(:); pha=pi*pha(:)/180;
      Hf=mag.*exp(-sqrt(-1)*pha);
      if isa(H_Sys,'double'), H=H_Sys; else, H=nyquist(H_Sys,w); end
      switch g_reltn
      case 1, HH=Hf./(1+H.*Hf); 
      case {2,3}, HH=ones(size(Hf))./(1+H.*Hf); if g_reltn==3, HH=1-HH; end
      end
      mag=abs(HH); pha=phase(HH')'*180/pi;   
   else   
      G_forward=G_forward*GP;
      switch g_reltn
      case 1, G_overall=feedback(G_forward,H_Sys*GP); %closed loop I/O
      case {2,3}, G_overall=feedback(1,G_forward*H_Sys*GP); %closed-loop sensitivity & comp sens
         %evaluate complementary sensitivity function