How to set up the right model for solving diffraction numerically? The diffraction of waves in free space is governed by the sourceless Helmholtz equation
$$\nabla^2u +k^2u=0$$
There is a well stablished solution called Kirchhoff´s diffraction integral (KDI)  and it´s deduction is straightforward from greens theorem, the solution is
$$ u(r) = \int_{S} \left [ u \space \partial_n \left(\frac{e^{iks}}{s} \right)  - \left(\frac{e^{iks}}{s} \right) \space \partial_n u\space \right] dS$$
Where $S$ is a close surface containing the aperture and an imaginary boundary at infinity that fulfills Sommerfeld radiation condition and $s$ is the distance between any point in the surface and the point $r$ inside the enclosed volume.
Now I am trying to solve the Helmholtz equation numerically using the Finite Element Method (FEM). For this, I am assuming that a scalar plane wave propagating in the z-direction is illuminating a 1D screen with an aperture at the center. To set up the model, I am assuming the following Neumann conditions for the dependent variable u(x,z):
Absorbing boundary condition (ABC) to model an open system (boundaries at infinity) and ABC for the diffraction obstacle too, $n\cdot \nabla u=-ik_0 u$
radiation condition for the aperture to model illumination by a plane wave, $n \cdot \nabla u= ik_0 ((1+n\cdot e)-u)$
where $n$ is the boundary normal and $e$ is the incoming wave propagation direction. The solution obtained with this model is similar to the one using the KDI, but there is a mismatch in the scale and phase of the solution

I have tried different boundary conditions, for example, Sommerfeld's and Kirchhoff's, the ones used to deduce the KDI, but the results are similar. The solutions differ in a factor of $\sqrt{s}$ in the 2D case.
I found that the numerical and theoretical solutions match in the 3D case, i.e. u(x,y,z) and the aperture being a surface instead of a line as in the 2D case.
The problem remains in how to reduce the KDI to 2D instead of 3D, normally when one consider a 2D diffraction problem, one assumes that, for example, the aperture extend to infinity in the $y$ direction, therefore it is not necessary to take the $y$ coordinate into account. I think that maybe the $\sqrt{s}$ factor comes from the reduction to 2D of the KDI. I tried doing the integration from $[-\infty,\infty]$ numerically, but it didn't work.
I will appreciate if someone knows what is wrong with the model or how to reduce de KDI to 2D rigorously.
Note: The boundary conditions were taken from Jian Ming's "The Finite Element Method in Electromagnetism".
 A: Short instruction. In 2D case one needs instead of $\frac{e^{ikr}}{r}$ to put a Bessel function of the second kind. The 2D version of KDI is the following:
$$
u(r)=\frac{1}{4} \left[ \int_S u \partial_n Y  - Y \partial_n u  \right]\quad(1)
$$
where
$$
Y=-Y_0(k r),
$$
where $Y_0(x)$ is a Bessel function of the second kind
$$
\pi Y_0(x)=\int_0^\pi \sin(x \sin{\theta})d \theta - \int_0^\infty e^{-x \sinh(t)}dt.
$$
For large $x$ one can use an approximation
$$
 Y_0(x)\approx-i\frac{1}{\sqrt{2\pi x}}e^{i(x-\frac{\pi}{4})}.
$$
You can notice that your phase shift is about $\pi/4$ everywhere, so this formula is quite likely applicable to your case.
Some background. For any two functions, $u$ and $u_1$, satisfying the Helmholtz equation, the following integral can be shown to be zero
$$
\int_S u \partial_n u_1  - u_1 \partial_n u = 0.\quad (2)
$$
Function $u_1=Y_0(k r)$ satisfies Helmholtz equation everywhere except at $r=0$. Picking an area that includes $r=0$, and cutting out a tiny disk out around $r=0$ ensures that the integral is zero over such area. The integral over the perimeter of the tiny disk gives us $u(r=0)$ (with some coefficient), while the integral over the big outside perimeter gives us the right side of the Eq. (1).
A: Another good method is the Fourier approach. First decompose your aperture function in terms of a set of plane waves propagating in different directions, then propagate each plane wave to the observation plane, then add them up again. Thus two Fourier transforms with a simple phase factor introduced in between them. Sorry I don't have time now to write out the details but I hope this might be helpful anyway.
