function busFUN(action)

%**************************************************************%
%This functions is intended for use with the busGUI.m file.    %
%It contains the code necessary to run the GUI's animation     %
%and controls.                                                 %
%                                                              %
%Copyright (C) 1997 by the Regents of the University of        %
%Michigan.                                                     %
%**************************************************************%

%Callback for the RUN button%
if action ==1

  m1=2500;
  m2=320;
  k1 = 80000;
  k2 = 500000;
  b1 = 350;
  b2 = 15020;

  nump=[(m1+m2) b2 k2];
  denp=[(m1*m2) (m1*(b1+b2))+(m2*b1) (m1*(k1+k2))+(m2*k1)+(b1*b2) (b1*k2)+(b2*k1) k1*k2];

  num1=[-(m1*b2) -(m1*k2) 0 0];
  den1=[(m1*m2) (m1*(b1+b2))+(m2*b1) (m1*(k1+k2))+(m2*k1)+(b1*b2) (b1*k2)+(b2*k1) k1*k2];

  numf=num1;
  denf=nump;
      
  z1=3+3.5i;
  z2=3-3.5i;
  p1=30;
  p2=60;
  numc=conv([1 z1],[1 z2]);
  denc=conv([1 p1],[1 p2]);

%Get the value of the gain from the slider%
  gainhandle = findobj('Tag','gainslider');
  k=get(gainhandle,'Value');

  numc=k*numc;
  numa=conv(conv(numf,nump),denc);
  
  poly1=conv(denp,denc);
  poly2=conv(nump,numc);

%Polyadd function%  
  if length(poly1)<length(poly2)
    short=poly1;
    long=poly2;
  else
    short=poly2;
    long=poly1;
  end
  mz=length(long)-length(short);
  if mz>0
    poly=[zeros(1,mz),short]+long;
  else
    poly=long+short;
  end
  
  dena=conv(denf,poly); 

%Run the simulation with a 0.1m step disturbance%     
  t=0:0.04:7;
  u=0.1*ones(size(t));
  [Y]=lsim(numa,dena,u,t);
  
  lt=length(t);
  changet=t;
  for nt=1:lt-1,
    changet(1:23)=-0.91:0.04:0;
    changet(nt+23)=t(nt);
  end
  t=changet;  
  
  lY=length(Y);
  change=Y;
  for nc=1:lt-1,
    change(1:23)=change(1:23)*0;
    change(nc+23)=Y(nc);
  end  
  Y=change;

%Check if response should be plotted separately%  
  phandle = findobj('Tag','plotbox');
  plotval=get(phandle,'Value');   
  if plotval == 1
    subplot(2,2,2)
    plot(t,Y)
  else
    subplot(2,2,2) 
    plot(t(1),Y(1), 'EraseMode', 'none')
  end
  
  axis([-1 7 -0.1 0.1])
  
  title('Bus response to a step disturbance of 0.1m')
  xlabel('Time (sec.)')
  ylabel('Bus oscillation (m)')
  hold on

