dyst.m

function [A,N,K,lambda,E,omega] = dyst(B,V,F,M)
%
%   A = dyst(B,V,F,M)
%
%  Show the normal dynamical modes of oscillation
%    of a two or three-dimensional structure
%
%  B - 2 column matrix giving the edges
%        the ith row of B says which two nodes are connected by edge i
%  V - matrix whose rows are (x,y) or (x,y,z) coordinates of nodes
%
% Optional arguments:
%
%  F - row vector of fixed nodes (none is default)
%  M - row vector of masses      (unit masses is default)
%
% Output:
%
%  A - incidence matrix for structure
%  N = null(A,'r') - basis for instabilities
%  K - stiffness matrix for structure
%  lambda - eigenvalues of K
%  E - eigenvectors of K
%  omega = sqrt(lambda) - fundamental frequencies
%
%   See also STRUT, STRUTMECH, MOTION, MOTIONS, MSPR,

if nargin < 2, error('Not enough input arguments.'); end
if nargin < 3, F = []; end

global buttonid;

buttonid = 0;

icons = ['text(.5,.5,''pause''   ,''Horiz'',''center'')'
         'text(.5,.5,''input''   ,''Horiz'',''center'')'
         'text(.5,.5,''stop''    ,''Horiz'',''center'')'];
callbacks = ['button(1)'
             'button(2)'
             'button(3)'];
presstype =  ['toggle'
              'flash '
              'flash '];

hf = figure(1); 
clf;

bg = btngroup(hf,'GroupID','b','ButtonId',['p';'i';'s'],...
            'IconFunctions',icons,'GroupSize',[1 3],...
	    'CallBack',callbacks,'Position',[.3 0 .4 .07],...
            'PressType',presstype);


pa = 25; tol = .00001;

rad = .08; dt = .5;

nedge = size(B,1);
[nvert, d] = size(V);
U = V';
G = setdiff([1:nvert],F);

A = zeros(nedge,d*nvert);
Nrepmat = NaN; Nrepmat = Nrepmat(ones(nedge,1));

minV = min(V); maxV = max(V);
dV = max(V) - min(V); dmax = max(dV);
delV = (dmax - dV)/2;
mm = [minV - delV - .5;maxV + delV + .5];


for i=1:nedge,
  m = B(i,1); n = B(i,2);
  x1 = V(m,:); x2 = V(n,:);
  unit = (x2 - x1)/norm(x2 - x1);
  vm = [d*(m-1)+1:d*m]; vn = [d*(n-1)+1:d*n];
  A(i,vm) = - unit;
  A(i,vn) =   unit;
end

if nargin < 4, M = ones(1,nvert); 
elseif size(M,2) ~= nvert
	error('Not enough masses');
end


nfix = size(F,2); nmove = nvert - nfix;
fix = fliplr(sort(F));
nmode = d*nmove;

for i=1:nfix
  m = fix(i);
  x1 = V(m,:); 
  vm = [d*(m-1)+1:d*m];
  A(:,vm) = [];
  M(m) = [];
end

Ma = repmat(M,d,1);  Ma = diag(Ma(:));
Mi = inv(sqrt(Ma));

N = null(A,'r');
n = size(N,2);
disp(' ')
if n == 0
  disp('Structure is rigid!')
  disp(' ')
else
  disp(['Number of motions = ',num2str(n)])
end

if size(B,2) == 3, C = diag(B(:,3)); else C = eye(nedge); end
K = Mi * A' * C * A * Mi;

[E,D] = eig(K);
lambda = diag(D);
[lambda,lperm] = sort(lambda);
zerol = abs(lambda) < tol;
lambda = (1 - zerol) .* lambda;
D = diag(lambda);
E = E(:,lperm); E(:,1:n) = N;
omega = sqrt(lambda);

disp('Incidence Matrix');
disp(' ');
disp(A);
disp('Stiffness Matrix');
disp(' ');
disp(K);

