#!/usr/bin/env python3
# -*- coding: utf-8 -*-

"""
Descriptif du fichier
"""

# Libraries import
import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import nquad,dblquad
import scipy

# Constants definition
D = 6e-5
lamb = 6e-7
k = 2*np.pi/lamb
a = 1e-5
d = a*10000*D

F = a**2/(D*lamb)

print(F)

# Function definitions



# Main program
X = np.linspace(-d,d,100)
T = lambda x,y : (1+0.5*np.sin(1*x/lamb))*(x>0)




def E(x,y,x1,y1):
    r = np.sqrt(D**2+(x-x1)**2+(y-y1)**2)
    return (np.exp(1j*r*k))*T(x,y)/lamb

def rE(x,y,x1,y1):
    r = np.sqrt(D**2+(x-x1)**2+(y-y1)**2)
    return(np.real(E(x,y,x1,y1)))

def iE(x,y,x1,y1):
    return(np.imag(E(x,y,x1,y1)))

def vE(x,y,x1,y1):
    return np.array([rE(x,y,x1,y1),iE(x,y,x1,y1)])



def I(x1,y1) :
    rarea = dblquad(lambda x, y: rE(x,y,x1,y1), -a,a,lambda x : -a,lambda x : a, epsrel = 1e6)
    iarea = dblquad(lambda x, y: iE(x,y,x1,y1),-a,a, lambda x : -a,lambda x : a, epsrel = 1e6)
    return np.abs(rarea[0] + iarea[0]*1j)**2



    

Y = [[I(x,0),print("A")] for x in X]

plt.plot(X/a,Y), plt.show()