%Draw the bus and ground%  
  l=1;
  h=0.5;

  x1=-l/2;
  y1=Y+0.3;
  x2=-l/2;
  y2=Y+0.3+h;
  x3=l/2;
  y3=Y+0.3+h;
  x4=l/2;
  y4=Y+0.3;
  x5=l/10;
  y5=Y+0.3;
  x6=l/10;
  y6=Y+0.4;
  x7=-l/10;
  y7=Y+0.4;
  x8=-l/10;
  y8=Y+0.3;

  st=0.4+Y;
  sb=0.29+Y;
  tire=0.2;

  gxr=1;
  gy=0.105;
  sxr=-1:0.04:7;
  sxl=-1;
  sy=0.205;

  wx1=-0.35;
  wy1=0.6+Y;
  wx2=-0.35;
  wy2=0.7+Y;
  wx3=-0.25;
  wy3=0.7+Y;
  wx4=-0.25;
  wy4=0.6+Y;
  
  subplot(2,2,4)
  cla
  A = plot([x1 x2 x3 x4 x5 x6 x7 x8 x1],...
         [y1(1) y2(1) y3(1) y4(1) y5(1) y6(1) y7(1) y8(1) y1(1)],...
         'k','EraseMode', 'background','LineWidth',[2]);
  hold on
         
  B = plot([wx1 wx2 wx3 wx4 wx1],[wy1(1) wy2(1) wy3(1) wy4(1) wy1(1)],...         
         'k','EraseMode', 'background','LineWidth',[2]);
         
  C = plot([wx1+0.2 wx2+0.2 wx3+0.2 wx4+0.2 wx1+0.2],[wy1(1) wy2(1) wy3(1) wy4(1) wy1(1)],...         
         'k','EraseMode', 'background','LineWidth',[2]);

  D = plot([wx1+0.4 wx2+0.4 wx3+0.4 wx4+0.4 wx1+0.4],[wy1(1) wy2(1) wy3(1) wy4(1) wy1(1)],...         
         'k','EraseMode', 'background','LineWidth',[2]); 

  E = plot([wx1+0.6 wx2+0.6 wx3+0.6 wx4+0.6 wx1+0.6],[wy1(1) wy2(1) wy3(1) wy4(1) wy1(1)],...         
         'k','EraseMode', 'background','LineWidth',[2]); 
                         
  L = plot([sxl sxr(1) sxr(1) gxr],[sy sy gy gy], 'g',...
                'EraseMode', 'background','LineWidth',[2]);

  I = plot([0 0],[st(1) sb(1)], 'k',...
                'EraseMode', 'background','LineWidth',[2]);
                
  radius = 0.09;
  arcstep = 18;
  j = 0:arcstep:(360-arcstep);  
  arcx = radius * cos((j+arcstep) * pi/180);  
  arcy = radius * sin((j+arcstep) * pi/180);  

  J = patch(arcx, arcy+tire(1), 'b', 'EraseMode', 'background');                                        
                                                                
  axis([-1 1 0 2])
  title('1/4 Bus suspension model animation')

%Check if frames should advance manually%  
  mhandle = findobj('Tag','manbox');
  man=get(mhandle,'Value');
  if man == 1
    'Press any key to advance the animation'
    pause
  end   

%Run the animation% 
  lY=length(Y);
  for z=2:lY-1,
    
    if man == 1
      pause
    end 

    if z == 24
      sb=0.39+Y;
      tire=0.3;
      y1=Y+0.4;
      y4=Y+0.4;
      y2=Y+0.4+h;
      y3=Y+0.4+h;
      y8=Y+0.4;
      y5=Y+0.4;
      y7=Y+0.5;
      y6=Y+0.5;
      st=0.5+Y;
      wy1=0.7+Y;
      wy2=0.8+Y;
      wy3=0.8+Y;
      wy4=0.7+Y;
    end

    set(L, 'XData', [sxl sxr(z) sxr(z) gxr]);
    set(J, 'YData', arcy+tire);
    set(A, 'YData', [y1(z) y2(z) y3(z) y4(z) y5(z) y6(z) y7(z) y8(z) y1(z)]);
    set(B, 'YData', [wy1(z) wy2(z) wy3(z) wy4(z) wy1(z)]);
    set(C, 'YData', [wy1(z) wy2(z) wy3(z) wy4(z) wy1(z)]);
    set(D, 'YData', [wy1(z) wy2(z) wy3(z) wy4(z) wy1(z)]);
    set(E, 'YData', [wy1(z) wy2(z) wy3(z) wy4(z) wy1(z)]); 
    set(I, 'YData', [st(z) sb(1)]);
    drawnow;
    
    if plotval == 0
      subplot(2,2,2)
      plot([t(z),t(z+1)],[Y(z),Y(z+1)], 'EraseMode','none')
    end 
    
  end

%Callback for the RESET button%
elseif action == 2
  
