kbeam = 1e3; mbeam = 2;
alpha = 0.01; beta = alpha;
Nnodes = 29; % number of nodes of the bridge (odd)
m = 2 * Nnodes;
[ K ] = bridge_stiffness_matrix( Nnodes, kbeam );
M = mbeam * speye( m, m );
C = alpha * M + beta * K;
plot_bridge( Nnodes );

%%

i_t = [ 3 : m - 1 ];
m_t = m - 3;
M_t = M( i_t, i_t ); C_t = C( i_t, i_t ); K_t = K( i_t, i_t );

%%

I_t = speye( m_t, m_t );
P_t = [ M_t, sparse( m_t, m_t ); sparse( m_t, m_t ), I_t ];
A_t = P_t \ [ - C_t, - K_t; I_t, sparse( m_t, m_t ) ];

%%

node_force = 15; % node in which the force is applied
b = zeros( m, 1 ); b( 2 * node_force, 1 ) = -1;
b_t = b( i_t, 1 );
y0_t = zeros( 2 * m_t, 1 );
t0 = 0; tf = 250;

t_ref = 25;
g_t = @( t ) P_t \ [ b_t * ( ( t / t_ref ) * ( t <= t_ref ) + ...
                               1 * ( t > t_ref ) ); zeros( m_t, 1 ) ];
Nh = 2500;
h = ( tf - t0 ) / Nh;
[ tv, uv ] = backward_euler_system_nhcc( A_t, g_t, y0_t, t0, tf, Nh );
d_t = uv( m_t + 1 : end, : );
d = zeros( m, Nh + 1 );
d( i_t, : ) = d_t( :, : ); % displacement vector (including the constraints)
inode = 15; % node in which we are interested to evaluate the displacement
u_inode_x = d( 2 * inode - 1, : ); u_inode_y = d( 2 * inode, : );
plot( tv, u_inode_x, '-b', tv, u_inode_y, '-r' );
grid on; xlabel('t'); ylabel('disp');
legend('disp x','disp y' );

%%

g_t = @( t ) P_t \ [ b_t * ( ( t / t_ref ) * ( t <= t_ref ) ); zeros( m_t, 1 ) ];

%%

omega = 0.25;
g_t = @( t ) P_t \ [ b_t * sin( omega * t ); zeros( m_t, 1 ) ];

%%
omega_natural = sort( sqrt( eigs( K_t, M_t, 5, 'SM' )' ) )
% omega_natural = 
%     0.4688    1.6652    2.3342    3.6036    5.5365

%%

for i = 1 : 10 : Nh + 1
    plot_bridge( Nnodes, d( :, i ) );
    pause( 0.1 )
end