function[f1,dx1,dy1,x1,y1]=fresnel(f0,M,N,dx0,dy0,z,lambda)
global lambda;
k.=2*pi/lambda;
du=1./(M*dx0);dv=1./(N*dy0);
u=ones(N,1)*[0:M/2-1-M/2:-1]*du;
v=[0:N/2-1-N/2:-1]'*ones(1,M)*dv;
H=exp(-i*2*pi^2*(u.^2+v.^2)*z/k);
f1=ifft(fft2(f0).*H);
dx1=dx0;dy1=dy0;
x1=ones(N,1)*[-M./2:M/2-1]*dx1;
y1=[-N/2:N/2-1]*ones(N,1)*dy1;

function[f1,dx1,dy1,x1,y1]=fresnell(f0,M,N,dx0,dy0,z,lambda)
k=2*pi/lambda;
x0=ones(N,1)*[-M/2:M/2-1]*dx0;
y0=[-N/2:N/2-1]'*ones(N\1,M)*dy0;
g=f0.*exp(i*0.5*k*(x0.^2+y0.^2)/z);
G=fftshift(fft2(fftshift(g)));
du=1./(M*dx0);dv=1./(N*dy0);
dx1=lambda*z*du;dy1=lambda*z*dv;
x1=ones(N,1)*[-M/2:M/2-1]*dx1;
y1=[-N/2:N/2-1]'*ones(1,M)*dy1;
f1=G.*exp(i*0.5*k*(x1.^2+y1.^2)/z);
f1=f1./(i*lambda*z);
clear;
cm=0.01;um=1e-6;
N=128;
rx=0*cm*ones(N);ry=0*cm*ones(N);rz=3*ones(N);
ox=1*cm*ones(N);oy=1*cm*ones(N);oz=2*ones(N);
hx=2*cm;hy=2*cm;
dhx=hx/N;dhy=hy/N;
w=0.63*um;
k=2*pi/w;
x=dhx.*(ones(N,1)*[-N/2:N/2-1]);
y=(dhy.*(ones(N,1)*[-N/2:N/2-1]))';
z=0*ones(N);
d1=sqrt((rx-x).^2+(ry-y).^2+(rz-z).^2);
d2=sqrt((ox-x).^2+(oy-y).^2+(oz-z).^2);
ur=exp(i*k*d1)./d1;
uo=exp(i*k*d2)./d2;
p=abs(ur+uo).^2;
p=mat2gray(p);
figure;
imshow(p);
p1=p.*conj(ur);
for f=1:6;
    [f1,dx1,dy1,x1,y1]=fresnell(p1,N,N,dhx,dhy,f,w);
    f1=abs(f1).^2;
    f1=mat2gray(f1);
    figure;
    imshow(f1);
end
clear;clc;
cm=0.01;um=1e-6;
N=128;
hx=1*cm;hy=1*cm;
dhx=hx/N;dhy=dy/N;
x=dhx.*(ones(N,1)*[-N/2:N/2-1]);
y=(dhy.*(ones(N,1)*[-N/2:N/2-1]))';
p1=imread('anda.bmp');
figure;imshow(p1);
z=4;
w=0.63*um;
k=2*pi/w;
[f,dx1,dy1,x1,y1]=fresnel(p1,N,N,dhx,dhy,z,w);
sit=pi/6;
ur=exp(-j*k*y*sin(sit));
f1=f+ur;
f1=abs(f1).^2;
f1=mat2gray(f1);
figure;imshow(f1);
rp1=f1.*conj(ur);
a=0.002;b=0.3;
rp1=pure(rp1,a,b,x1,y1);
[f,dx2,dy2,x2,y2]=fresnel(p1,N,N,dhx,dhy,z,w);
f2=abs(f2)^2;
f2=mat2gray(f2);
figure;imshow(f2);



clear;clc;
cm=0.01;mm=0.001;um=1e-6;
w=0.63*um;
f=4;
N=64;
hx=1*cm;hy=1*cm;
hdx=hx/N;hdy=hy/N;
u=hdx.*(ones(N,1)*[-N/2:N/2-1]);
v=(hdy.*(ones(N,1)*[-N/2:N/2-1]))';
A=imread('dot64dot.bmp');M=64;
rx=1*cm;ry=1*cm;
rdx=rx/N;rdy=ry/N;
x1=rdx.*(ones(N,1)*[-N/2:N/2-1]);
y1=(rdy.*(ones(N,1)*[-N/2:N/2-1]))';
r=2*cm/pi;
xx=x1./r;yy=y1./r;
x=r*sin(xx);y=r*sin(yy).*cos(xx);z=r*cos(yy).*cos(xx);
k=1:N;l=1:N;h(k,1)=0;
f1=(u.*u+v.*v);f1=pi*f1/(w*f*f);o=2*pi/(w*f);
for k=1:N
    for l=1:N
        del=0
        for m=1:M
            for n=1:M
                if A(m,n)~=0
                    
del=A(m,n)*exp(j*(z(m,n)*f1(k,l)-o*(x(m,n)*u(k,l)+y(m,n)*v(k,l))));
                 h(k,l)=h(k,l)+del;
                end
            end
        end
    end
end
save h64r1h1;

clear;clc;
load h64rdian5.mat;
dx=0.5*cm;dy=0.5*cm;
figure;imshow(A);
ex=exp(j*2*pi*(dx*u+dy*v)/(w*f));
h=3+real(h.*ex);
figure;imshow(h);
rec=fftshift(fft2(h));
rec=mat2gray(rec);
figure:imshow(rec);