function [Ma_return] = Ma_calcs(lb_Q_above_base, rb_Q_above_base, ...
    lb_Q_with_base, rb_Q_with_base, t_steps, lb_cell, rb_cell) 
% Function organizes and makes calculations for mass flux along
% channel-floodplain interface, given discharge outputs

lb_length = cell2mat(lb_cell(13,1));
rb_length = cell2mat(rb_cell(13,1));

% (kg) calculate mass flux of water transferred to left WITH BASE
lb_Ma_raw = lb_Q_with_base * 1000; 
% (kg) calculate mass flux of water transferred to right WITH BASE
rb_Ma_raw = rb_Q_with_base * 1000;

% (kg) calculate mass flux of water transferred to left NO BASE
lb_Ma = lb_Q_above_base * 1000; 
% (kg) calculate mass flux of water transferred to right NO BASE
rb_Ma = rb_Q_above_base * 1000;

[num_lb,~] = size(lb_Ma);
[num_rb,~] = size(rb_Ma);

lb_Ma_per_l = zeros(num_lb,t_steps);
rb_Ma_per_l = zeros(num_rb,t_steps);

% Calculate mass flux per unit length
for i = 1:num_lb
    lb_Ma_per_l(i,:) = lb_Ma(i,:) ./ lb_length(i,1);
end 

for j = 1:num_rb
    rb_Ma_per_l(j,:) = rb_Ma(j,:) ./ rb_length(j,1);
end

% Calculate mass flux (water) crossing fp boundaries at each time step
lb_Ma_t = sum(lb_Ma);
rb_Ma_t = sum(rb_Ma);

% Calculate positive and negative flux at left bank fp interface for
% each time step
lb_Ma_pos = zeros(num_lb,t_steps);
lb_Ma_neg = zeros(num_lb,t_steps);

for Ma_i = 1:num_lb
    for Ma_j = 1:t_steps
        if lb_Ma(Ma_i,Ma_j) >= 0
            lb_Ma_pos(Ma_i,Ma_j) = lb_Ma(Ma_i,Ma_j);
            lb_Ma_neg(Ma_i,Ma_j) = 0;
        elseif lb_Ma(Ma_i,Ma_j) < 0
            lb_Ma_pos(Ma_i,Ma_j) = 0;
            lb_Ma_neg(Ma_i,Ma_j) = lb_Ma(Ma_i,Ma_j);
        end
    end
end

lb_Ma_t_pos = sum(lb_Ma_pos);
lb_Ma_t_neg = sum(lb_Ma_neg);

% Calculate positive and negative flux at right bank fp interface for each
% time step
rb_Ma_pos = zeros(num_rb,t_steps);
rb_Ma_neg = zeros(num_rb,t_steps);

for Ma_i = 1:num_rb
    for Ma_j = 1:t_steps
        if rb_Ma(Ma_i,Ma_j) >= 0
            rb_Ma_pos(Ma_i,Ma_j) = rb_Ma(Ma_i,Ma_j);
            rb_Ma_neg(Ma_i,Ma_j) = 0;
        elseif rb_Ma(Ma_i,Ma_j) < 0
            rb_Ma_pos(Ma_i,Ma_j) = 0;
            rb_Ma_neg(Ma_i,Ma_j) = rb_Ma(Ma_i,Ma_j);
        end
    end
end

rb_Ma_t_pos = sum(rb_Ma_pos);
rb_Ma_t_neg = sum(rb_Ma_neg);

Ma_return = cell(1,6);
Ma_return(1,1) = {lb_Ma_raw};
Ma_return(1,2) = {rb_Ma_raw};
Ma_return(1,3) = {lb_Ma};
Ma_return(1,4) = {rb_Ma};
Ma_return(1,5) = {lb_Ma_per_l};
Ma_return(1,6) = {rb_Ma_per_l};
Ma_return(1,7) = {lb_Ma_t};
Ma_return(1,8) = {rb_Ma_t};
Ma_return(1,9) = {lb_Ma_t_pos};
Ma_return(1,10) = {lb_Ma_t_neg};
Ma_return(1,11) = {rb_Ma_t_pos};
Ma_return(1,12) = {rb_Ma_t_neg};

end