function [f] = v2f(r1, r2, theta, beta, a, s)
% Integrand for the two-jump interaction potential:
%   function [f] = v2f(r1, r2, theta, s, beta, a).
% Return value f is a matrix in r1, r2, theta, beta, a, and s.
% Integrate by summing over s.

	[ignored, nr1]    = size(r1);
	[ignored, nr2]    = size(r2);
	[ignored, ntheta] = size(theta);
	[ignored, nbeta ] = size(beta);
	[ignored, na]     = size(a);
	[ignored, ns]     = size(s);

	f = zeros(nr1, nr2, ntheta, nbeta, na, ns);

	for ir1 = 1 : nr1
		sr1 = r1(ir1);

		for ir2 = 1 : nr2
			sr2 = r2(ir2);

			for itheta = 1 : ntheta
				stheta = theta(itheta);

				for ibeta = 1: nbeta
					sbeta = beta(ibeta);

					for ia = 1: na
						sa = a(ia);

						for is = 1 : ns
							ss = s(is);

							f(ir1, ir2, itheta, ibeta, ia, is) = ...
								v2f_scalar(sr1, sr2, stheta, sbeta, sa, ss);
						end
				end
			end
		end
	end
end

%   r 1    2         3     4
% s
% 0.1 0.1982 0.0038 0.0000 0.0000
% 0.2 0.1535 0.0184 0.0005 0.0000
% 0.3 0.1230 0.0257 0.0019 0.0000
