%%animutube
%
%Original version written by Aage Nielsen.
%
%This program makes a dynamic simulation of the movements of a fluid in 
%an u-shaped tube. The principle is a very simple example of the principle 
%in the article of O'Brien and Hodgins. My hope is that understanding this 
%program will make it easier for you to understand the article.

%%
%In the next command I give value to a constant. It is fictive, as the 
%tube is fictive. Think of it as containing the gravity, measures of 
%the tube etc.

k = 0.1

%%
%The next constant c, is the cross-sectional area of the tube.

c = 1

%%
%The constant d, is a damping constant, which is a measure of the 
%resistance against the movements in the tube. We suppose, that the 
%damping is d * the flow.

d = 0.1

%%
%The number r determines the number of timesteps, and accordingly the 
%number of pictures, I intend to make.

r = 200

%%
%Now I give the length of the time intervals. The 'pause' causes 
%the program to stop and wait for you to press any key.

deltaT = 0.5
pause

%%
%Now I make a system of shelves, in mathematics called a matrix. I will 
%denote the horizontal shelves as rows and the vertical as columns. 
%I will use this matrix to store the results of my calculations. From 
%the beginning I put 0's everywhere. When I end a command with ";" the 
%result will not be visible. It would fill a lot on the screen. In the 
%next command I print the first two rows.

M = zeros(r,6);
M(1:2,1:6)
pause

%%
%In the next 6 commands I calculate the numbers corresponding to the 
%starting conditions of the fluid in the tube. I put the results in the 
%first row of the matrix.
%
%In column no 1 I write the height of the fluid in the right branch and 
%in column no 2 I write the height in the left branch.
%
%In column no 3 I write the acceleration in the first timeinterval, 
%which I assume is constant in the timeinterval and equal to k times 
%the height-difference. I chose the sign, such that a downwards 
%acceleration in the right branch of the tube is positive.
%
%In column no 4 I write the initial flow, which is 0.
%
%In column no 5 I write the flow by the end of the first time interval, 
%calculated as deltaT * c * acceleration
%
%In column no 6 I write deltaV, that is the amount of fluid which has 
%been moved in the first timeinterval. See O & H formula (7). It is 
%positive, if fluid is moved from the right to the left branch.

M(1,1) = 4.00;
M(1,2) = 1.00;
M(1,3) = k*(M(1,1)-M(1,2));

M(1,4) = 0;
M(1,5) = deltaT*c*M(1,3);

M(1,6) = deltaT*(M(1,4)+M(1,5))/2;

%%
%We look at the result. I write the first three rows and all 6 columns.

M(1:3,1:6)
pause

%%
%The next block of commands is a loop which calculates the numbers which 
%I want in row no 2 to r in the matrix. In row no i, I want the numbers 
%which give the state in the tube at time i, that is after i-1 time 
%intervals.
%
%In the row no i, column 1 and 2 I put the heights of fluid in the two 
%branches calculated as the heights one timeinterval earlier -/+ the mowed 
%fluid.
%
%In column no 3 the acceleration is calculated. Now a damping term 
%is included.
%
%In column no 4 the flow at the beginning of the timeinterval is 
%calculated. It is simply the flow at the end of the timeinterval one 
%step earlier.
%
%In column no 5 I calculate the flow at the end of the timeinterval; 
%see (6) in O & H.
%
%In column no 6 I calculate deltaV; see (7) in O & H.

for i = 2:r
    M(i,1) = M(i-1,1)-M(i-1,6)/c;
    M(i,2) = M(i-1,2)+M(i-1,6)/c;
    M(i,3) = k*(M(i,1)-M(i,2))-d*M(i-1,5);
    M(i,4) = M(i-1,5);
    M(i,5) = M(i,4)+deltaT*c*M(i,3);
    M(i,6) = deltaT*(M(i,4)+M(i,5))/2;
end

%%
%Let us now take a look on the first 3 rows

M(1:3,1:6)
pause

%%
%We also take a look on the first 20 pair of heights.

M(1:20,1:2)
pause

%%
%and the last 21 pairs

M(r-20:r,1:2)
pause

%%
%It is not necessary for you to understand the details of the last 
%commands in the program. Just think of them as something which I do
%to transform the r pairs of heights to an  animation with r pictures 
%showing the heights as bars representing the fluid heights.

N = [M(:,2) M(:,1)];

figure(1)

for t=1:r
    bar(N(t,1:2), 0.5);
    axis([0 3 0 4.5])
    drawnow
end

