function [h,hlines] = plotFieldMap(x, x_step, y, y_step, r, M, B_0, magFunc, wavelength, fieldPlotChoice, supMagParams)
%% [h,hlines] = plotFieldMap(x, x_step, y, y_step, r, M, B_0, magFunc, wavelength, fieldPlotChoice, supMagParams)
%
% Calcul et tracassage des lignes de champ et de la cartographie du champ
% Les neutrons se propagent selon x et l'axe du fil est selon z
%
% Arguments :
% ===========
% - x : +- la longeur de la zone calculée par rapport au centre du fil(en m)
% - x_step: pas de calcule selon x (en m)
% - y : +- la largeur de la zone calculée par rapport au centre du fil(en m)
% - y_step : pas de calcule selon y (en m)
% - r : rayon du fil (en m)
% - M : aimantation à saturation du fil (en T)
% - B_0 : champ magnétique appliqué (en T)
% - @magFunc: fonction de r_wire, M, B_0, X et Y calculant le champ
%                            r_wire -> nombre
%                            M      -> vecteur xyz
%                            B_0    -> vecteur xyz
%                            X      -> matrice position de dimention 2 ou 3
%                            Y      -> matrice position de dimention 2 ou 3
%                       et doit retourner [U,V,W] ou U, V et W sont des
%                       matrices de même dimension et taille que X et Y et
%                       représentent les composantes du champ selon x, y et
%                       z respectivement.
% - wavelength : Longeur d'onde en m
% - fieldPlotChoice : permet de choisir ce qui est tracé sur la colormap
%   = 1 -> Valeur absolue du champ
%   = 2 -> Critère d'adiabaticité, Omega_L/Omega_B
% - supMagParams : Liste de paramètres suplémentaires transmit à magFunc
%
% Exemple :
% =========
%
%  [h,hlines] = plotFieldMap(400e-6, 0.5e-6, 400e-6,0.53e-6, 1e-4, [0.1,0,0], [0.001,0.00,0.0], @wirePlanField_vector, 4e-10, 1);
%

% Unité d'affichage µm
unitFactor = 1e6;

% On compètes les arguments optionnels
if nargin < 9
    wavelength = 4e-10;
end
if nargin < 10
    fieldPlotChoice = 1;
end
if nargin < 11    
    supMagParams = [];
end


% Création des grilles de de coordonnée X,Y et Z
[X,Y,Z] = meshgrid(-x:x_step:x, -y:y_step:y,0:1);


% Calcul du champ
[U, V, W] = feval(magFunc, r, M,B_0, X, Y, supMagParams);


% Norme du champ
B_abs = (U.^2 + V.^2 + W.^2).^(0.5);
if fieldPlotChoice == 2
    omegaLB = ratioOmegaLB(U, V, W, B_abs, wavelength, x_step, 2);
end

X = X*unitFactor;
Y = Y*unitFactor;
Z = Z*unitFactor;


%% plot des lignes de champ
% peu planter à cette étape si le pas est trop grand.
% Erreur du genre -> "Attempted to access startgrid(12,0); index must be a 
%                     positive integer or logical.")
h = figure;
%Y2 = Y;
%Y2(:) =0;
hlines = streamslice(X,Y,Z,U,V,W,[],[],0);
set(hlines,'LineWidth',1,'Color','black')


%% plot de la color map
hold on;
if fieldPlotChoice == 1 % On trace la valeurs absolue du champ
    slice(X,Y,Z,B_abs,[],[],0);
elseif fieldPlotChoice == 2 % Ou bien on trace le rapport OmegaL/OmegaB
    omegaLB(isnan(omegaLB)) = 1e-8 ;
    omegaLB(isinf(omegaLB)) = 1e10 ;
    slice(X(:,2:end-1,:),Y(:,2:end-1,:),Z(:,2:end-1,:),log10(omegaLB),[],[],0);
end
shading interp
daspect([1 1 1])
%colormap(hot)
colormap(jet)
title('3D field map','FontSize',20,'FontWeight','bold')
ylabel('y (µm)','FontSize',16)
xlabel('x (µm)','FontSize',16)
xlim([-x,x]*unitFactor)
ylim([-y,y]*unitFactor)
%zlabel('z (mm)','FontSize',16)


%view(90,90)
end