% ========= Sample Matlab FEA eigensolution =================
% Previously:
n_freq  = 3  ;   % number of frequencies and modes to compute
show_v  = 3  ;     % number of eigenvectors to list <= n_freq

% identify the free and fixed degrees of freedom
[Fixed] = get_fixed_from_ebc_flags (n_d, EBC_flag)   ; % EBCs
[Free]  = get_free_from_fixed (n_d, Fixed)    ; % to be found

% ===========================================================
%   PHASE 5: SOLVE EIGENVALUE SYSTEM, PRETTY PRINT & SAVE
% ===========================================================
fprintf ('Phase 5, solve modified system eigenproblem \n')

% Solve for smallest eigenalues and eigenvectors
fprintf ('\n'); fprintf ('Beginning eigen-solution \n'); 
[V, D] = eigs (S(Free, Free), M(Free, Free), n_freq, 'sm'); 
V = real(V); F_sq = real (diag(D))        ; % real frequency^2 

%*************************************************************
% NOTE: The D is not in increasing order. They must be sorted.
%       Order(1:n_freq) is a vector subscript that gives the 
%       lowest n_freq eigenvalues (frequences squared) as the 
%       vector F_sq(Order(1:n_freq))
%*************************************************************
[F_sq, Order] = sort (F_sq); % reorder eigenvalues^2, lo to hi

%*************************************************************
% NOTE: Now the eigenvalues are in the correct order, but not 
%       the eigenvectors. The Order array is needed to extract
%       the columns with lowest n_freq eigenvectors desired !
%*************************************************************

fprintf ('\n')        ; % list the eigenvalues and frequencies
fprintf ('NOTE: Rigid body motions give zero frequencies \n');
kount = 0                       ; % number of true eigenvalues
for k = 1:n_freq ; % list eigenvalues only --> --> --> --> --> 
  kount = kount + 1           ; % number of output eigenvalues
  fprintf('Eigenvalue %i = %9.4e \n', k, F_sq(k)); % reference 
  Freq (k)  = sqrt(F_sq(k))         ; % frequency, radians/sec
  Hertz (k) = Freq (k) / (2*pi)      ; % frequency, cycles/sec
  fprintf ('Natural frequency (Hz, rps) %i %9.4e %9.4e \n',...
            kount, Hertz (k), Freq (k) )      ; % print values
end ; % for k frequencies for output   <-- <-- <-- <-- <-- <-- 

fid = fopen('mode_value_vector.txt', 'w')    ; % e-v save file 
%*************************************************************
%NOTE: file mode_value_vector.txt has one line per eigenvector
%      component preceeded by its eigenvalue_number, and its
%      eigenvalue.
%      It is used by other plot programs to show mode shapes.
%*************************************************************
fprintf ('\n') ; kount = 0      ; % number of true eigenvalues
V_Full = zeros (n_d, 1)               ; % allocate work vector
for K = 1:n_freq ; % list requested eigenvalues ===> ===> ===>
  kount  = kount + 1            ; % number of true eigenvalues
  actual = Order(K)                 ; % where V is this column 
  V_Free = V(:, actual)             ; % free dof for this mode
  V_Full (Free) = V_Free (:)       ; % expand to original size
  %   Now V_Full has both the EBCs and the computed mode shape
  biggest = max ( max ( abs (V_Full) ) )        ; % scale them

  %B % send to mode surface plot with specific name
  %B mode_shape_surface (x, y, V_Full, nodes, Hertz(K), K, 1)

  if (K <= show_v); fprintf ('\n'); % list show_v eigenvectors
    fprintf ('Frequency (Hz) %i = %9.4e \n', kount, Hertz(K));
    fprintf ('Eigenvector (mode shape): \n') ; % vector header
    fprintf ('node, component_values \n')    ; % vector values
    for k = 1:n_m ; % loop over all nodes  --> --> --> --> -->
       fprintf ('%i  %9.4e  \n', k, V_Full(k)) ;
    end; % for k nodes <-- <-- <-- <-- <-- <-- <-- <-- <-- <--
  end                          ; % if list show_v eigenvectors

% Save V_Full for mode shape plots
  Hz_K = Hertz (K)                      ; % data for save file
  for k = 1:n_m ; % loop over all nodes    --> --> --> --> -->
    fprintf (fid,'%i  %g  %g  \n', ...
             kount, Hz_K, V_Full(k)) ;
  end ; % for save node k values   <-- <-- <-- <-- <-- <-- <--

  if ( kount == n_freq )     ; % have finished requested range
    break                       ; % break out of this for loop
  end                ; % if finished requested number of modes
end ; % for K modes of intetest  <=== <=== <=== <=== <=== <===
% ================== END PHASE 5 =============================

% L-shaped membrane test
% Natural frequencies for 5 requested modes 
% Mode  Hertz   Scaled 
% 1  4.284e+02  3.704e-01