clear all ;close all;clc;
% wavelength
lambda=0.000632;
%the distance from grating to the image plane,in general it's focus length
f=1000;
%modulate factor
m=1;
f0=2;% grating period
L=5;% grating size
theta=0;%initialize phase
[x,y]=meshgrid(linspace(-2.1,2.1,1500));
%xx=1.6*lambda*f/1;
%[x,y]=meshgrid(-xx:xx/30:xx,-xx:xx/30:xx);
a=L/(lambda*f);
i1=(L.^2/(2*lambda*f)).^2;
i2=sinc(a.*y).^2;
i3=sinc(a.*x).^2;
i4=sinc(a.*(x+f0*lambda*f)).^2.*exp(-i*theta).^2;
i5=sinc(a.*(x-f0*lambda*f)).^2.*exp(i*theta).^2;
I=(i3+0.25.*i4+0.25.*i5);

figure(1);
plot(x/(lambda*f),I/max(max(I)),'r');
figure(2);
I=I.*i1.*i2;
NCLevels=255;
Ir=(I/1)*NCLevels;
colormap(gray(NCLevels));
image(Ir);
%mtf=fftshift(fft2(I));figure;plot(x,abs(mtf))