I am trying to simulate (MonteCarlo integrate) the following scenario, see sketch below (as prequel of a somewhat bigger simulation).
Assume a small (point) L'Ambertian emitter, i.e. the radiant intensity is distributed as cosine $I=I_0 cos(\theta)$. In a distance $d$ and coaxially, there is a circular aperture of width $a$. How does the total power $\Phi$ transferred from the source to the aperture decrease with distance? We know, that at large distance one should obtain an inverse square law $\Phi \propto 1/d^2$.
This is how I tried to code:
i) Generate N samples (i.e. rays, actually just their emission angles) with distribution $cos(\theta)$.
ii) Count only those samples/rays, that have a value $\theta < \theta_{max} = \arctan(\frac{a}{2d})$
iii) Do that for all distances and plot the counted samples from ii) as function of distance, see image below.
The problem; Only when I square the count values I obtain an inverse square dependence. Actually, it even reproduces a textbook expression for small $d$.
Question: It seems that my sketch is only 2Dimensional, while the inv.-square law refers to 3D space. But, I cannot rationalize why I suddenly have to square the counts. Do I have to setup the sampling differently to get the right dependence directly? I would appreciate if someone could point out (formally) where this squaring comes from.
Thank you!