Huygens Fresnel Diffraction integral using dblquad in python

Question

I am attempting to create a python function to assist in calculating the following numerical integration of the Huygens Fresnel integral in the form of def fresnel(x_,y_,x,y,z) where x_ and y_ are $x'$ and $y'$ coordinates.

$$E(x,y,z)=\frac{1}{i\lambda}\iint_{-\infty}^{+\infty}E(x',y',0)\frac{ze^{ikr}}{r}dx'dy'$$

where

$$r=\sqrt{(x-x')^2+(y-y')^2+z^2}$$

$k=\frac{2\pi}{\lambda}$is the wavenumber, and $\lambda$ is the wavelength of the lightsource, and it is assumed that light field is uniform, $E(x',y',0)=1$

The few examples I have seen online use

scipy.integrate.dblquad(func, a, b, gfun, hfun, 
args=(), epsabs=1.49e-08, epsrel=1.49e-08)

to perform the double integration, however, I've been struggling to implement the gfun and hfun callables of the dblquad function as I'm still quite new to this.

Answer 1

The limits of integration on dblquad can be functions (callable type) or floats. In your case you need the limits to be infinite.

Following is an example integrating a Gaussian function over $\mathbb{R}^2$, the result should be $\pi$.

from numpy import exp, pi, inf
from scipy.integrate import dblquad


def gaussian(x, y):
    return exp(-x**2-y**2)


inte = dblquad(gaussian, -inf, inf, -inf, inf)

print(inte)
print(pi)

And the result is

(3.141592653589777, 2.5173086244657047e-08)
3.141592653589793