function [phi,wn] = eigenModes(geometry,BC,Nmodes)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%                               GOAL

%  Compute eigen frequencies and mode shapes of a beam with different
%  configurations and geometries

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%                               INPUT

% geometry: structure with fields defined in the Example.m file
% BC : 1,2,3 or 4 : Boundaries conditions (cf Example.m)
% Nmodes : Number of modes wished

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%                              OUPUT
%                         ----------------
% Phi : Matrix of mode shapes is a [Nmodes x Nz] vector.
% wn is a  vector of size [1 x Nmodes ] in rad/sec.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%  check input

if ~isnumeric(BC),
    error('BC must be an integer between 1 and 4');
end

if max(ismember([1,2,3,4],BC))==0,
    error('BC is unknown. Please, choose BC = 1,2,3 or 4');
end

if ~isnumeric(Nmodes),
    error('Nmodes must be an integer');
end


name = [{'L'},{'E'},{'nu'},{'rho'},{'I'},{'y'},{'m'}];
for ii=1:numel(name),
    if ~isfield(geometry,name{ii}),
        error([' The field ',name{ii},' is missing in the structure "geometry".'])
    end
end

%%
y= geometry.y;
Nz = numel(y); % number of discrete elements

% volume V and mass m of the beam
m = geometry.m;
L = geometry.L;

% get the non trivial solution of f
Ndummy=1:Nmodes^2; % the number is arbitrary fixed as the square of the number of Nmodes.
h =fsolve(@modeShape,Ndummy);
% small values of h come from numerical errors
h(h<0.1)=[];
% Analytically, many solutions are identical to each other, but numerically, it is not the
% case. Therefore I need to limit the prceision of the solutions to 1e-4.
h = round(h*1e4).*1e-4;
% the uniques solution are called beta:
beta = unique(h);
beta = beta(1:Nmodes);



%modes shapes calculations
if Nmodes>numel(beta),
    error(' Mode shape required is too high \n\n')
end
wn = beta.^2.*sqrt(geometry.E*geometry.I./(m.*geometry.L^4)); % in rad
phi = zeros(Nmodes,Nz);

for ii=1:Nmodes,
    switch BC,
        
        
        case 1,% pinned-pinned
            phi(ii,:)=sin(beta(ii)*y./L);
            
            phi(ii,:)=phi(ii,:)./max(abs(phi(ii,:)));
            
            
        case 2 % clamped-free
            
            phi(ii,:)=(cosh(beta(ii)*y./L)-cos(beta(ii)*y./L))+...
                ((cosh(beta(ii))+cos(beta(ii))).*(sin(beta(ii)*y./L)-sinh(beta(ii)*y./L)))./...
                (sin(beta(ii))+sinh(beta(ii)));
            
            phi(ii,:)=phi(ii,:)./max(abs(phi(ii,:)));
            
            
        case 3 % clamped-clamped
            
            phi(ii,:)=sin(beta(ii)*y./L)-sinh(beta(ii)*y./L)-...
                ((sin(beta(ii))-sinh(beta(ii))).*(cos(beta(ii)*y./L)-cosh(beta(ii)*y./L)))./...
                (cos(beta(ii))-cosh(beta(ii)));
            
            phi(ii,:)=phi(ii,:)./max(abs(phi(ii,:)));
            
            
            
        case 4 % clamped-pinned
            
            phi(ii,:)=sin(beta(ii)*y./L)-sinh(beta(ii)*y./L)-...
                ((sin(beta(ii))-sinh(beta(ii))).*(cos(beta(ii)*y./L)-cosh(beta(ii)*y./L)))./...
                (cos(beta(ii))-cosh(beta(ii)));
            
            phi(ii,:)=phi(ii,:)./max(abs(phi(ii,:)));
    end
    
end



% Function to solve is different for each BC:

    function [f] = modeShape(y)
        switch BC,
            case 1,% pinned-pinned
                f=sin(y).*sinh(y);
            case 2 % clamped-free
                f=cos(y).*cosh(y)+1;
            case 3 % clamped-clamped
                f=cos(y).*cosh(y)-1;
            case 4 % clamped-pinned
                f=cos(y).*sinh(y)-sin(y).*cosh(y);
        end
    end

end