Numerical simulation for density of strings in a ball

Question

I am trying to numerically reproduce the Mathematica result for the density of strings in a sphere as explained here.

1. The program is to generate random, uniformly distributed pair points on a unit sphere. This pair forms a line segment.

2. Use shells of constant thickness, incrementally increasing in radius, to calculate the density of lines at that radius by finding the portion of all line segments that lie in the shell.

3. Plot the density as a function of r.

Code:

noOfTrials = 500;
points = RandomPoint[Sphere[], {noOfTrials, 2}];
(*gives the length of a line from pt1 to pt2 that lies inside a 
 sphere at origin of radius r*)
dFunc[pt1_, pt2_, r_] := Module[
  {midPt, h},  
  midPt = (pt1 + pt2)*0.5;
  h = Norm[midPt];
  If[h >= r, 0, 2 N@Sqrt[r^2 - h^2]]
]  
(*simulation*)
sim = Module[
{dr, l, dlTotal},
dr = 10^-2;
 Table[
 dlTotal = 
  Plus @@ Table[(dFunc[#1, #2, r + dr] - dFunc[#1, #2, r]) & @@ 
     pair, {pair, points}];
 {r, dlTotal/(4 Pi r^2 dr)}
 , {r, 0.01, 1, dr}
 ]
];  
(*plot*)
ListPlot[sim]

Problem

1. Even for up to 10000 pts, the plot doesn't seem to tend to a constant value near 1/$\pi$ as it should as seen from the answer mentioned before

I've separately checked and rechecked the dFunc, which I presumed to be the only source of error, but found it to be correct; tweaked dr value; Checked the geometry -- but nothing seems wrong.

Answer 1

Update after comments: It turned out that the problem with OP's code was not dividing by the number of lines.

Here is a straightforward, alternative way to compute the same thing:

distance[Point[{pt1_, pt2_}]] := Norm[pt1 - pt2]
distance[_] := 0

massInsideSphere[line_, r_] := 
 distance@RegionIntersection[Sphere[{0, 0, 0}, r], line]

density[lines_, r_] := Module[{mass1, volume1, mass2, volume2},
  mass1 = Total[massInsideSphere[#, r] & /@ lines];
  volume1 = Volume[Ball[{0, 0, 0}, r]];
  mass2 = Total[massInsideSphere[#, r + 0.01] & /@ lines];
  volume2 = Volume[Ball[{0, 0, 0}, r + 0.01]];
  (mass2 - mass1)/((volume2 - volume1) Length[lines])
  ]

nrOfTrials = 10000;
lines = InfiniteLine /@ RandomPoint[Sphere[], {nrOfTrials, 2}];

estimates = Table[{r, density[lines, r]}, {r, 0.1, 0.9, 0.1}];
ListPlot[estimates]

This agrees approximately with the theoretical result $\rho = \frac{1}{\pi}$.

1./Pi

0.31831