%Reset the Bus to initial conditions%
  subplot(2,2,4)
  l=1;h=0.5;x1=-l/2;y1=0.3;
  x2=-l/2;y2=0.3+h;x3=l/2;
  y3=0.3+h;x4=l/2;y4=0.3;
  x5=l/10;y5=0.3;x6=l/10;
  y6=0.4;x7=-l/10;y7=0.4;
  x8=-l/10;y8=0.3;st=0.4;
  sb=0.29;tire=0.2;gxr=1;
  gy=0.105;sxr=-1;sxl=-1;
  sy=0.205;wx1=-0.35;wy1=0.6;
  wx2=-0.35;wy2=0.7;wx3=-0.25;
  wy3=0.7;wx4=-0.25;wy4=0.6;

  cla
  A = plot([x1 x2 x3 x4 x5 x6 x7 x8 x1],...
         [y1 y2 y3 y4 y5 y6 y7 y8 y1],...
         'k','EraseMode', 'background','LineWidth',[2]);
  hold on
         
  B = plot([wx1 wx2 wx3 wx4 wx1],[wy1 wy2 wy3 wy4 wy1],...         
         'k','EraseMode', 'background','LineWidth',[2]);
         
  C = plot([wx1+0.2 wx2+0.2 wx3+0.2 wx4+0.2 wx1+0.2],[wy1 wy2 wy3 wy4 wy1],...      
         'k','EraseMode', 'background','LineWidth',[2]);

  D = plot([wx1+0.4 wx2+0.4 wx3+0.4 wx4+0.4 wx1+0.4],[wy1 wy2 wy3 wy4 wy1],...      
         'k','EraseMode', 'background','LineWidth',[2]); 

  E = plot([wx1+0.6 wx2+0.6 wx3+0.6 wx4+0.6 wx1+0.6],[wy1 wy2 wy3 wy4 wy1],...      
         'k','EraseMode', 'background','LineWidth',[2]); 
                         
  L = plot([sxl sxr sxr gxr],[sy sy gy gy], 'g',...
         'EraseMode', 'background','LineWidth',[2]);

  I = plot([0 0],[st sb], 'k',...
         'EraseMode', 'background','LineWidth',[2]);
                
  radius = 0.09;
  arcstep = 18;
  j = 0:arcstep:(360-arcstep);  
  arcx = radius * cos((j+arcstep) * pi/180);  
  arcy = radius * sin((j+arcstep) * pi/180);  

  J = patch(arcx, arcy+tire, 'b', 'EraseMode', 'background');                                   
                                                                
  axis([-1 1 0 2])
  title('1/4 Bus suspension model animation')

%Clear the time response plot%
  subplot(2,2,2)
  cla
  axis([-1 7 -0.1 0.1])
  title('Bus response to a step disturbance of 0.1m')
  xlabel('Time (sec.)')
  ylabel('Bus oscillation (m)')

%Callback for the gain slider%
elseif action == 3

%Display current value of the slider%
  gainhandle = findobj('Tag','gainslider');
  gainval=get(gainhandle,'Value');
  curhandle = findobj('Tag','curtext');
  set(curhandle,'String',sprintf('%0.2g',gainval)); 
  
  if  abs(gainval)<1e6
    k=sprintf('%0.1g',gainval);
  else  
    k=sprintf('%0.2g',gainval);
  end
    
  K=eval(k);
  
  Ks=-1e6:1e5:1e7;

  lk=length(Ks);
  for n1=1:lk,
    if Ks(n1) == K
      n2=n1;
    elseif abs(K)<1e5
      n2=11;    
    end
  end
  
  m1=2500;
  m2=320;
  k1 = 80000;
  k2 = 500000;
  b1 = 350;
  b2 = 15020;

  nump=[(m1+m2) b2 k2];
  denp=[(m1*m2) (m1*(b1+b2))+(m2*b1) (m1*(k1+k2))+(m2*k1)+(b1*b2) (b1*k2)+(b2*k1) k1*k2];

  num1=[-(m1*b2) -(m1*k2) 0 0];
  den1=[(m1*m2) (m1*(b1+b2))+(m2*b1) (m1*(k1+k2))+(m2*k1)+(b1*b2) (b1*k2)+(b2*k1) k1*k2];

  numf=num1;
  denf=nump;
          
  z1=3+3.5i;
  z2=3-3.5i;
  p1=30;
  p2=60;
  numc=conv([1 z1],[1 z2]);
  denc=conv([1 p1],[1 p2]);
  
  numk=conv(nump,numc);
  denk=conv(denp,denc);    
      
  R=rlocus(numk,denk,Ks);

%Plot the root locus%
  subplot(2,2,3)
  cla
  rlocus(numk,denk)
  title('Open-loop Root Locus')
  hold on
  RL=plot(real(R(n2,:)),imag(R(n2,:)),'r+');
 
end