function [lcalc]=F0F1_ef_CC_novect(lcalc,lsim,layer,SF,sequence)
% Code for Elzo'formalism, fix energy. No vectorization so as to ...
% keep some clarity in the code. Same results than vectorized code...
% but slower.

wave_k =5.0677e-4*lsim.energy(1);
kz=zeros(1,length(layer)+1);
eps_xy=zeros(1,length(layer)+1);
eps_xz=zeros(1,length(layer)+1);
eps_yz=zeros(1,length(layer)+1);
eps_mag=zeros(1,length(layer)+1);
eps=zeros(length(layer)+1,1)+1;
u=zeros(3,length(layer)+1);
A=zeros(4); A_phi=zeros(4);
B=zeros(4); B_phi=zeros(4); 
P=zeros(4);P_phi=zeros(4);

thick=zeros(1,length(layer)+1);
rough=zeros(1,length(layer)+1);

for ll=1:length(layer)+1
    
  if ll>1
    cappa = 830.9471/lsim.energy(1)^2*layer(ll-1).density(2);
    eps(ll) = 1 - cappa*SF.cf(ll,1);
    eps_mag(ll) = cappa*layer(ll-1).MMS(2)*SF.mf(ll,1);
    u(:,ll)  = [sind(layer(ll-1).phi(2))*cosd(layer(ll-1).gamma(2));...
          sind(layer(ll-1).phi(2))*sind(layer(ll-1).gamma(2));...
          cosd(layer(ll-1).phi(2))];
    eps_xy(ll) = - 1i *u(3,ll)*eps_mag(ll);
    eps_xz(ll) =   1i *u(2,ll)*eps_mag(ll);
    eps_yz(ll) = - 1i *u(1,ll)*eps_mag(ll);
    thick(ll)=layer(ll-1).thickness(2);
    rough(ll)=layer(ll-1).roughness(2);  
  end
end
for ii=1:length(lcalc.reflectivity(1,:))
    
    s=eye(4); s_phi=eye(4);
    alphay_0 = cos(lcalc.reflectivity(1,ii));
    
    for kk=1:1:length(sequence)
        
        k=sequence(kk);
                
        alphay = cos(lcalc.reflectivity(1,ii))/sqrt(eps(k)); 
        alphaz = sqrt(1-alphay^2);

        kz(k) = wave_k * sqrt(eps(k)) * alphaz;

n_right_down = sqrt( eps(k) - eps_mag(k)*( u(2,k)*alphay +u(3,k)*alphaz));
n_left_down  = sqrt( eps(k) + eps_mag(k)*( u(2,k)*alphay +u(3,k)*alphaz));
n_right_up   = sqrt( eps(k) - eps_mag(k)*( u(2,k)*alphay -u(3,k)*alphaz));
n_left_up    = sqrt( eps(k) + eps_mag(k)*( u(2,k)*alphay -u(3,k)*alphaz));
        
alphay_right_down  = alphay_0/n_right_down;
alphaz_right_down  = sqrt(1-alphay_right_down^2);
alphay_left_down   = alphay_0/n_left_down;
alphaz_left_down   = sqrt(1-alphay_left_down^2) ;
alphay_right_up    = alphay_0/n_right_up;
alphaz_right_up    = sqrt(1-alphay_right_up^2);
alphay_left_up     = alphay_0/n_left_up;
alphaz_left_up     = sqrt(1-alphay_left_up^2);
        
A(1) = -1 -1i*eps_xy(k)*alphaz_right_down  -1i*eps_xz(k)*alphay_right_down;
A(2) =  - eps_xy(k) + 1i*alphaz_right_down -1i*eps_yz(k)*alphay_right_down;
A(3) =  -1i*n_right_down*A(1);
A(4) =  -alphaz_right_down*n_right_down*A(1);
        
A(5) =  1 -1i*eps_xy(k)*alphaz_left_down -1i*eps_xz(k)*alphay_left_down;
A(6) =  eps_xy(k) + 1i*alphaz_left_down  -1i*eps_yz(k)*alphay_left_down;
A(7) =  1i*n_left_down*A(5);
A(8) =  -alphaz_left_down*n_left_down*A(5);
        
