function FreeFreeBarFreq (n_e, show)
% n_e = number of elements ; show is one exact mode number to display
% File: FreeFreeBarFreq.m     E*Uxx + w^2*Rho*U = 0, E=1=Rho
% Ref: J.E. Akin, "FEA w Error Estimators" Elsevier, 2005
%.....................................................................
%      x_1=0   x_2    x_3     x_4   x_5= L
%       *------(1)-----*------(2)------*     dU(0)/dx=0=dU(L)/dx
%       T_1    T_2    T_3     T_4     T_5
% Connectivity list:    e     i      j     k
%                       1     1      2     3  % L3
%                       2     3      4     5  % L3
%.....................................................................
% clear all
if ( nargin == 1 ) ; 
  show = 0 ; % flip exact mode shape by using a negative mode number
elseif ( nargin == 0 ) ; 
  show = 0 ; n_e = 2 ; 
end % if from argument count
n_n = 3                     ;  % number of nodes per element
n_g = 1                     ;  % number of DOF per node
%            Set given constants 
n_i = n_n*n_g               ;  % number of DOF per element
n_d = n_e*(n_n - 1) + 1     ;  % system degrees of freedom (DOF)
L = 1.0 ; L_e = L/n_e       ;  % domain and element lengths
E = 1.0 ; Rho  = 1          ;  % stiffness data
X = [0:1:(n_d-1)]*L/(n_d-1) ; % merge lengths

%         Constant quadratic element square matrices
K_E = E   * [ 7, -8,  1, ; ...              % axial
             -8, 16, -8, ; ...
              1, -8,  7 ] / (3*L_e)       ;
M_e = Rho * [4, 2, -1 ; ...                 % mass
             2, 16, 2 ; ...
            -1,  2, 4 ] * L_e/ 30         ;

%         Assemble two n_d by n_d square matrix terms
S = zeros (n_d, n_d) ;  M = zeros (n_d, n_d) ; T(1:n_d) = 0;
for j = 1:n_e                           ; % loop over elements
  rows = [1:n_n] + (j - 1)*(n_n - 1) ;   % get connectivity 
  S (rows, rows) = S (rows, rows) + K_E ; % add to stiffness
  M (rows, rows) = M (rows, rows) + M_e ; % add to mass      
end % for (Note no EBCs for free-free modes)

% Solve for eigenvalues and eigenvectors, general form
[Vec, Diag] = eig(S,M) ; % eigenvector columns, eigenvalue^2 diagonal

% End finite element calculations. ---- Begin graphics ----
 
clf ; hold on ; grid on                    % clear & keep
axis([X(1) X(n_d) -1. 1.]) % non-dimensionalize mode shapes
title ('Quadratic Free-Free Bar Frequencies and Modes')  
xlabel (['Position: (', int2str(n_e), ' elements)']) ; 
ylabel ('Normalized axial displacement (one exact mode in red)')

%  ------------------------------------------------------
%  graph of mode shape in 3 node quadratic line elements
%  ------------------------------------------------------
% e_x, e_y   = coordinates & data values at 3 element nodes
% x_el, y_el = interpolated parametric values

 n_poly = ceil (75/n_e) ; % curve segments
 e_x(3)=0 ; e_y(3)=0 ; H(3)=0 ; % pre-allocation

%  Loop over mode shapes
  root = pi*sqrt(E/(L^2*Rho))       ; % scaling
  Limit = min([ n_d, 3])            ; % limit # of plots <<==========
  for j = 1:Limit                   ; % mode number
    J = n_d + 1 - j                 ; % reverse order
    scale = max(abs(Vec(1:n_d, J))); T = Vec(1:n_d, J) / scale;
    omega = sqrt(real(Diag(J, J))) ;

%    Loop over all elements
    for it = 1:n_e  ;
    
%      Extract nodal coordinates & values. Plot nodes.
      rows = [1, 2, 3] + (it - 1)*(n_n-1); % get connectivity 
      e_x  = X (rows) ; e_y  = T (rows) ;  % gather element data
      plot (e_x, e_y, 'ko')                % nodal symbols
      if (it == 1) % one summary per mode
        e_text = sprintf (' --(Mode %g) Omega = %g', j, omega); 
        text (e_x(2), e_y(2), e_text)     % Plot element number
        plot (e_x(2), e_y(2), 'k*')       % Plot element number
      end % if

%     Loop over local points on the quadratic element
      for k = 1: (n_poly + 1) % points in parametric space
%     get parametric interpolation functions
        Z    = 2*(k - 1)/n_poly - 1 ;   % on -1 to 1
        H(1) = 0.5*(Z*Z - Z) ; H(2) = 1. - Z*Z ;
        H(3) = 0.5*(Z*Z + Z) ;
        x_el (k) = dot(H,e_x) ; y_el (k) = dot(H,e_y) ; % true value
      end % for k parametric points in element
      if ( mod(j, 2) == 0 ) % then alternate colors
        plot (x_el, y_el, 'k-') % parametric curve through 3 nodes
      else
        plot (x_el, y_el, 'b-') % parametric curve through 3 nodes
      end % if color choice
      if ( j == abs(show) ) % show one exact mode
        exact = sign(show)*cos(x_el*(j-1)*root) ; % Thomson "Mech. Vib."
        plot (x_el, exact, 'r-') % exact mode shape at parametric pts
        if ( it == 1 ) % then for first element only
          omega = (j-1)*pi ;
          e_text = sprintf ('----(Mode %g) Exact = %g', j, omega); 
          text (x_el(9), exact(9), e_text)     % Plot 9-th pt
        end % if first element
      end % if show exact

    end % for it over all elements
  end % for j-th mode
% end of function FreeFreeBarFreq (n_e, show)