#JCourtin 02/2016
#modified 12/2016

from matplotlib import pyplot as plt
from numpy import *
from cmath import *

"""
On note x la position sur le diaphragmme et X celle sur l'écran :

1 - on crée l'ensemble de émetteurs sur le diaphragmme (source de Fresnel)
2 - on crée les positions de calcul sur l'écran

3 - En chacun des points de l'écran on fait la somme de Fresnel des ondes émises par toutes les sources secondaires du diaphragmme. (c'est long !)


RQ : les tableaux ou diamètres d'ouverture du diaphragmme sont tous impairs
     car il y a une cellule centrale et un nombre symétrique de chaque coté.



a - Le TD se propose de tester le code sur le cas de la diffraction par une ouverture rectangulaire,
    Vous devez en particulier retourver les interfranges connus selon x et selon y.

b - Vous pouvez ensuite essayer d'autres résolutions :
    NX,NY pour le diaphragmme et MX MY pour l'écran

c - Enfin vous pouvez tester d'autres formes de diaphragmmes.



"""


############################   Diaphragmmes   ##################################


### Fonctions diaphragmme 


#ouverture cirulaire
def Disque(NX,NY, diam):   #Attention a est le diamètre en pixel => impair
    a=diam//2
    fente=zeros((NY,NX))
    XC=NX/2;YC=NY/2
    
    for i in range(NY):
        for j in range(NX):
            dist=((i-YC)**2+(j-XC)**2)**0.5
            if (dist < a):
                fente[i,j]=1.

    return fente


#ouverture fente infini
def fenteInfinie(NX,NY, ax):
    
    XC=NX//2 # NX = 2*XC + 1 impair
    
    fente=zeros((NY,NX))
    fente[:, NX//2 - ax//2 : NX//2 + ax//2]=1.
    
    return fente
    
    
    #fixer les bornes pour que l'ouverture soit au centre
    #NX,NY, ax, ay sont tous des entiers impairs.
    
def ouvertureRectangulaire(NX,NY, ax, ay):
    XC=NX//2 # NX = 2*XC + 1 impair
    YC=NY//2 # NY = 2*YC + 1 impair
    
    fente=zeros((NY,NX))
    fente[#bornes ???]=1.
    
    return fente




#Trous de Young
def doubleDisque(Nx, Ny, diam):
    a=diam//2
    if (2*a > Nx//3):
        return "PB"
    fente=zeros((NY,NX))
    
    XC1=NX//6; XC2=5*NX//6 #il y a deux centres ....
    
    YC=NY/2
    
    for XC in [XC1, XC2]:
        for i in range(NY):
            for j in range(NX):
                dist=((i-YC)**2 + (j-XC)**2)**0.5
                if (dist < a):
                    fente[i,j]=1.

    return fente

################################################################################


#                              *** Partie TD *** 

## 0 caracteristiques physiques
LOnde = 500.e-9 #lambda
k0=2*pi/LOnde   #vecteur onde

D = 10.   #metre  -> distance ecran
L = 0.01  #metre  -> echelle physique du diaphragme (!!! il y a 2 directions !!!)


## 1 Construction du diaphragmme

#taille matrice
NX=10+1     #nombre impair => 1 cellule au centre
NY=20+1     #Ne pas modifier ces valeurs dans un premier temps

dx=L/NX   #taille pixel (physique)
dy=L/NX   #même échelle


ax=11      #fente -> ouverture fente en pixel (!!! diametre ou largeur !!!)
ay=21

InterX=LOnde*D/(ax*dx)    #ODG en mètre de l'interfrange
InterY=LOnde*D/(ay*dy)  


#Choix du diaphragmme :

#ouverture fente infinie : a impair
#----------------------------------
#fente = fenteInfinie(NX, NY, ax)


#ouverture rectangulaire :
#-------------------------
fente=ouvertureRectangulaire(NX, NY, ax, ay)

#ouverture circulaire :
#----------------------
#fente = Disque(NX,NY, ax)


#trous de young :
#----------------
#fente = doubleDisque(NX, NX, ay)



#positions physiques x et y des cellules sur le diaphragme (faire un schéma)
x=[(-(NX)/2 + k)*dx for k in range(NX)]
y=[(-(NY)/2 + k)*dy for k in range(NY)]   #corriger les bornes !


#Objet numpy pour le diaphragme : matrice de tuple (y, x, fente(x,y))
#
#soit la nature de l'objet :   liste de liste (matrice) contenant en tout point
#                              un tuple (y, x, 0. ou 1.)

diaph = array([[(y[i], x[j], fente[i, j]) for j in range(NX)] for i in range(NY)])




## 2 Construction de l'écran

#Caracteristiques de l'image à produire :

MX=60      #toujours impair => 1 cellule au centre
MY=60      #Ne pas modifier ces valeurs dans un premier temps

#reduire à 30x30 si le temps de calcul est trop long.  [ t < 10 secondes ]

image =  ???  #A compléter: on veut initialiser la matrice avec des zeros

#La taille physique de l'écran est fixée en terme de nombre de fois l'interfrange :

TX=6*InterX   # position de calcul sur l'écran
TY=6*InterX   #même dimension

dX=TX/MX
dY=TY/MY

X=[#Compréhension de liste --> ??? for k in range(MX)]  #positions sur l'écran
Y=[#idem sur Y  : vérifier les bornes !!!

## 3 Calcul de la diffraction

for k in range(MY):
    for l in range(MX):
                
        #indiquer la bonne opération et expliquer la formule.
        #Quel est le rôle du if ?
        
        amplitude = #quelleOpération?([diaph[i,j,2]*exp(1j*k0/D*(Y[k]*diaph[i,j,0] + X[l]*diaph[i,j,1])) for i in range(NY) for j in range(NX) if diaph[i,j,2]!=0.])
    
        #Puis calcul de l'intensité résultante.
        image[k,l]=??  # A compléter


#Représentation de l'intersité sur l'écran
plt.close()
plt.matshow(image)
plt.show()

"""
#Intersité sur l'écran pour y=0 (coupe transverse)
plt.close()
plt.plot(X,image[0])   #Vérifier que vous retrouver les bons interfranges.
plt.show()

"""