# Huygens Fresnel Diffraction integral using dblquad in python

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.

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

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