A(9) =  -1 -1i*eps_xy(k)*alphaz_right_up -1i*eps_xz(k)*alphay_right_up;
A(10) =  -eps_xy(k) - 1i*alphaz_right_up -1i*eps_yz(k)*alphay_right_up;
A(11) =  -1i*n_right_up*A(9);
A(12) =  alphaz_right_up*n_right_up*A(9);
        
A(13) = 1 -1i*eps_xy(k)*alphaz_left_up -1i*eps_xz(k)*alphay_left_up;
A(14) =  eps_xy(k) - 1i*alphaz_left_up -1i*eps_yz(k)*alphay_left_up;
A(15) =  1i*n_left_up*A(13);
A(16) =  alphaz_left_up*n_left_up*A(13);
        
A_phi(1) = -A(5);
A_phi(2) = eps_xy(k) +1i*alphaz_left_down +1i*eps_yz(k)*alphay_left_down;
A_phi(3) =  A(7);
A_phi(4) =  -A(8);
        
A_phi(5) =  -A(1);
A_phi(6) =-eps_xy(k) +1i*alphaz_right_down +1i*eps_yz(k)*alphay_right_down;
A_phi(7) =  A(3);
A_phi(8) =  -A(4);
        
A_phi(9) =  -A(13);
A_phi(10) =  -eps_xy(k) -1i*alphaz_left_up -1i*eps_yz(k)*alphay_left_up;
A_phi(11) =  A(15);
A_phi(12) =  -A(16);
        
A_phi(13) = -A(9);
A_phi(14) = -eps_xy(k) -1i*alphaz_right_up +1i*eps_yz(k)*alphay_right_up;
A_phi(15) =  A(11);
A_phi(16) = -A(12);
        
A=A/sqrt(2)/eps(k);
A_phi=A_phi/sqrt(2)/eps(k);
        
        if k>1
            j=sequence(kk-1);
            F= A\B;
            
            F_phi= A_phi\B_phi;
            rugosn = exp(-((-kz(k)+kz(j))^2)*rough(k)^2 /2);
            rugosp = exp(-(( kz(k)+kz(j))^2)*rough(k)^2 /2);
            
            R=[rugosn rugosn rugosp rugosp;
               rugosn rugosn rugosp rugosp;
               rugosp rugosp rugosn rugosn;
               rugosp rugosp rugosn rugosn];
            
            F_rough=F.*R; F_rough_phi = F_phi.*R;
                   
            if k < length(layer)+1

phase=wave_k*thick(k);
P(1) = exp( 1i*phase*n_right_down*alphaz_right_down);
P(6) = exp( 1i*phase*n_left_down*alphaz_left_down);
P(11)= exp(-1i*phase*n_right_up*alphaz_right_up);
P(16)= exp(-1i*phase*n_left_up*alphaz_left_up);
P_phi(1)  = P(6);  P_phi(6)  = P(1);
P_phi(11) = P(16); P_phi(16) = P(11);
               
 s = P*F_rough*s; s_phi = P_phi*F_rough_phi*s_phi;
            
            else
                M = F_rough*s;
                M_phi = F_rough_phi*s_phi;
            end         
        end
        
        B=A;
        B_phi=A_phi;
        
    end
    lcalc.M_final(:,:,ii)=(M(3,3)*M(4,4)-M(4,3)*M(3,4))\...
        ([-1 1;-1i -1i]*...
        [-M(4,4)*M(3,1)+M(3,4)*M(4,1)...
         -M(4,4)*M(3,2)+M(3,4)*M(4,2);...
          M(4,3)*M(3,1)-M(3,3)*M(4,1)...
          M(4,3)*M(3,2)-M(3,3)*M(4,2) ]*[-1 1i;1 1i]*0.5) ;
    
    lcalc.M_final_phi(:,:,ii)=(M_phi(3,3)*M_phi(4,4)-M_phi(4,3)*M_phi(3,4))\...
       ([-1 1;-1i -1i]* [-M_phi(4,4)*M_phi(3,1)+M_phi(3,4)*M_phi(4,1) ...
        -M_phi(4,4)*M_phi(3,2)+M_phi(3,4)*M_phi(4,2);
        M_phi(4,3)*M_phi(3,1)-M_phi(3,3)*M_phi(4,1) ...
        M_phi(4,3)*M_phi(3,2)-M_phi(3,3)*M_phi(4,2) ]*[-1 1i;1 1i]*0.5);      
end
end