Why do I obtain the wrong total power when integrating this Fraunhofer diffraction irradiance pattern for a circular aperture? I am looking for some help understanding something. I am working through my optics textbook (Hecht, 4th Edition), and am following the example of the far-field Fraunhofer diffraction pattern caused by a plane wave passing through a circular aperture. The situation is illustrated in the schemtaic below:

The book presents the general expression for the "optical disturbance" at a distance $R$ when using an arbitrary aperture, which is given as
$$
\tilde{E} = \frac{\mathcal{E}_A\; e^{i(\omega t - kR)}}{R}\iint e^{ik(ux+vy)/R} dS \tag{1}
$$
where the integration is over the aperture in the source plane. The book then goes on to calculate for the specific geometry of a circular aperture to arrive at the following result:
$$
\tilde{E}(\rho,t) = \frac{\mathcal{E}_A\; e^{i(\omega t - kR)}}{R} 2\pi a^2\frac{R}{ka\rho}J_1\Bigg(\frac{ka\rho}{R}\Bigg), \tag{2}
$$
where $J_1$ is the first order Bessel function of the first kind, and $\rho$ is the radial coordinate in the image plane (because the pattern is circularly symmetry). Finally, the book calculates the irradiance in the image plane using
$$
I(\rho) = \frac{1}{2}\tilde{E} \tilde{E}^* = \frac{2\mathcal{E}_A^2A^2}{R^2} \Bigg[ \frac{J_1(ka\rho/R)}{ka\rho/R} \Bigg]^2. \tag{3}
$$
My question is really regarding the scaling factors, because I am not sure what $\mathcal{E}_A$ really is. The book describes this as the "source strength per unit area", which is a bit vague to me.
In my schematic, I want the total power through the aperture to be the input to the problem, and am prescribing this to be $P=1$mW, so the irradiance at the aperture plane I know is constant at
$$
I_{in}=\frac{P}{\pi a^2}=\frac{1\textrm{mW}}{\pi (1\textrm{mm})^2}=318\;\textrm{W}/\textrm{m}^2. \tag{4}
$$
How do I use this number to calculate $\mathcal{E}_A$? I should also be able to calculate the electric field amplitude at the aperture using the relation (assuming constant phase across aperture)
$$
I_{in} = \frac{1}{2}c\epsilon_0|E|^2 \rightarrow E = \sqrt{\frac{2I_{in}}{c\epsilon_0}} = 490\;\textrm{V/m}, \tag{5}
$$
but it seems that this is not equal to $\mathcal{E}_A$?
Eventually, I would like to be able to check that when numerically integrating Eq. (3) over the entire $(x,y)$ plane, I obtain the correct 1mW of power (assuming no losses). To test this, I have created the following Matlab script, but it can be seen that when I integrate the final diffraction pattern I don't recover my original 1 mW of power. What am I misunderstanding?
Thank you!
P = 1e-3;              % Power through aperture [W]
a = 1e-3;              % Aperture radius [m]
I_in = P/(pi*a^2);     % Intensity at aperture [W/m^2]

c = 3e8;               % speed of light
e0 = 8.85e-12;         % vacuum permittivity
source_field = sqrt(2*I_in/(c*e0)); % 490 V/m

epsilonA = 1/lambda*sqrt(P/(2*A))*source_field;    % "Source strength per unit area" in Hecht language

A = pi*a^2;      % Aperture area [m^2]
R = 1000;        % Distance to screen [m]
lambda = 780e-9; % Wavelength of light [m]
k = 2*pi/lambda; % Wavenumber of light


% Image plane coordinates
x = linspace(-5,5,1000);y = x; % [m]
dx = x(1) - x(2);
dy = y(1) - y(2);
[X, Y] = meshgrid(x, y);
rho = sqrt(X.^2 + Y.^2); % radial coordinate

% Calculate diffraction pattern (Hecht, Optics)
I = 2 * epsilonA^2 * A^2 / R^2 * (besselj(1,k*a*rho/R)./(k*a*rho/R)).^2; 

% Check total power in diffraction pattern
power_image_plane = sum(sum(I))*dx*dy; % gives wrong value of 2.3e-13 Watts


 A: I will follow Born&Wolf [1] Sections 8.3.3 and 8.5.2 according to which, see eq (10)
