# 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.

-
I am not sure I understand what you mean. A(x0,x0) is infinity, isn't it? The elementary wave is already normed in a way way that each wave has the power of 1 (integral of |A| of a circled centered at the source of the wave). That's why I wrote I expected a combines power of, as w=2 (but the slit radiates in two direction, so it reduces to 1). However, I do not understand why I should scale my result with that again. Could you expand your answer, per chance using the latex syntax? –  mcandril Sep 7 '12 at 12:56
Latex is a problem (for me). But I shall try to be clearer. In your formula for the amplitude, x0 is a point (secondary source) located in the slit aperture. You need to know what is the energy I0 of the incoming beam before it is radiated. Once you know it you scale the amplitude of each radiated wave by the square root of I0. So the problem is to determine I0. To do it you sum over x0 the squareroot of the modulus of A(x0,x0). –  Shaktyai Sep 7 '12 at 13:15
Re reading your post, I am not sure I understand your geometry. Is it linear or circular ? –  Shaktyai Sep 7 '12 at 13:18
A(x0,x0) is not infinity... –  Shaktyai Sep 7 '12 at 14:54
Sorry, I made a mistake defining my elementary wave - or better: Here I took it from the calculations and it differs from my normal version (where x and x0 are the vectors for the center and the current position). I should have checked that again :) Concerning the geometry: The first part is linear, while the code is circular. I should have made that more clear. I am working linear ATM, for the full intensity, I changed it to circular, to keep everything finite. I will check everything when I am back in my office on Monday, as I need a licence server for Mathematica. Thanks! –  mcandril Sep 8 '12 at 17:49

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}}$$.

This Green function has all sorts of implications for Huygen's construction. The main thing lies in the formation of wakes - wave amplitudes that linger behind the circular wavefront. These wakes do not exist in 1 and 3 dimensions. See Baker and Copson, Mathematical Theory of Huygens Principle, 1950, especially Chapter 1.

-
Can you maybe provide a more detailed reference to Baker and Copson (title, year)? –  Bernhard Dec 18 '12 at 10:07
Sorry, I thought I had provided it, apparently I forgot to paste what I had copied. –  rlgordonma Dec 18 '12 at 13:51