¿Cómo simular una dependencia de la ley del cuadrado inverso del flujo radiante de los rayos muestreados?

Estoy tratando de simular (integrar MonteCarlo) el siguiente escenario, vea el boceto a continuación (como precuela de una simulación algo más grande).

Supongamos un emisor L'Ambertiano pequeño (puntual), es decir, la intensidad radiante se distribuye como coseno$I=I_0 cos(\theta)$. en la distancia$d$y coaxialmente, hay una abertura circular de ancho$a$. ¿Cómo es la potencia total?$\Phi$transferida de la fuente a la apertura disminuye con la distancia? Sabemos que a gran distancia se debe obtener una ley del inverso del cuadrado$\Phi \propto 1/d^2$.

Así es como traté de codificar:

i) Generar N muestras (es decir, rayos, en realidad solo sus ángulos de emisión) con distribución$cos(\theta)$.

ii) Cuente solo aquellas muestras/rayos que tengan un valor$\theta < \theta_{max} = \arctan(\frac{a}{2d})$

iii) Haga eso para todas las distancias y trace las muestras contadas de ii) en función de la distancia, vea la imagen a continuación.

El problema; Solo cuando elevo al cuadrado los valores de conteo obtengo una dependencia cuadrada inversa. En realidad, incluso reproduce una expresión de libro de texto para pequeños$d$.

Pregunta : Parece que mi boceto es solo bidimensional, mientras que la ley del cuadrado inverso se refiere al espacio 3D. Pero, no puedo racionalizar por qué de repente tengo que cuadrar las cuentas. ¿Tengo que configurar el muestreo de manera diferente para obtener la dependencia correcta directamente? Agradecería si alguien pudiera señalar (formalmente) de dónde viene esta cuadratura.

¡Gracias!