function [BL_ratio] = ratioOmegaLB(Bx, By, Bz, B, wavelength, x_step,x_axis)
%% Calcule le repport omegaL/OmegaB 
%[BL_ratio] = ratioOmegaLB2(Bx, By, Bz, B, x_step,x_axis)
%
% Prend en arguments :
%  -> Bx, By et Bz, 3 matrices de même taille comprenant les composantes du champ
%     magnétique selon x, y et z. Si le champ est décrit en 1D ou 2D il
%     faut remplacer Bz (et By si 1D) par 0 
%  -> B, matrice de même dimension que Bx et décrivant la norme du champ.
%     Elle peut aussi être calculé au coeur de la fonction en remplacant
%     B par 0.
%  -> x_step, le pas entre chaque points en x.
%  -> a_axis, argument optionnel. A spécifié si la variation en x se
%     trouve dans le second indices des matrices (par exemple si X,Y et Z
%     sont générés par meshgrid). Il faut alors spécifié la dimension
%     dans la quelle x varie. Peut aussi être utilisé pour calucler
%     les gradients selon d'autres directions que x.
%
% retourne une matrice de même taille que Bx, excepté selon la direction de
% variation (voir x_axis), selon laquelle la longeur est réduite de 2.

    %% Initialisation des paramètres
    
    % Remplace By par une matrice de taille adaptée si aucune composante
    % selon y n'a été donnée (By = 0).
    if By == 0;
        By = Bx;
        By(:) = 0;
    end
    
    % Même chose, mais selon z
    if Bz == 0;
        Bz = Bx;
        Bz(:) = 0;
    end
        
    % Génère la norme si elle n'a pas encore été calculée
    if B == 0;
        B = (Bx.^2 + By.^2 + Bz.^2).^0.5;
    end
    
    % Nombre de dimensions
    nb_dim = ndims(B);

    % Inverse les axes, utile si X,Y,Z sont généré par meshgrid 
    if nargin >= 6;    
        % Genère un vecteur de permutation en fonction du nombre de
        % dimensions
        P = 1:nb_dim;
        P(1)= x_axis;
        P(x_axis) = 1;
        
        Bx = permute(Bx, P);
        By = permute(By, P);
        Bz = permute(Bz, P);
        B = permute(B, P);
    end
    
    % tous les angles sont calculés en radian
    %% carto pulsation de Larmor (avec gamma_n/2pi = 29.1646943 MHz T-1)
    omega_L=abs(B)*2*pi*29.1646943e6;
    
    %% calcul de la "fréquence" du changement de champ
    %wavelength = 0.4e-9; % longeur d'onde en m
    masse_neutron = 1.67493e-27; %masse du neutron en Kg
    h_plank = 6.6260688e-34; % constante de plank en J.s
    neutron_speed = h_plank / (masse_neutron * wavelength); % en m/s

    Bx_norm=Bx./B;
    By_norm=By./B;
    Bz_norm=Bz./B;

    
    Bxyz_norm = cat(nb_dim+1, Bx_norm, By_norm, Bz_norm);
    if nb_dim == 1; % cas d'une matrice 1D
        dThetaB = dot(Bxyz_norm(1:end-2,:),Bxyz_norm(3:end,:),2);
    elseif nb_dim == 2; % cas d'une matrice 2D
        dThetaB = dot(Bxyz_norm(1:end-2,:,:),Bxyz_norm(3:end,:,:),3);
    else % Cas d'une matrice 3D
        dThetaB = dot(Bxyz_norm(1:end-2,:,:,:),Bxyz_norm(3:end,:,:,:),4);
        
    end
    
    % Pour corriger des bug d'arrondis numérique entrainant l'apparition de
    % valeur très légèrement supérieur à 1 et résultant en des solutions
    % complexes par la suite.
    index = find(dThetaB > 1);
    dThetaB(index) = 1;

    
    dThetaB = acos(dThetaB);
    
    omega_B = dThetaB * (neutron_speed / (2 * x_step));

    %% Calcul du ratio et mise en forme de la matrice de sortie
    if nb_dim == 1; % cas d'une matrice 1D
        BL_ratio=omega_L(2:end-1)./omega_B;
    elseif nb_dim == 2; % cas d'une matrice 2D
        BL_ratio=omega_L(2:end-1,:)./omega_B;
    else % Cas d'une matrice 3D
        BL_ratio=omega_L(2:end-1,:,:)./omega_B;
    end

    % Les axes sont réordonnée de la même manière qu'en entré
    if nargin >= 6 && x_axis == 2;
        BL_ratio = permute(BL_ratio, P);
    end
end 