function vineplot
%
%  Function: vineplot
%
%  Usage:  vineplot
%
%  Description:  Modified verision of era.m to allow for vine plots 
%  
%
%  Reference for ERA:
%    "An Eigensystem Realization Algorithm for
%    Modal Parameter Identification and Model Reduction," 
%    Juang, J-N, and Pappa, R.S., Journal of Guidance, Control, and
%    Dynamics, Vol. 8, No. 5, pp. 620-627, 1985.
%

%
%  Version SWD970907
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% This matlab source code was originally     %
% developed as part of "DIAMOND" at          %
% Los Alamos National Laboratory. It may     %
% be copied, modified, and distributed in    %
% any form, provided:                        %
%  a) This notice accompanies the files and  %
%     appears near the top of all source     %
%     code files.                            %
%  b) No payment or commercial services are  %
%     received in exchange for the code.     %
%                                            %
% Original copyright is reserved by the      %
% Regents of the University of California,   %
% in addition to Scott W. Doebling, Phillip  %
% J. Cornwell, Erik G. Straser, and Charles  %
% R. Farrar.                                 %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


global MODES DATA

%  Define variables and allocate space for matrices

fprintf(1,' Running Vine Plots ... ')
fprintf(1,'   Forming Hankel matrices...')
   makehank
fprintf(1,'done\n')

nxvec = [MODES.eraparam.nxlow:MODES.eraparam.dnx:MODES.eraparam.nxhigh];

nresp = size(DATA{1}.respDOF,1);
nref = size(DATA{1}.refDOF,1);

% Start loop over model order

for i = 1:length(nxvec),
   
   nx = nxvec(i);
   
   % Form A, B, C, D, V, W
   
   U0N = MODES.eraparam.U0(:,1:nx);
   S0N = MODES.eraparam.S0(1:nx,1:nx);
   V0N = MODES.eraparam.V0(:,1:nx);
   SNsqrt = sqrt(S0N);
   SNinv = sqrt(inv(S0N));
   
   A = SNinv * U0N' * MODES.eraparam.H1 * V0N * SNinv;
   B = SNsqrt * V0N(1:nref,:)';
   C = U0N(1:nresp,:) * SNsqrt;
   D = shiftdim(MODES.eraparam.Y(1,:,:),1);
   
   V = U0N * SNsqrt;
   W = SNsqrt *  V0N';
   
   % Compute emac
   [emac_era,msv_era] = emaccomp(A,V,W,MODES.eraparam.q,MODES.eraparam.d,MODES.eraparam.fs_out);
   
   % Run era2mode
   
   era2mode(A,B,C,MODES.eraparam.fs_out,emac_era,msv_era,DATA{1}.respDOF,DATA{1}.refDOF);
   
   % Fill in freq and cmi vectors for plotting
   temp1 = length(MODES.Freq);
   freq_s(1:temp1,i) = MODES.Freq;
   cmi_s(1:temp1,i) =  MODES.cmi;
   
   % Plot the line
   for j = 1:temp1
      x1 = MODES.Freq{j};
      y1 = nx;
      len = MODES.cmi{j}*MODES.eraparam.dnx;
      line([x1 x1],[y1 y1+len]);
   end
   drawnow
end

fprintf(1,'Done.\n');

return
