function [U,V,W] = wireVortexField_vector_V4(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.
% 
% Ajoute par rapport à la V3 un paramètre z_ratio controlant la proportion
% initial de l'aimantation selon z
%
% $Revision: 16 $
% $Author: lenertz $
% $Date: 2015-07-30 17:37:06 +0200 (jeu., 30 juil. 2015) $


    r_core = supMagParams(1);
    Delta = supMagParams(2);
    mixe_ratio = supMagParams(3);
    M1 = mixe_ratio*M;
    M2 = (1-mixe_ratio)*M;
    
% [U1,V1,W1] = wirePlanField_vector(r, M*mixe_ratio, B_0*mixe_ratio, X, Y, supMagParams);
% 
% [U2,V2,W2] = wireVortexField_vector_V3(r, M*(1-mixe_ratio), B_0*(1-mixe_ratio), X, Y, supMagParams);
%   
% 
% U = U1 + U2;
% V = V1 + V2;
% W = W1 + W2;



    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;
    U_temp = U;
    V_temp = U;
    B_abs = U;
    ratio = U;
    
    % On cherche les positions dans et hors du fil
    position_r = (X.^2 + Y.^2).^(1/2);
    in_position = find(position_r <= r);
    out_position = find(position_r > r);
    
    % Dans le fil
    U(in_position) = M1(1);
    V(in_position) = M1(2);
    W(in_position) = M1(3);
    
    
    % Hors du fil
    % On calcule d'abord dans un repère ou le champ M1xy est selon x
    U_temp(out_position) = (X(out_position).^2 - Y(out_position).^2) .* norm(M1(1:2))*(r^2) ./ ((X(out_position).^2 + Y(out_position).^2).^2);
    V_temp(out_position) = 2 .* X(out_position) .* Y(out_position) .* norm(M1(1:2))*(r^2) ./ ((X(out_position).^2 + Y(out_position).^2).^2);
    % puis on procède à un changement de base
    theta = atan2(M1(2),M1(1)); % angle entre le champ M1 dans le plan xy avec l'axe x
    U(out_position) = cos(theta) * U_temp(out_position) + sin(theta) * V_temp(out_position);
    V(out_position) = cos(theta) * V_temp(out_position) - sin(theta) * U_temp(out_position); 
    
    
    
    %% vortex V3
    
     % 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);
    M_abs = norm(M2);
    % Inverse le sens du vortex si une composante de champ est appliqué selon -z
    if B_0(3) < 0;
        M_abs = -M_abs;
    end
    
    shell_position = position_r <= r & position_r >= r_core;
    core_position = position_r < r_core & position_r > 0;
 
    % Dans le fil
    ratio(core_position) = exp((position_r(core_position)/r_core - 1)/Delta);
    W(core_position) = M_abs .* ratio(core_position) + W(core_position);
    W(shell_position) = M_abs + W(shell_position);
    
    U(core_position) = -M_abs .* sin(acos(ratio(core_position))) .* Y(core_position) ./ position_r(core_position) + U(core_position);
    V(core_position) =  M_abs .* sin(acos(ratio(core_position))) .* X(core_position) ./ position_r(core_position) + V(core_position);
    
    
    
    
    
    B_abs(in_position) = (U(in_position).^2 + V(in_position).^2 + W(in_position).^2).^(0.5);
    U(in_position) = norm(M) * U(in_position) ./ B_abs(in_position);
    V(in_position) = norm(M) * V(in_position) ./ B_abs(in_position);
    W(in_position) = norm(M) * W(in_position) ./ B_abs(in_position);
    
    
    
    
    
    
    
    % On ajoute le champ appliqué
    U = U + B_0(1);
    V = V + B_0(2);
    W = W + B_0(3);



end