#!/usr/bin/env python
# coding: utf-8

# ## Fresnelova difrakcija na dvostrukoj pukotini
# #### Kod je pisan u Python programskom jeziku koristeći Jupyter notebook

# In[1]:


import numpy as np
import matplotlib.pyplot as plt
import math

from scipy.fftpack import fft      #paket za FFT
from scipy.fftpack import ifftshift,  fftshift

from mpl_toolkits.mplot3d import Axes3D   #paket za 3d plot
plt.style.use('dark_background')          # izbor boja za plot
plt.rcParams['image.cmap'] = 'plasma'


# In[2]:


def Fresnel_input(lambd,w,z):
    N=1048   #rezolucija
    L=5.    # širina zaslona na kojoj se javlja difrakciona slika

    lamb=lambd*10**(-6) #(mm)
#    w  #širina otvora (mm)
#    z    #udaljenost otvor - zaslon (mm)
    k=2*np.pi/lamb  # talasni broj 

    delta_src=L/N     #prostorna frekvencija

    #definisanje ravni ovora
    xv=np.linspace(-N/2,N/2-1,N)
    yv=np.linspace(-N/2,N/2-1,N)
    [x1,y1]=np.meshgrid(xv,yv)
    x1=x1*delta_src  #(mm)
    y1=y1*delta_src  # (mm) 
 
    #definisati dvostruku pukotinu
    X=np.where(np.abs(x1)<w/4,0,0)
    Y=np.where(np.abs(y1)<10*w,1,0)
    X[np.logical_and(x1<-0.05*w, x1>-0.15*w)]=1.0
    X[np.logical_and(x1<0.15*w, x1>0.05*w)]=1.0

    u1=X #elekticno polje na pukotini
    ap1=u1
    
# Fourierova transformacija  izračunavanje integrala 
    g=u1*np.exp(1j*np.pi/(lamb*z)*x1**2)
    G=fftshift(np.fft.fft(g))
    h=np.exp(1j*k*z)/1j*lamb*z*np.exp(1j*np.pi/(lamb*z)*x1**2)
    u2=h*G
    identity=np.ones((N,N))
    u2=u2*identity
    
     # Intenzitet zračenja (normiran na jedinicu)
       
    I=1./(lamb*z)**2*(np.abs(u2))**2
    I_1d=np.zeros(N, dtype='complex128')
    for a in range (1,N):
        I_1d[a]=I[0][a]
        
    
    I_1d=I_1d/I_1d.max()
        
    return(u2,I_1d,y1,ap1) #OUTPUT: električno polje na zaslonu, intenzitet, koordinate ravni otvora i polje na otvoru


# In[3]:


lamb=632.
w=1.
z=20.

u1,I1,x1,ap1=Fresnel_input(lamb,2*w,30.)
u2,I2,x2,ap2=Fresnel_input(lamb,w/4,30.)
u3,I3,x3,ap3=Fresnel_input(lamb,0.6*w,200)

I1=I1/I1.max()
I2=I2/I2.max()
I3=I3/I3.max()


# In[5]:


# PLOTANJE relativne veličine otvora, difrakcione slike i intenziteta zračenja
fig2=plt.figure(figsize=(11,10)) 

ax1=fig2.add_subplot(4,3,1)
plt.imshow(np.abs(ap1))
plt.xlim(324,724)
plt.ylim(460,588)
plt.axis('off')

ax1=fig2.add_subplot(4,3,2)
plt.imshow(np.abs(ap2))
plt.xlim(324,724)
plt.ylim(460,588)
plt.axis('off')

ax1=fig2.add_subplot(4,3,3)
plt.imshow(np.abs(ap3))
plt.xlim(324,724)
plt.ylim(460,588)
plt.axis('off')


ax1=fig2.add_subplot(4,3,4)
plt.title('(A) z = 30 mm, w = 2 mm')
plt.imshow(np.abs(u1))
plt.xlim(324,724)
plt.ylim(460,588)
plt.axis('off')
plt.xlabel('(A)')

ax2=fig2.add_subplot(4,3,5)
plt.title('(B) z = 30 mm, w = 0.25 mm')
plt.imshow(np.abs(u2))
plt.xlim(324,724)
plt.ylim(460,588)
plt.xlabel('(B)')
plt.axis('off')

ax3=fig2.add_subplot(4,3,6)
plt.title('(C) z = 200 mm, w = 0.6 mm')
plt.imshow(np.abs(u3))
plt.xlim(324,724)
plt.ylim(460,588)
plt.xlabel('(C)')
plt.axis('off')

ax4=fig2.add_subplot(2,1,2)
plt.title('Fresnelova difrakcija - dvostruka pukotina')

ax4_1=ax4.plot(x1,np.abs(I1),'r-')
ax4_2=ax4.plot(x2,np.abs(I2),'b-')
ax4_3=ax4.plot(x2,np.abs(I3),'y:')

ax4.text(-0.6,0.5,'(A)', color='red',fontsize=16)
ax4.text(0.4,0.5,'(B)', color='blue',fontsize=16)
ax4.text(0.05,0.,'(C)', color='yellow',fontsize=16)

plt.xlabel('X(mm)')
plt.ylabel('Relativni intenzitet')
fig2.savefig('Fresnel_2slit.png')
plt.show()


# In[ ]: