function [ub] = flowfield_Schwenk(B,C,Do,s,alphaa,F2,Cf_ch,truncate_conv)

% 11/19/2014 - Did speed tests on for computing integral two additional
% ways: 1) over-allocating the integrands vector, padding the extras with
% zeros so that Matlab doesn't need to keep creating a new one each
% iteration. 2) Computing the convolution in a for loop and truncating when
% the additional contribution falls below some threshold. Both methods were
% longer. The for loop was exceptionally longer.

% compute parameters and nondimensionalize variables
s = s/B;
Ro = min(abs(1./C)); % minimum radius of curvature along reach
C = C*Ro; 
beta = B/Do; % (half-width) aspect ratio
lam_o = -2*beta*Cf_ch; % characteristic exponent describing upstream influence
nu_o = (B/Ro);
A = (alphaa + 2); % accounts for secondary flow, see Sun 1996 among many others

% upstream boundary condition
ub(1) = 0;
% downstream boundary condition


if truncate_conv == 0
    % compute convolution integral in loop
    int_term = NaN(numel(C),1);
    for sidx = 1:numel(s)
        if sidx == 1
            int_term(sidx) = 0;
        else
            integrands = C(1:sidx).*exp(lam_o*(s(sidx)-s(1:sidx)));
            S_axis = s(1:sidx);
            m = numel(S_axis);
            int_term(sidx) = diff(S_axis)' * (integrands(1:m-1) + integrands(2:m))/2;
        end
    end

else
    % compute convolution integral in loop, but truncate convoluation after
    % some threshold distance
    conv_int_thresh = 100; % multiples of B: threshold distance to compute convolution integral (saves computation time)
    int_term = NaN(numel(C),1);
    for ds_idx = 1:numel(s)
        if ds_idx == 1
            int_term(ds_idx) = 0;
        else
            s_ds = s(ds_idx);
            if s_ds > conv_int_thresh
                s_us = s_ds - conv_int_thresh;
                [~, us_idx] = min(abs(s-s_us)); % find the index that is conv_int_thresh*B away from the node in question
            else
                us_idx = 1;
            end

            integrands = C(us_idx:ds_idx).*exp(lam_o*(s(ds_idx)-s(us_idx:ds_idx)));
            S_axis = s(us_idx:ds_idx);
            m = numel(S_axis);
            int_term(ds_idx) = diff(S_axis)' * (integrands(1:m-1) + integrands(2:m))/2;
        end
    end
end

ub = (ub(1)+nu_o*(C(1))).*exp(lam_o*s) + nu_o.*(-C - lam_o/2.*(F2 + A).*int_term); % equation 2 in Schwenk 2014



end % function
