# Conservation of energy with Huygens-Fresnel principle

I am currently experimenting with Huygens-Fresnel principle. I am trying to simulate the propagation of a beam, emerging from an aperture slot of width w. I assume the slot to be long and therefore work in two dimension.

My elementary wave is given by $$A_E(x,x_0)=\frac{\exp\left(2\pi i/\lambda 2\sqrt{(x-x_0)^2+a^2} \right)}{\left(2\pi\sqrt{(x-x_0)^2+a^2}\right)^{(1/2)}}$$ where a is the distance to the screen, x is the position on the screen and {x0,0} is the center of the elementary wave. This wave should carry an intensity of 1.

According to the Huygens-Fresnel principle, I should be able to calculate the amplitude in any place {x,a} behind the slot via $$A(x)=\int_{-w/2}^{w/2} A(x,x_0)\text{d}x_0$$ $$\approx\sum_{n=0}^{w/\Delta w} A(x,n*\Delta w) \Delta w$$ This seems to work quite well (result does not change notably for smaller step width) for $\Delta w<\lambda/5$.

I am using in Mathematica in 2 dimensions. This should numerically calculate the image on a 180° screen in a distance of 10 from the center of the aperture:

r = 10
L = .1
SlotRange = {-1, 1}
SlotStepW = L/10
ImageRange = {-Pi/2,Pi/2}
ImageStepW = (ImageRange[[2]] - ImageRange[[1]])/1000
beam = Table[
Sum[
N[
Exp[
I (2 Pi/L Sqrt[(r Sin[Phi] -
x0)^2 + (r Cos[Phi])^2])]/(2 Pi Sqrt[(r Sin[Phi] -
x0)^2 + (r Cos[Phi])^2])^(1/2)
]*SlotStepW,
{x0, SlotRange[[1]], SlotRange[[2]], SlotStepW}
],
{Phi, ImageRange[[1]], ImageRange[[2]], ImageStepW}
];
ListPlot[Abs[#]^2 & /@ beam, DataRange -> ImageRange, PlotRange -> All]
energy = (Plus @@ (Abs[#]^2 & /@ beam))*ImageStepW*r


Basically, this works quite well, as said. For smaller distances to the screen, the intensity distribution has a (oscillating) flat top, and in the distance it becomes roughly Gauss shaped. The full energy is always around 0.032. I would expect, however, something like 1 (elemental waves with intensity 1 over a range of length 2, gives two, but I am just getting half the intensity, as the situation is symmetrical around the aperture slot, so there is a backwards beam with another intensity 1).

Probably I am just doing some small mistake, but I cannot see it. However, I need proper amplitudes, as in the next step I would like to start new waves from the screen, continuing the beam from there. I initially thought I could just multiply the amplitudes on the screen to the elemental wave above, to get an intensity and phase adjusted new wave, but that doesn't seem to be the case ...

-

I have found what is missing:

The amplitude of the diffracted wave you use is missing the Fresnel-kirchhoff $K(x)$ factor. (Born and wolf page 413 and 423) You should use instead:

$$-\frac{i}{2\lambda} \frac{\exp(iks)}{s}(1+\cos(X))$$

instead of $$\frac{\exp(iks)}{s}$$

With $$X= \pi-(k_0,s)$$

-
 Oh, great. After I did not really get forward with the Jackson Electrodynamics, that book would have been my next try. It is already lying on my desk in my office for Monday :) – mcandril Sep 8 '12 at 17:33

You need to scale the amplitude of the incoming wave by the square root of the incoming wave energy= sqrt(sum (modulus(A(x0,x0), x0=-w/2..w/2)). Unless you do it, the incoming energy is not 1, that is why the diffracted energy is not one as well.

Do the normalisation by using the same approximate formula that you use to compute the integral.

-
The problem is that you are using the wrong Green function for the 2D Helmholtz operator. The Green function of this operator that satisfies conditions at infinity (outgoing waves only) is a Hankel function $H_0(k x)$. Note that the asymptotic behavior of this function is
$$H_0(k x) \approx \frac{\exp{(i k x)}}{\sqrt{k x}}$$.