disp('Eigenvalues');
disp(' ');
disp(lambda');
disp('Eigenvectors');
disp(' ');
disp(E);
disp('Frequencies');
disp(' ');
disp(omega');

j = 1; r = .3;

U0 = U;
U0(:,F) = [];
U0 = U0(:);

rst = 1; j = 1;

ax = axes('Position',[.15 .15 .75 .75]);

Ux = U(1,:); X = [ Ux(B) Nrepmat]'; X = X(:);
Uy = U(2,:); Y = [ Uy(B) Nrepmat]'; Y = Y(:);
if d==3 
Uz = U(3,:); Z = [ Uz(B) Nrepmat]'; Z = Z(:);
end	

if d==2 
	hpl=plot(X,Y,'EraseMode','background');
	else
	hpl=plot3(X,Y,Z,'EraseMode','background');
end

axis(mm(:));

popupStr = reshape(['comb -----  ' num2str(1:nmode,'%5i')  '  size '],5,nmode+3)';
popupStr(logical([0;0;zerol]),5) = '*';

popupHndl=uicontrol( ...
      'Style','text', ...
      'HorizontalAlignment','left',...
      'Units','normalized', ...
      'Position',[.93 .5 .07 .07], ...
      'String','mode');
popupHndl=uicontrol( ...
      'Style','popup', ...
      'Units','normalized', ...
      'Position',[.93 .465 .07 .07], ...
      'String',popupStr);

jmode=2;
set(popupHndl,'Value',jmode)
j = 1;

if zerol(j) == 0, fq = 'unstable'; else fq = num2str(omega(j)); end
      title(['Mode:  ',num2str(j),'       Size:  ',num2str(r),...
	         '       Frequency:   ',fq],...
			 'Units','normalized','Position',[.5 1.05]);

stopit = 0;

while stopit == 0 
  switch buttonid
   case 1  
	  buttonid = 0;
      while buttonid == 0, pause(.2); end 
	  buttonid = 0;
   case 2
      rst = 1;
	  disp(['Currently j = ',num2str(j),'  r = ',num2str(r),...
	           '  dt = ',num2str(dt)]);
      keyboard
	  buttonid = 0;
   case 3
	  stopit = 1;
	  buttonid = 0;
  end

 if jmode~=get(popupHndl,'Value'),
    jmode=get(popupHndl,'Value');
    rst = 1;
    if or(jmode == 1,jmode == 2),
	   jo = input('Input mode combination :');
    	if size(jo,2) ~= 1 & size(jo,2) ~= nmode 
     	  disp(['Invalid: combination must be a row vector with ',...
               num2str(nmode),' entries.']);
     	  rst = 0;
		else 
     	  j = jo;
		end
       jmode=2;
	   set(popupHndl,'Value',jmode)
	else
    if jmode == nmode+3
	   r = input('Input size of mode :');
       jmode=2;
	   set(popupHndl,'Value',jmode)
	else
	   j = jmode - 2;
	end
	end
 end
			
			
 if rst == 1
    if size(j,2) == 1 
      if zerol(j) ~= 0, fq = 'unstable'; else fq = num2str(omega(j)); end
      title(['Mode:  ',num2str(j),'       Size:  ',num2str(r),...
	         '       Frequency:   ',fq],...
			 'Units','normalized','Position',[.5 1.05]);
      dj = zeros(nmode); dj(j,j) = r;
    else
      title(['Combination mode: ',num2str(j)]);
      dj = diag(j);
    end
  Ej = E * dj; 
  s = 0;
  rst = 0;
 end

S = U0 + Ej * (sin(omega*s) + zerol*s);
S = reshape(S,d,nmove);
U(:,G) = S;

Ux = U(1,:); X = [ Ux(B) Nrepmat]'; X = X(:);
Uy = U(2,:); Y = [ Uy(B) Nrepmat]'; Y = Y(:);
if d==3 
Uz = U(3,:); Z = [ Uz(B) Nrepmat]'; Z = Z(:);
end	



if d==2 
	set(hpl,'xdata',X,'ydata',Y);
	else
	set(hpl,'xdata',X,'ydata',Y,'zdata',Z);
end
drawnow

s = s + dt;

k = 0;
for i=1:1000*pa, k = 1-k; end

end