function [U,V,W] = wireVortexField_vector_V2(r, M, B_0, X, Y, supMagParams)
%%% Calcule du champ magnétique pour une aimantation sous forme de vortex
% au sein du fil en tout point de l'espace. Cette version 2 ajoute un coeur
% au vortex.
% X et Y peuvent être données sous forme de matrices
% supMagParams = [h, Delta] où h et Delta sont les paraètres du vortex
% - h : proportion de la composante vertical de l'aimantation à l'infini. 
%       0 -> l'aimantation est dans le plan xy
%       1 -> l'aimant est selon z
%       Attention, h n'est pas totalement linéaire
% - Delta : Taille caractéristique du coeur du vortex
%
% $Revision: 14 $
% $Author: lenertz $
% $Date: 2015-07-27 18:22:01 +0200 (lun., 27 juil. 2015) $

    % Paramètres à externaliser
    % h normalisé entre 0 et 1
    if nargin < 6;
        h = 0;
        Delta = 1e-6;
    else
        h = supMagParams(1);
        if numel(supMagParams) < 2;
            Delta = 1e-6;
        else
            Delta = supMagParams(2);
        end
    end
   
    % TODO h pourrait être lié à l'application de B_0 (entre -0.0175 et 1,1321)
    % 
    h = h*1.1495 - 0.0174;
    

<<<<<<< .mine
    h = 0.9; % Pourrait être lié à l'application de B_0
    Delta = 1e-6; % longeur caractéristique (~largeur de parois de domaine)
    %delta_x = 0; % décentrage du vortex
=======
    %Delta = 5e-6; % longeur caractéristique (~largeur de parois de domaine)
    
    % TODO décentrage du vortex
    %delta_x = 0; % 
>>>>>>> .r16
    %delta_y = 0;


    B_abs = norm(M);
    % Inverse le sens du vortex si une composante de champ est appliqué selon -z
    if B_0(3) < 0;
        B_abs = -B_abs;
    end
    % fonction de rotation du reférenciel à rajouter
    %angle = atan2(norm(cross([1,0,0],H)),dot([1,0,0],H));
    
    % if x is given as an array an not y or the opposit
    if numel(X) ~= numel(Y); 
        if numel(Y) == 1;
            y = Y;
            Y = X;
            Y(:) = y;
        elseif numel(X) == 1;
            x = X;
            X = Y;
            X(:) = x; 
        end
    end
    
 
    
    
    
    % Initialisation des matrices de champs
    U = X;
    U(:)=0;
    V = U;
    W = U;
    
    % On cherche les positions dans et hors du fil
    position_r = (X.^2 + Y.^2).^(1/2);
    in_position = find(position_r <= r & position_r > 0);
    %out_position = find(position_r > r);
    
    A1 = 4.108 - 7.628/(2.989 - h);
    A2 = 1.387 +0.0717/(1.173 - h);
    A3 = 0.230 + 0.505/(1.144 - h);
   
    Rho = U;
    Theta = U;
    Rho(in_position) = position_r(in_position) / (Delta * A3);
    Theta(in_position) = A1 * (sinh(Rho(in_position))) ./ (A2 + sinh(Rho(in_position)));
    Theta(isnan(Theta)) = A1 ;
    
    % Dans le fil
    U(in_position) = -B_abs .* sin(Theta(in_position)) .* Y(in_position) ./ position_r(in_position);
    V(in_position) =  B_abs .* sin(Theta(in_position)) .* X(in_position) ./ position_r(in_position);
    W(in_position) =  B_abs .* cos(Theta(in_position));
    
    % On ajoute le champ appliqué
    U = U + B_0(1);
    V = V + B_0(2);
    W = W + B_0(3);

    
end