(* Code for calculating the probability distribution function of the disorientation angles of two cubes, its mean, standard deviation, and median, based on J. K. Mackenzie, Second Paper on Statistics Associated with the Random Disorientation of Cubes, Biometrika, Vol. 45, No. 1-2 (1958), pp. 229-240. *)

d1 = Pi/4; 
d2 = Pi/3;
d3 = 2 * ArcTan[Sqrt[2]*(Sqrt[2]-1)];
d4 = 2 * ArcTan[(Sqrt[2]-1) * Sqrt[5 - 2*Sqrt[2]]];
p1[d_] := (24/Pi) * (1-Cos[d]);
p2[d_] := (24/Pi) * (3 * (Sqrt[2]-1) * Sin[d] - 2 * (1-Cos[d]));
p3[d_] := (24/Pi) * ((3 * (Sqrt[2]-1) + 4/Sqrt[3]) * Sin[d] - 6 * (1-Cos[d]));
X[d_] := (Sqrt[2]-1)/Sqrt[1-(Sqrt[2]-1)^2 * Cot[d/2]^2];
Y[d_] := (Sqrt[2]-1)^2/Sqrt[3 - Cot[d/2]^2];
p4[d_] := (24/Pi) * ((3 * (Sqrt[2]-1) + 4/Sqrt[3]) * Sin[d] - 6*(1-Cos[d])) - 
          (288 * Sin[d]/Pi^2) * (2*(Sqrt[2]-1) * ArcCos[X[d]*Cot[d/2]]+(1/Sqrt[3])*ArcCos[Y[d]*Cot[d/2]])+
          (288 * (1-Cos[d])/Pi^2) * (2 * ArcCos[X[d] * (Sqrt[2]+1)/Sqrt[2]] + ArcCos[Y[d] * (Sqrt[2] + 1)/Sqrt[2]]);
p[d_] := Which[d <= d1, p1[d], d1 < d <= d2, p2[d], d2 < d <= d3, p3[d], d <= d4, p4[d]];

RealDigits[NIntegrate[d * p[d], {d, 0, d4}, WorkingPrecision -> 110], 10, 100][[1]] (* mean *)

RealDigits[Sqrt[NIntegrate[d^2*p[d], {d, 0, d4}, WorkingPrecision -> 110] - NIntegrate[d*p[d], {d, 0, d4}, WorkingPrecision -> 110]^2], 10, 100][[1]] (* standard deviation *)

RealDigits[x /. FindRoot[Integrate[p[d], {d, 0, x}] == 1/2, {x, 0.5}, WorkingPrecision -> 110], 10, 100][[1]] (* median *)