Providing truly ubiquitous connectivity requires development of low Earth orbit (LEO) satellite Internet, whose theoretical study is lagging behind network-specific simulations. In this paper, we derive analytical expressions for the downlink coverage probability and average data rate of a massive inclined LEO constellation in terms of total interference power's Laplace transform in the presence of fading and shadowing, ergo presenting a stochastic geometry based analysis. We assume the desired link to experience Nakagami m fading, which serves to represent different fading scenarios by varying integer m, while the interfering channels can follow any fading model without an effect on analytical tractability. To take into account the inherent nonuniform distribution of satellites across different latitudes, we model the LEO network as a nonhomogeneous Poisson point process with its intensity being a function of satellites' actual distribution in terms of constellation size, the altitude of the constellation, and the inclination of orbital planes. From the numerical results, we observe optimum points for both the constellation altitude and the number of orthogonal frequency channels; interestingly, the optimum user's latitude is greater than the inclination angle due to the presence of fewer interfering satellites. Overall, the presented study facilitates general stochastic evaluation and planning of satellite Internet constellations without specific orbital simulations or tracking data on satellites' exact positions in space.