for k=1:length(x1), % calculer Theta et Delta
 
    distance=sqrt((x1-x1(k)).^2+(y1-y1(k)).^2+(z1-z1(k)).^2);
    voisin = find(distance<4*hmesh) ;
    
    %             x1_local = zeros(1,size(voisin));
    %             y1_local = zeros(1,size(voisin));
    %             z1_local = zeros(1,size(voisin));
    
    %Matrices de Passage P1 et P2
    Theta(k) = angle(x1(k)+1i*y1(k)) ;
    Delta(k) = 0 ; %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
    
    xyz1_local = P2*P1*[x1(voisin)';y1(voisin)';z1(voisin)'];
    xyz1_klocal = P2*P1*[x1(k)';y1(k)';z1(k)'];
    
    x1_local = xyz1_local(1,:)';
    y1_local = xyz1_local(2,:)';
    z1_local = xyz1_local(3,:)';
    
    x1_klocal = xyz1_klocal(1,:);
    y1_klocal = xyz1_klocal(2,:);
    z1_klocal = xyz1_klocal(3,:);
    
    %vecteur normal dans la base tourné
    A = [ones(length(voisin),1),x1_local-x1_klocal,y1_local-y1_klocal,(x1_local-x1_klocal).*(y1_local-y1_klocal), ...
        (x1_local-x1_klocal).*(x1_local-x1_klocal), (y1_local-y1_klocal).*(y1_local-y1_klocal) ];
    B = z1_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 (originale)
    vect_n(k,:)= P1'*P2'*[vect_n_local(k,1); vect_n_local(k,2); vect_n_local(k,3)];
    
end