% hangingBridge.m
% Solution to problem 2 of control 2, Computacion I, 2016/02/08. 
% Model of a hanging bridge.
% J.M.Soler, Feb.2016

% Set parameters
clear
n = 7;     % number of bridge nodes, including end ones
m = 1e2;   % mass of each node (kg)
g = 9.8;   % gravitational acceleration (m/s^2)
k = 5e2;   % spring constant (N/m)
L = 100;   % total bridge length (m)
dv = 5;    % average initial velocities (m/s)
dt = 0.1;  % time interval (s)
tmax = 10; % time span (s)

% Find the equilibrium bridge shape by solving the n equations
%   x(1) = 0
%   -k*(x(i)-x(i-1)) - k*(x(i)-x(i+1)) = 0,    i=2:n-1
%   x(n) = L                 
% whose solution is simply x(i)=(i-1)*L/(n-1). Equally
%   y(1) = 0
%   -k*(y(i)-y(i-1)) - k*(y(i)-y(i+1)) = m*g
%   y(n) = 0
% for the n unknowns y(i). We write these eqs. as a*y=b
a = zeros(n,n);
b = zeros(n,1);
for i = 2:n-1
    a(i,i) = -2*k;
    a(i,i-1) = +k;
    a(i,i+1) = +k;
    b(i) = m*g;
end
a(1,1) = 1;
a(n,n) = 1;
y = a\b;
y = y';                        % row vector
x = linspace(0,L,n);           % uniformly separated x coordinates

% Check that forces are zero for debugging
% f = bridgeForces( [x;y], m, k )

% Plot cable and nodes
plot(x,y,x,y,'o');
axis equal                     % to show real shape
axis([0,L,-50,10])             % consistent with movie window
xlabel('x (m)')
ylabel('y (m)')
pause(1.0)                     % show equil. shape 1s before movie starts

% Set initial conditions, with random velocities
% y = zeros(1,n);              % alternative initial conditions
r0 = [x;y]
v0 = zeros(2,n);
v0(:,2:n-1) = dv*randn(2,n-2);

% Integrate trajectory
myForce = @(r)bridgeForces( r, m, k );
mi = m*ones(1,n);
r = verlet( myForce, mi, r0, v0, dt, tmax);

% Plot trajectory as a movie
t = (0:dt:tmax);               % simulation times
nt = numel(t);                 % number of time frames
for it = 1:nt                  % loop on movie frames
    x = r(1,:,it);             % x coordinates of nodes
    y = r(2,:,it);
    plot(x,y,x,y,'o');
    axis equal
    axis([0,L,-50,10])         % this is to keep window still during movie
    xlabel('x (m)')
    ylabel('y (m)')
    title(['time = ',num2str(t(it)),' s'])   % show time of each frame
    pause(dt);                 % pause for the real time interval
end
