function [lcalc]=F0F2_Parratt_qf(lcalc,lsim,layer,SF,sequence)

twok = 1.01355e-3*lcalc.reflectivity(1,:); % 2*2PI/12398.4 a line

factor = 830.9471./lcalc.reflectivity(1,:).^2;% 12398*12398*0.0000540556 a line
theta_0=asin(lsim.Q*6199.2./lcalc.reflectivity(1,:)); %a line
alfa0y=cos(theta_0);
lcalc.Q(:,:,1) = [twok.* sin(theta_0);twok.* sin(theta_0)]; % 1 is vacuum

for k=1:length(layer)
    eps_para = 1 - factor.*layer(k).density(2).*(SF.cf(k+1,:) + layer(k).AS(2)*SF.cf2(k+1,:));
    eps_perp = 1 - factor.*layer(k).density(2).*SF.cf(k+1,:);
    lcalc.Q(:,:,k+1 ) = [ twok .* sqrt(eps_perp - alfa0y.^2) ; twok .* sqrt((eps_para - alfa0y.^2)./(eps_perp.*eps_para))]; % array jj+1 to account for vacuum layer
    lcalc.exp_phase(:,:,k+1) = [exp(1i*lcalc.Q(1,:,k+1)*layer(k).thickness(2)) ; exp(1i*lcalc.Q(2,:,k+1)*layer(k).thickness(2))];%
end

for k=length(sequence)-1:-1:1 % jj is position of the layer ;  sequence(jj) is type of the layer
    
    r = ((lcalc.Q(:,:,sequence(k))-lcalc.Q(:,:,sequence(k+1))) ./ (lcalc.Q(:,:,sequence(k))+lcalc.Q(:,:,sequence(k+1))))...
        .* exp(-0.5*(lcalc.Q(:,:,sequence(k)).*lcalc.Q(:,:,sequence(k+1)))*layer(sequence(k+1)-1).roughness(2)^2);
    lcalc.R(:,:,k) = (r + lcalc.R(:,:,k+1) .* lcalc.exp_phase(:,:,sequence(k+1))) ./ (1 + r .* lcalc.R(:,:,k+1) .* lcalc.exp_phase(:,:,sequence(k+1)));
 
end

lcalc.M_final(1,1,:) = lcalc.R(1,:,1);
lcalc.M_final(2,2,:) = lcalc.R(2,:,1);

end