import matplotlib.pyplot as plt
import numpy as np
from scipy import integrate

def Fresnel2dreal(yp, xp, y, x, k, z):
    kernel = np.cos((k/(2*z))*((x-xp)**2+(y-yp)**2))
    return kernel

def Fresnel2dimag(yp, xp, y, x, k, z):
    kernel = np.sin((k/(2*z))*((x-xp)**2+(y-yp)**2))
    return kernel

c = 3e8
e0 = 8.85e-12

def plot1D(aperture, z, k, screen_range):

    y = 0

    xp1=yp1=-aperture/2
    xp2=yp2=aperture/2
    
    xs = np.linspace(-screen_range/2, screen_range/2, num=1000)
    intensities = []

    for x in xs:
        realpart, realerror = integrate.dblquad(Fresnel2dreal, xp1, xp2, yp1, yp2, args=(y, x, k, z), epsabs=1e-11, epsrel=1e-11)
        imagpart, imagerror = integrate.dblquad(Fresnel2dimag, xp1, xp2, yp1, yp2, args=(y, x, k, z), epsabs=1e-11, epsrel=1e-11)

        I = c*e0*(realpart**2+imagpart**2)
        intensities.append(I)

    return xs, intensities

ax = plt.axes()
xs, intensities = plot1D(2e-4, 0.005, 8.377e6, 0.002)
ax.plot(xs, intensities)
# xs, intensities = plot1D(2e-5, 0.05, 8.377e6, 0.015)
# ax.plot(xs, intensities)
plt.show()