for k=1:length(x2),
     
    distance=sqrt((x2-x2(k)).^2+(y2-y2(k)).^2+(z2-z2(k)).^2);
    voisin = find(distance<4*hmesh) ;
    %     x2_local = zeros(1,size(voisin));
    %     y2_local = zeros(1,size(voisin));
    %     z2_local = zeros(1,size(voisin));
    
    %Matrices de Passage P1 et P2
        
    Theta(k) = angle(x2(k)+1i*y2(k)) ;
    Delta(k,1) = atan(sqrt(x2(k).^2+y2(k).^2)/z2(k)); %result in radians
    P1=[cos(Theta(k)), sin(Theta(k)), 0; -sin(Theta(k)), cos(Theta(k)), 0; 0, 0, 1]; %Rotation dans Z
    P2=[cos(Delta(k)), 0, -sin(Delta(k)); 0, 1, 0; sin(Delta(k)), 0, cos(Delta(k))]; %Rotation dans Y
    
    xyz2_local = P2*P1*[x2(voisin)';y2(voisin)';z2(voisin)'];
    xyz2_klocal = P2*P1*[x2(k)';y2(k)';z2(k)'];
    
    x2_local = xyz2_local(1,:)';
    y2_local = xyz2_local(2,:)';
    z2_local = xyz2_local(3,:)';
    x2_klocal = xyz2_klocal(1,:);
    y2_klocal = xyz2_klocal(2,:);
    z2_klocal = xyz2_klocal(3,:);
    
    
    %vecteur normal dans la base tourné
    A = [ones(length(voisin),1),x2_local-x2_klocal,y2_local-y2_klocal,(x2_local-x2_klocal).*(y2_local-y2_klocal), ...
        (x2_local-x2_klocal).*(x2_local-x2_klocal), (y2_local-y2_klocal).*(y2_local-y2_klocal) ];
    B = z2_local ;
    xxXX = A\B ;  %Backslash or matrix left division Matlab
    vect_n_local(k,:) = [-xxXX(2) , -xxXX(3) , 1] / sqrt(xxXX(2)^2+xxXX(3)^2+1) ; %La valeur entre la magnitud
    
    %vecteur normal dans la base global
    vect_n(k,:)= P1'*P2'*[vect_n_local(k,1); vect_n_local(k,2); vect_n_local(k,3)];
    
end