% This program caculate the phase function of the input image of a 2f system.

function [phase]=PhaseSL(Spinput,Froutput)
input=Spinput;
[m,n]=size(input);
a=Froutput;
E2=sqrt(a);
b=Spinput;
E1=sqrt(b);
iter=1;         % 迭代次數(shù)
iterations=100;
E_1=rand(m);
while(iter<iterations+1)
    E_plane_2=fftshift(fft2(E1));
    E_2=E_plane_2;
    I_2=E_2.*conj(E_2);
    E_plane_2_new=E2.*exp(i.*angle(E_plane_2));
    E_plane_1=fftshift(ifft2(E_plane_2_new));
    phase=angle(E_plane_1);
    E_1=E1.*exp(i.*angle(E_plane_1));
    iter=iter+1;
end

% figure
% imshow(mat2gray(phase));

% save 'phase.mat',phase;

% figure(2)
% c=E1.*exp(i.*phase);
% d=fft2(c);
% I=d.*conj(d);
% imshow(mat2gray(log(1+I)))