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

lambda = 12398.42./lcalc.reflectivity(1,:);
twok = 12.56637./lambda;      
factor = 0.00000540556*lambda.^2;% /2 for sqrt approximation

alfa0y=cosd(lsim.angle(1));
lcalc.Q(:,:,1) = [twok;twok]* sind(lsim.angle(1)); % 1 is vacuum

%% LOOP FOR INDICES CALCULATIONS
for jj=1:length(layer)
    eps_para = 1 - factor.*layer(jj).density(2).*(SF.cf(jj+1,:) + layer(jj).AS(2)*SF.cf2(jj+1,:));
    eps_perp = 1 - factor.*layer(jj).density(2).*SF.cf(jj+1,:);
    lcalc.Q(:,:,jj+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(:,:,jj+1) = [exp(1i*lcalc.Q(1,:,jj+1)*layer(jj).thickness(2)) ; exp(1i*twok.*sqrt(eps_perp.*(1-alfa0y^2./eps_para))*layer(jj).thickness(2))];
end
   
%% MAIN LOOP FOR REFLECTIVITY CALCULATION
for jj=length(sequence)-1:-1:1 % jj is position of the layer ;  sequence(jj) is type of the layer
    
    r = ((lcalc.Q(:,:,sequence(jj))-lcalc.Q(:,:,sequence(jj+1))) ./ (lcalc.Q(:,:,sequence(jj))+lcalc.Q(:,:,sequence(jj+1))))...
        .* exp(-0.5*(lcalc.Q(:,:,sequence(jj)).*lcalc.Q(:,:,sequence(jj+1)))*layer(sequence(jj+1)-1).roughness(2)^2);
    lcalc.R(:,:,jj) = (r + lcalc.R(:,:,jj+1) .* lcalc.exp_phase(:,:,sequence(jj+1))) ./ (1 + r .* lcalc.R(:,:,jj+1) .* lcalc.exp_phase(:,:,sequence(jj+1)));
 
end
%% RESULTS

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

end