$$U(\mathcal{P})= 2\pi\mathcal{C} \int_0^a J_0(k\rho w)\rho d\rho \tag{8.5.10}\label{10}$$where $$\mathcal{C}=\frac{1}{\lambda R}\sqrt{\frac{\mathcal{\dot E}}{\mathcal{D}}} \tag{8.3.42}\label{42}$$ and $\mathcal{\dot E}$ is defined as the total integrated power (energy) through (in) the circular aperture whose area is $\mathcal{D}=\pi a^2$. In the center of the screen the intensity follows from the Fourier transform relationship, see Parseval's theorem, (8.3.40) and is equal to
$$I_0 = \frac{\mathcal {\dot E} \mathcal {D}}{\lambda^2 R^2}=\mathcal{C}^2 \mathcal D^2 \tag{8.3.44}\label{44}$$
Now the diffracted power (intensity) over the screen at a point $\mathcal{P}$ that is at an angle $\gamma$ from the center line and a distance $R$ from the aperture is Airy's formula eq 8.5.14:
$$I(\mathcal{P}) = |U(\mathcal{P})|^2=I_0\left|\frac{2J_1(kaw)}{kaw} \right|^2 \tag{8.5.14}\label{14}.$$
Here $\gamma \approx w = \frac{\ell}{R}$
Now let us integrate $\eqref{14}$ over the whole screen. The area element on the *screen* in polar coordinates $\ell,\alpha$ is $\ell d\ell d \alpha$ but because of the circular symmetry of $I(\mathcal P)$ the integration over the angle is $2\pi$ and we get the total power deposited on the screen as
$$ K =2\pi\int_0^{\alpha_{max}} I_0\left|\frac{2J_1(kaw)}{kaw} \right|^2 \ell d\ell $$
Now let $x=kaw = ka \frac{\ell}{R}$ or $ \ell = \frac{Rx}{ka}$ then
$$ K =2\pi \frac{I_0 R^2}{(ka)^2}\int_0^{x_{max}} \left|\frac{2J_1(x)}{x} \right|^2x dx \\= 4\pi \frac{I_0 R^2}{(ka)^2} 2\int_0^{x_{max}}  \frac{J_1^2(x)} {x} dx $$
The integral can be expressed explicitly
$$2\int_0^{x_0}  \frac{J_1(x)^2} {x} dx  = 1- J_0^2(x_0)-J_1^2(x_0) \tag{8.5.12}\label{12}$$
and for $x_0 \to \infty$ the RHS of $\eqref{12}$ goes to $1$, and since $x=ka\ell/R$ this limit is the same as $k\to \infty$; therefore with $\eqref{44}$ we get
$$\begin{align}
K\approx 4\pi \frac{I_0 R^2}{(ka)^2} \\
= 4\pi\frac{\mathcal {\dot E} \mathcal {D} R^2}{\lambda^2 R^2 (ka)^2}
= 4\pi\frac{\mathcal {\dot E} \mathcal {D}}{4\pi^2 a^2} \\= \frac{\mathcal {\dot E} \mathcal {D}}{\pi a^2}= \mathcal {\dot E}
\end{align}$$
That is $K\approx \mathcal{\dot E}$, which shows that the diffracted power is the same as the one passing through the aperture within the paraxial (far-field) approximation.

The dimensions are consistent. If we measure the intensity $I_0$ in units of $\rm{W} \rm{m}^{-2}$, which is dissipated power per unit area, then from $\eqref{44}$ $[\mathcal C] = \sqrt { \rm{W} m^{-2} } \rm{m}^{-2} = \sqrt{\rm{W}} \rm{m}^{-3}$. Similarly, using $\eqref{42}$ we get $[\mathcal {\dot E}] = [\lambda^2 R^2 I_0 /\mathcal {D}] = \rm{m}^2 \rm{m}^2 \rm{W} \rm{m}^{-2} \rm{m}^{-2}= \rm{W} $ that represents the total power incident on the aperture, as it should be.
[1] https://archive.org/details/PrinciplesOfOptics/page/n433/mode/2up
