function [out]=travel3(M,B_0,r_wire, X, y, wavelength, magFunc, startPolar, angPhase, supMagParams)
%% Résolution de l'équa diff de la polarisation du neuutron selon x 
%
% Version 3 de travel. Résoud l'équation différentielle régissant le moment
% magnétique du neutron en fonction de x dans un champ H. Le champ
% magnétique est est donnée par une fonction de xy (entre autre). peremtant
% de calculer le champ en tout point de l'espace (voire doc
% computPolarization).
%
% Equation : 
% dM/Dx = gamma M x B / v
% pour les neutrons gamma = 29e6 (29MHz/T) (29kHz/mT)
%
% INPUT: M un vecteur d'aimantation [Mx,My,Mz] dans le fil
%        Chaque point correspond à 1 nm dans la paroi. L'épaisseur du film magnetique 
%        e est donc donnée par d=size(M(:,1))
%        B0 est le vecteur du champ de polarisation des neutrons 
%        (constant dans l'espace). Typiquement B0 = [0.001,0,0] T.
%        t angle de tilt sur la paroi (input en degrés)
%        lamb est la longueur d'onde du neutron exprimée en A
%        plotYN ets un flag (0 ou 1) pour indiquer s'il faut plotter la précession 
%        du neutron dans la couche
% NOTE: on travaille en mètres parce que les échelles pertinentes vont du nm au
% cm à cause de l'incidence rasante
%
%
% $Revision: 14 $
% $Author: lenertz $
% $Date: 2015-07-27 18:22:01 +0200 (lun., 27 juil. 2015) $


%% initialisation des paramètres
% quelques constantes de la physique utilisées dans la procédure de precession
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);
gamma_n =183.3e6;    % rapport gyromagnetique du neutron

% paramètres optionnels
if nargin < 8; startPolar = 1; end
if nargin < 9; angPhase = 0; end

    
    



% distance d'intégration, 0 est le centre du fil et
% polarisation du neutron incident (le long du champ au départ)
% en fonction du mode de calcule du champ
% if MagMod == 1 || MagMod == 3;
%     % On calcule sur une plus grande distance à cause des champ de fuite
%     x_0 = -50 * r_wire; % position de départ
%     x_f = 50*r_wire; % position final
%     if MagMod == 1;
%         [U,V,W] = wirePlanField_vector(r_wire, M, B_0, x_0, y);
%     elseif MagMod == 3;
%         [U,V,W] = wire3Field(r_wire, M, B_0, x_0, y);
%     end    
% elseif MagMod == 2;
%     x_0 = -2 * r_wire; % position de départ
%     x_f = 2*r_wire; % position final
%     [U,V,W] = wireVortexField_vector(r_wire, M, B_0, x_0, y);
% end

% vecteur champ magmétique en x_0
[u,v,w] = feval(magFunc, r_wire, M, B_0, X(1), y, supMagParams);
mu_0 = [u,v,w];
mu_0 = mu_0 / norm(mu_0);

%% Traitement de la polarisation et phase de départ
if startPolar == 1;
elseif startPolar == -1;
    mu_0 = -mu_0;
else
    
    polarAng = acos(startPolar);
    H_x0 = mu_0; % on retient la direction du champ
    
    % On cherche un deuxième vecteur non parallèle à mu_0
    if sum(cross(mu_0, [1,0,0])) == 0; % mu_0 déja selon x
        v2 = [0,1,0];
    else
        v2 = [1,0,0];
    end
    
    % On cherche le vecteur normal au plan formé par mu_0 et v2
    v3 = cross(mu_0, v2); % pas grave s'il n'est pas normé, rotGeneral s'en charge
    
    % application de la polarisation (rotation de mu_0 d'un angle polarAng 
    % autour d'un axe perpendiculaire)
    mu_0 = rotGeneral(mu_0,v3,polarAng);
    % Application de la phase (rotation de mu_0 autour du champ)
    mu_0 = rotGeneral(mu_0,H_x0,angPhase);    

end


%% calcul de la precession lors du trajet
OPTIONS = odeset('RelTol',1e-4,'AbsTol',1e-10);
[x,m] = ode45(@(x, mu) precessXY(x,mu, y, r_wire, M, B_0, neutron_speed, gamma_n, magFunc, supMagParams),X ,mu_0,OPTIONS);



%% visualisation des data sur un trajet 
% TODO: à mettre à jour car plus supporté par la version actuelle du code !!!!
mx = m(:,1);
my = m(:,2);
mz = m(:,3);


aux=size(mx);
%out=[mx(aux(1)),my(aux(1)),mz(aux(1))];
out = [x,m];
