function [Yr,bars,L,N,H]=trussCantilever(N,Aval)
% Yr = trussCantilever(N,Aval);
%
% Usage:
%    Y=trussCantilever(N,Aval)
%    [Yr,bars,L,N,H]=trussCantilever(N,Aval)
% 
% Inputs:
%
%    N: number of horizontal elements (recommend range of [1,400])
%
%    Aval: cross section area for each rod (recommend range of [10 1000])
%
% Outputs:
%
%    Yr: 
%
%

% Copyright 2015 The MathWorks, Inc


%% Defaults

if nargin < 2
    % Ten horizontal truss elements in bottom cord 
    N = 10;
    
    % Cross sectional area of 100
    Aval = 100;
end


%% Truss Parameters
% Physical parameters of the material and truss

% Horizontal Length of the truss
L = 1;

% Vertical Height of the truss
H = 0.025;

% Density of truss bar material
rho = 1;

% Modulus of Elasticity of truss bar material
Eval = 1e-1;

% Rayleigh, "alpha,beta" damping coefficients
dampingCoeffAlpha = 0.01;
dampingCoeffBeta = 0.01;


%% Force Parameters
%
% The following downward force is applied to the truss
%
%      max(0,(applyTime-t)*downwardForceMag/applyTime)*sin(2*pi*freq*t)
% 
% The frequency of the force is 1. The magnitude is linearly lowered from
% _downwardForceMag_ to zero over the interval |[0,applyTime]|. The
% simulation time is simTime. Because of the damping the vibration will die
% out.

downwardForceMag = 0.02;
freq = 1;
applyTime = 20;
simTime = 100;


%% Map reduced Dofs to actual Dofs
% The truss is fixed to the wall on the left hand side. Therefore the two
% nodes there are not mobile and can be eliminated. Create a map between
% the "reduced" degrees of freedom and the actual degrees of freedom

numDofs = 2*2*(N+1)-2; % -2 because the cantilever is pointed upwards at the end
groundDofs = [1,2,2*(N+1)+1,2*(N+1)+2]; % Degrees of Freedom that will be eliminated

%% Add bars along with geometry
% bars are #of Bars X [Area,E,length of bar,angle of bar w.r.t. horizontal
% axis, node1 (from node), node2 (to node)]
bars = zeros(2*N+2*(N-1),6);
for n = 1:N
    % upper bars
    lelem = L/N;
    bars(n,:) = [Aval,Eval,lelem,0,n,n+1];
    % diagonal bars
    lelem = sqrt((L/N)^2+H^2);
    bars(N+n,:) = [Aval,Eval,lelem,pi/4,N+1+n,n+1];
end
for n = 1:N-1
    % lower bars
    lelem = L/N;
    bars(2*N+n,:) = [Aval,Eval,lelem,0,N+1+n,N+1+n+1];
    % vertical bars
    lelem = H;
    bars(2*N+N-1+n,:) = [Aval,Eval,lelem,pi/2,N+1+n+1,n+1];
end

%% Assemble all bars into global stiffness and mass matrices

% Use a dense matrix for insertion efficiency
K = zeros(numDofs);
M = zeros(numDofs);

% Mass element matrix
% the mass matrix is a lumped matrix with half the mass going to one
% node and half to the other node
unitMassMatrix = diag([1/2 1/2 1/2 1/2]);
melem = rho*lelem*unitMassMatrix; 

for ii=1:size(bars,1)
    % extract parameters for stiffness and mass matrices
    Aelem = bars(ii,1);
    Eelem = bars(ii,2);
    lelem = bars(ii,3); 
    telem = bars(ii,4);
    kelem = localStiffness(Aelem,Eelem,lelem,telem); % stiffness matrix
    n1 = bars(ii,5); 
    n2 = bars(ii,6);
    % convert to reduced dofs
    ix = [n1*2-1,2*n1,n2*2-1,n2*2];
    % element "stamping"
    K(ix,ix) = K(ix,ix) + kelem; % Add in element's contribution to K
    M(ix,ix) = M(ix,ix) + melem; % Add in element's contribution to M
end

% Zero out rows and columns of ground DOFs in stiffness and mass matrix
K(groundDofs,:) = 0;
K(:,groundDofs) = 0;
M(groundDofs,:) = 0;
M(:,groundDofs) = 0;

% Add ones to diagonal of ground DOFs to freeze them
groundDofDiag = sub2ind(size(M),groundDofs,groundDofs);
K(groundDofDiag) = 1;
M(groundDofDiag) = 1;

% convert stiffness and mass matrices to sparse
Kr = sparse(K);
Mr = sparse(M);

% Force vector iz zeros with magnitude on end Y Dof in Y
F = zeros(size(Kr,1),1);
F(2*(N+1)) = downwardForceMag;

%% Solve the system
% Transform M*d2x/dt2 + damping*dx/dt + K*x = F into dydt = A*y + f form
% Convert 2nd order ODE to first order ODE with y = dx/dt transformation.
% The state space is therefore [position,velocity] and in that order

[LM,UM] = lu(Mr);
A = [-UM\(LM\Kr),-UM\(LM\(dampingCoeffAlpha*Kr+dampingCoeffBeta*Mr));sparse(numDofs,numDofs),speye(numDofs)];
f = [UM\(LM\F);zeros(numDofs,1)];
myEvaluator = @(t,y) A*y + max(0,(applyTime-t)*f/applyTime)*sin(2*pi*freq*t);
initialCondition = zeros(numDofs*2,1); % initial conditions for ode23t

% Yr is #Time points X [position, velocity]
[~,Yr] = ode23t(myEvaluator,[0 simTime],initialCondition);

% Remove the velocity components from Yr
Yr = Yr(:,1:numDofs);

end