function [v] = v2fp(s, x, y, beta, a)

	nmx2  = norm(x)^2;
	nmy2  = norm(y)^2;
	nmxy2 = norm(x-y)^2;

	A = nmx2 / 8 / beta;
	B = nmy2 / 8 / beta;

	T1 = 1.0 / sqrt(8*pi*beta);
	T2 = exp(nmxy2/8/beta);
	T1_a_T2 = T1 * a * T2;
	T3 = 1.0 ./ ((s.*(1-s)).^(1.5));
	T4 = exp(-A./s);
	T5 = exp(-B./(1-s));
	v = exp((A-A.*s+B.*s)./((s-1).*s)) .* ...
		sqrt((1-s).*s) .* ...
		(2*A.*(s-1).^2 - s.*(6*s.^2+(2*B-9).*s+3)) ./ ...
		(2 * (s-1).^4 .* s.^4);
	v = v * T1_a_T2;
	return
end
