83 result = 3.0/(P*(P-1.0)*(P-2.0)*(P-3.0)) *
84 (H*(H-1.0)*(H-2.0)*(H-3.0)*r*r*(r-1.0)*(r-1.0) +
85 2.0*H*(H-1.0)*(H-2.0)*(Pa*r*(1.0-r)*(1.0-r)+Na*r*r*(1.0-r)) +
86 H*(H-1.0)*(Pa*(Pa-1.0)*(1.0-r)*(1.0-r)+
87 4.0*Na*Pa*r*(1.0-r)+Na*(Na-1.0)*r*r) +
88 2*H*(Pa*Na*(Na-1.0)*r+Pa*(Pa-1.0)*Na*(1.0-r)) +
89 Pa*(Pa-1.0)*Na*(Na-1.0));
91 result /= (6.0*r*r*(1. - r) *(1. - r));
93 return std::max(0.0,result);