function hv = Hv(a,b,D,t)

     F = 1 - D;
     l = (a*b) / (a*F + b*D);

     c2x = 2*l*l*(((b*b)*D+(a*a)*F)/(a*a*b*b)) - 1;

     G = (a*D+b*F)*(a*D+b*F) - (a*a*D + b*b*F);
     A = a*F + b*D;

     vn = (2*l*G/(A*A*A)).*(1-exp(-A.*t)) + (l*c2x).*t;
     dvn = (2*l*G/(A*A)).*exp(-A.*t) + (l*c2x);

     hv = (0.5 .* t) .* (dvn ./ vn);
