function [Xrec , iterres , iterchange] = ER(Diagonal,Y_abs_noisy,imagetype,Index,X0,itermax)

%% parameters
i = sqrt(-1);
N11 = Index(2,1);
N22 = Index(2,2);
N1start = Index(3,1);
N1end = Index(3,2);
N2start = Index(4,1);
N2end = Index(4,2);
P = size(Diagonal,3);

%% initiation
Xk = X0;
iterres = []; % record residual
iterchange = []; % record ||X{k}-X{k-1}||/||X{k-1}||

Y_abs_noisy_fro = 0; % Euclidean norm of Y_abs_noisy
for k = 1 : P
Y_abs_noisy_fro = Y_abs_noisy_fro + norm(Y_abs_noisy(:,:,k),'fro')^2;
end
Y_abs_noisy_fro = sqrt(Y_abs_noisy_fro);

for iter = 1 : itermax
TemMatrix = Xk;
%% compute Pf{Xk}
XkP = Xk;
for k = 1 : P
FXk = fft2(Diagonal(:,:,k).*XkP);
FXkP = Y_abs_noisy(:,:,k).*(FXk./abs(FXk));
XkP = ifft2(FXkP)./Diagonal(:,:,k);
end
%% Po Pf {Xk}
Xk = zeros(N11,N22);
Xk(N1start:N1end,N2start:N2end) = projection(XkP(N1start:N1end,N2start:N2end),imagetype);
%% record residual
intensityerr = 0;
for k = 1 : P
intensityerr = intensityerr + norm(abs(fft2(Diagonal(:,:,k).*Xk))-Y_abs_noisy(:,:,k),'fro')^2;
end
intensityerr = sqrt(intensityerr) / Y_abs_noisy_fro;
iterres(iter) = intensityerr;
%% record the change of Xk
iterchange(iter) = norm(TemMatrix-Xk,'fro')/norm(TemMatrix,'fro');
fprintf('%6.0f residual = %6.4f p change = %6.6f p \n',iter,100*intensityerr,100*iterchange(iter))
end

Xrec = Xk;