You pointed out one error in your own answer, here's another. I haven't figured out the adjustments, but presumably they cancel.

You've done

∫₀¹ ∫₀¹⁻ᶜ (2π r)³ dr dc

However,

- the `dr` integral is assuming that the radii are uniformly likely in [0, r]

- the `dc` integral is assuming that the centers are uniformly likely in [0, c]

You need to wait these integrals by the conditional probability distributions.