%MARF algorithm for the generation of Holograms of arbitrary designed beam shape
close all
clear all
clc

m=0.5;%mixing parameter when iterating,if m=1,then it is totally GS algorithm   混合参数进行迭代时,如果m = 1,那么它是完全GS算法
row=1280;%resolution of DOE in x direction
column=1024;%resolution of DOE in y direction
[x,y]=meshgrid(1:row,1:column);
R=0.0003;%curvature of quadratic phase profile  二次相轮廓的曲率
Alpha=0.5;%aspect ratio of quadratic phase profile  二次相轮廓的长宽比
Phase=zeros(column,row);%initial phase  初始相位i
K1=zeros(column,row);

for j=1:row
    for k=1:column
        K1(k,j)=4*R*(Alpha*(j-row/2)^2+(1-Alpha)*(k-column/2)^2);%initial quadratic phase  最初的二次相
        Phase(k,j)=mod(K1(k,j),2*pi);
    end
end
figure,imshow(mod(Phase.*255./(2*pi),256),[0 255]);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%初始相位图计算

F=0.04;         %矩形光斑区域占相位图比例
w=F*row;
h=w*30/50;         %矩形长w,宽h
for j=1:row
    for k=1:column
        if j>row/2-w/2&j<row/2+w/2&k>column/2-h/2&k<column/2+h/2
            Photo(k,j)=255;
        else
            Photo(k,j)=0;
        end
    end
end
A=Photo;%target pattern
goals=double(A);%target pattern
figure(),surfc(x,y,goals), title('Target intensity distribution'), colormap jet, axis equal, axis tight, view([45, 45]), shading interp

SigRegion=((x-(row+1)/2).^2+(y-(column+1)/2).^2)./(2*w).^2;  %信号区域尺寸是2w大小
SigRegion(SigRegion<=1)=1;
SigRegion(SigRegion>1)=0;

A0=exp(-((x-row/2).^2+(y-column/2).^2)./(750^2));%Gaussian incident光斑大小该如何设置？
Efficiency=[];%diffraction efficiency
Accuracy=[];%accuracy of the generated laser pattern

tic    %tic/toc配合使用来计时
h=waitbar(0,'Caculating ...');
index=0;
loop=100;%times of iterative
for k=1:loop
    f=exp(1i.*Phase).*A0;
    f1=fftshift(fft2((fftshift(f))));
    f1=f1.*sqrt(sum(sum(abs(A0).^2))/sum(sum(abs(f1).^2)));
%     Factor=sqrt(sum(sum(abs(A0).^2))/sum(sum(abs(f1).^2)));
%     disp(Factor)
    f2=(m*sqrt(goals.*SigRegion)+(1-m)*abs(f1.*~SigRegion)).*exp(1i.*(angle(f1)));
    f3=fftshift(ifft2(ifftshift(f2)));
    Phase=angle(f3);
    Efficiency=[Efficiency sum(sum(abs(f1.*goals./255).^2))/sum(sum(abs(f1).^2))];
    fTarget=abs(f1).*goals./255;
    goals_u=1-((max((fTarget(:))))^2-(min((fTarget(fTarget~=0))))^2)./((max((fTarget(:))))^2+(min((fTarget(fTarget~=0))))^2); %均匀性计算
    Accuracy=[Accuracy goals_u];
    index=index+1;
    waitbar(index./loop); %进度条

end
close(h); %关闭进度条
toc

Intensity=abs(f1).^2./sum(sum(abs(f1).^2))*100000;
figure(),surfc(x,y,Intensity), title('Simulated intensity distribution'), colormap jet, axis equal, axis tight, view([45, 45]), shading interp;
figure(),imshow(abs(f1),[])
figure,plot(Efficiency,'*'),title('efficiency');%Efficiency varies with iterative
figure,plot(Accuracy,'*'),title('accuracy');%Accuracy varies with iterative
Phase=Phase-min(min(Phase));
Phase=mod(Phase,2*pi);
hologram=round(Phase.*255./(2.*pi));
figure(),imshow(hologram,[0 255]);
imwrite(uint8(hologram),'MRAF.bmp');
% Addition=zeros(1080,420);
% hologramStandard=[Addition hologram Addition];
% figure,imshow(hologramStandard,[]);
% imwrite(uint8(hologram),'E:\CGH\MRAF\letter\APL\apl30点\holo\L1920_03.bmp');