A question to understand the deconvblind codes proposed by the Richardson Lucy algorithm

I met a question after reading the codes of Richardson Lucy algorithm within "deconvblind", parts of codes as folows:

…
% 2.b  Make core for the LR estimation

    CC = corelucy(Y,psf2otf(B,sizeI),DAMPAR22,wI,READOUT,1,idx,[],[]);
    % 2.c Determine next iteration image & apply positivity constraint
    J{3} = J{2};
    H = psf2otf(P{2},sizeI);
   *scale = real(ifftn(conj(H).*fw)) + sqrt(eps);*
    *J{2} = max(Y.*real(ifftn(conj(H).*CC))./scale,0);*
    clear scale;
    J{4} = [J{2}(:)-Y(:) J{4}(:,1)];
    clear Y H;
    % 2.d Determine next iteration PSF & apply positivity constraint + normalization
    P{3} = P{2};
    H = fftn(J{3});
    *scale = otf2psf(conj(H).*fw,sizePSF) + sqrt(eps);
    P{2} = max(B.*otf2psf(conj(H).*CC,sizePSF)./scale,0);*
    clear CC H;
    sumPSF = sum(P{2}(:));
    P{2} = P{2}/(sumPSF + (sumPSF==0)*eps);
    ….

As shown in the first Bold part

scale = real(ifftn(conj(H).*fw)) + sqrt(eps);
    J{2} = max(Y.*real(ifftn(conj(H).*CC))./scale,0);

it seems that the computation of J{2} isn't conformity with the updated formula of the ramped Richardson Lucy algorithm or original Richardson Lucy algorithm; the updated formula of ramped Richardson Lucy algorithm as follows:

r_k = h_{k+1}*f_k;

U_k=-2/T^2[glog(r_k/g) – r_k +g];

bar(U_k) = min(U_k,1);

f_{k+1}=(1/sum(h_{k+1}))f_k x h_{k+1}*CC

Someone may give me a favor? Or someone may give me some suggestions or data for understanding the codes of this part? Thank you very much.