function phv = plot3D_rph2_hv_a_tau(a,b,w2,tau)
% Compute the Variability Index versu w2 and tau
% a and tau are row vectors

n = length(tau);

w1 = 1 - w2;

l = (a.*b) ./ (a.*w2 + b.*w1);

C2x = 2*(l.*l).*(((a.*a).*w2 + (b*b).*w1)./((a.*a)*b*b)) - 1;

A = (2*l.*((a*w1+b*w2).^2 - ((a.*a).*w1 + b*b*w2)))./(a*w2 + b*w1).^3;

for i=1:1:n
vn = A.*(1-exp(-1*(a*w2+b*w1)*tau(i))) + l.*C2x*tau(i);
dvn = A.*(a*w2+b*w1).*exp(-1*(a*w2+b*w1)*tau(i)) + l.*C2x;

hv(:,i) = (0.5 * tau(i)) * (dvn ./ vn)';
end

hv;

mesh(tau, a, hv);
xlabel('\tau (s)','FontSize',14);
ylabel('a','FontSize',14);
zlabel('H_v','FontSize',14);
