function [xt,yt,zt] = channel_waterfall( channel_width, height, u0, N )

% channel_waterfall( channel_width, height, u0, N ) returns a series of 
% N trajectories for particles exiting a square channel of specificied width
% and with a maximum velocity u0 and falling height distace. All units are
% meters or meters/sec.

% The starting locations of the particles. All start at x = 0
yvec = -( channel_width / N ) * (0:N);
zvec = -( channel_width / N ) * (0:N);
[ y, z ] = meshgrid( yvec, zvec );
x = zeros( size( y ) );

% The angle between the horizontal and the channel
slope_angle = 0;

% velocities in the x, y, z directions (respectively)
[ u, v, w ] = channel_velocity( y, z, width, slope_angle, u0 );


% Time step
delta_t = .01;

figure
hold on

% Calculate the individual trajectories
for i = 1:(N+1)
    for j = 1:(N+1)
        % Calculate the trajectory
        pos = [ x( i, j ) y( i, j ) z( i, j ) ];
        vel = [ u( i, j ) v( i, j ) w( i, j ) ];
        [ x_traj, y_traj, z_traj, t ] = freefall( pos, vel, height, delta_t );

        % Determine the size of the output matrices by length of trajectory
        % vectors.
        if j == 1
            len = length( x_traj );
            xt = zeros( N+1, len );
            yt = zeros( N+1, len );
            zt = zeros( N+1, len );
        end
        
        % store the trajectory
        xt( j, : ) = x_traj';
        yt( j, : ) = y_traj';
        zt( j, : ) = z_traj';
    end

    % plot all the trajectories with the same initial z
    surf( xt, yt, zt )
    shading interp
end

view( -15, -50 )
xlabel( 'x' )
ylabel( 'y' )
zlabel( 'z' )