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