function [out]=travel2(M,B_0,r_wire,y,lamb,plotYN, MagMod, angPhase)
% simulation de la precession d'un neutron dans un champ magnétique
% losqu'il traverse ou passe à coté d'un fil magnétique uniformément
% aimanté.
% 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: 11 $
% $Author: lenertz $
% $Date: 2015-07-10 14:32:19 +0200 (ven., 10 juil. 2015) $


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








% 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
mu_0 = [U,V,W];
%[U, V, W] = wirePlanField_vector(r_wire, M, B_0, x_0, y);
mu_0 = mu_0 / norm(mu_0);

%Tilt de mu_0 (dépolarisation partielle du fasceau)
% si H // x
%mu_0 = [1,0,0];
angleEcart = 20; % en °
%anglePhase = 0; % en °

%mu_0 = [cos(angleEcart*pi/180), 0, sin(angleEcart*pi/180)]; % selon x
%mu_0 = [0, cos(angleEcart*pi/180), sin(angleEcart*pi/180)]; % selon y
%mu_0 = [-sin(angleEcart*pi/180), 0, cos(angleEcart*pi/180)]; % selon z
%mu_0 = rotGeneral(mu_0,[0,1,0],angPhase*pi/180);



% 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, MagMod),[x_0:r_wire/100:x_f],mu_0,OPTIONS);
%[x,m] = ode113(@(x, mu) precessXY(x,mu, y, r_wire, M, B_0, neutron_speed, gamma_n, MagMod),[x_0:r_wire/100:x_f],mu_0,OPTIONS);






% visualisation des data (à mettre à jour !!!!)
mx = m(:,1);
my = m(:,2);
mz = m(:,3);

if plotYN==1
    figure(1)
    hold off
    plot(x,mx,'r')
    hold on
    plot(x,my,'g')
    plot(x,mz,'b')
    legend('Sx','Sy','Sz');
    %%%%%%%%%%%%%%
    figure(2)
    clf
    [X,Y,Z]=sphere(20);
    mesh(X,Y,Z)
    hidden off
    hold on
    plot3(mx,my,mz,'.','MarkerSize',20) %nonuniform
    xlabel('Sx');ylabel('Sy');zlabel('Sz');
    %%%%%%%%%%%%%%%
    figure(3)
    clf
    mesh(0.01*X,0.01*Y,0.01*Z);
    aux=size(mx);
    for i=1:aux(1),
        %[i 0 0]
        %[i+mx(i) my(i) mz(i)]
        h=line([x(i) x(i)+mx(i)],[0 my(i)],[0 mz(i)]);
        if (i<(aux(1)/4)) | (i>(aux(1)*3/4))
            set(h,'Color','b');
        else 
            set(h,'Color','r');
        end
    end
    axis([0 x(aux(1))*neutron_speed -1 1 -1 1])
    xlabel('Sx');ylabel('Sy');zlabel('Sz');
end

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