Consider a light source which emits a Gaussian beam with total power $P_0$, initial width $d$, and divergence half-angle $\theta$ (in radians). The beam radius, measured from the waist of the beam, can be given as $$\widetilde{w}(z) = w_0 \sqrt{1+\left(\frac{z}{z_R}\right)^2},$$ where $z$ is the distance measured from the waist, $w_0$ is the waist radius $$w_0 = \frac{ \lambda}{\pi n \theta},$$ and $z_R$ is the Rayleigh range $$z_R = \frac{\pi w_0^2 n}{\lambda}.$$ Let $z_0$ be the distance of the origin of the beam measured from the waist, i.e., $$\widetilde{w}(z_0) = \frac{d}{2}.$$ Then the radius of the beam can be given as $$w(z) = \widetilde{w}(z+z_0),$$ where $z$ is the distance from the origin, not from the waist.
The intensity profile of the beam can be given as $$ \mathcal{I}(r, z) = \frac{2P_0}{\pi w_0^2} \left( \frac{w_0}{w(z)} \right)^2 \mathrm{exp}\left(-\frac{2r^2}{w(z)^2}\right),$$ where $r$ is the radial distance from the center axis of the beam, and $z$ is the distance from the origin of the beam.
At a distance of $z = R$ from the origin, there is a disc-shaped detector with diameter $D$. I am trying to calculate the measured intensity at the detector using the Gaussian beam model.
The power through a surface $S$ can be calculated as the surface integral $$\int_{S} \mathcal{I} dS,$$ so the power at the detector can be calculated as an integral in polar coordinates, $$ P:=\int_{r \le D/2} \mathcal{I}(r, R) dA = \int_{r = 0}^{D/2} \int_{\theta = 0}^{2\pi} \mathcal{I}(r, R)r d\theta dr.$$ This evaluates to $$P = P_0 \left(1- \mathrm{exp} \left( -\frac{D^2}{2 w(R)^2}\right)\right).$$
The measured intensity is the power through the detector divided by the area of the detector, i.e., $$\mathcal{I}_{measured} = \frac{P}{\pi (D/2)^2} = \frac{4P_0 \left(1- \mathrm{exp} \left( -\frac{D^2}{2 w(R)^2}\right)\right)}{ \pi D^2}.$$
However, I do not see how this is in line with the inverse square law, i.e., the measured intensity should be proportional to $R^{-2}$. With this model, I get results which are in completely different orders of magnitude. Alternatively, one would calculate the measured intensity simply as $$\frac{P_0}{2\pi (1-\cos{\theta}) R^2} \propto R^{-2}.$$