function [x_traj,y_traj,z_traj] = mult_channel_waterfall( num_channels, width, height, u0, N, edge_shape )

% mult_channel_waterfall( num_channels, width, height, u0, N, edge_shape ) 
% plots the waterfall out of a waterfall with num_channels, each width by
% width. The maximum velocity (i.e. center of the pipe) is u0. There are NxN
% trajectories from each channel. The return values are the trajectories
% from one channel, [x_traj,y_traj,z_traj]. The function handle edge_shape,
% descripes the overall shape of the edge of the waterfall. This is
% evaluated at the center of each channel, which are assumed to have a
% square end. The particles fall from z=0 to z=-height. Units are meters or
% meters/sec, as appropriate.

% location of the center of the channels
channel_centers_y = (.5:num_channels) * width;
channel_centers_x = feval( edge_shape, channel_centers_y, width, u0 );

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

% The starting locations of the particles. All start at x =
% channel_centers_x
yvec = -( width / N ) * (0:N);
zvec = -( width / N ) * (0:N);
[ y_def, z_def ] = meshgrid( yvec, zvec );
x_def = zeros( size( y_def ) );

% velocity field, x,y,z components
[ u, v, w ] = channel_velocity( y_def, z_def, width, slope_angle, u0 );

figure
hold on

% Time step
delta_t = .01;

% Calculate the individual trajectories
for k = 1:num_channels
    % specify x, y, z based on the channel we're in
    x = x_def + channel_centers_x( k );
    y = y_def + channel_centers_y( k );
    z = z_def;
    for i = 1:(N+1)
        for j = 1:(N+1)
            % compute a 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 the surface of all the trajectories with the same initial z.
        surf( xt, yt, zt )
        shading interp
    end
end

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