
%%% Evaluation de la polarisation du faisceau au travers d'un fil

%% Paramètre du système %%
magFunc = @wirePlanField_vector

% largeur de mesure 
largeur = 300e-6;
% pas d'intégration
y_step = 5e-6;
% rayon du fil
r = 100e-6; % (m)
% longeur d'onde 
wavelength = 4e-10; %(m)
% dispertion en longueur d'onde (largeur à mis hauteur de la gaussienne)
wl_sigma = wavelength*0.07; %(m)
% pas de mesure en longeur d'onde
wl_step = wl_sigma;
% Sur conbien de sigma on mesure la divergence en longeur d'onde
wl_mesure = 2;
% largeur du faiseau
beam_width = 100e-6; %(m)

% Magnetic field
H = [0.00,0.010,0.000];
M_sat = 0.6; %(Tesla)
M = M_sat*H/norm(H); 


disp(strcat(['largeur = ', num2str(largeur),' m']))
disp(strcat(['y_step = ', num2str(y_step),' m']))
disp(strcat(['r = ', num2str(r),' m']))
disp(strcat(['wavelength = ', num2str(wavelength),' m']))
disp(strcat(['wl_sigma = ', num2str(wl_sigma),' m']))
disp(strcat(['wl_step = ', num2str(wl_step),' m']))
disp(strcat(['wl_mesure = ', num2str(wl_mesure),' sigma']))
disp(strcat(['beam_width = ', num2str(beam_width),' m']))
disp(strcat(['H = ', num2str(H),' T']))
disp(strcat(['M_sat = ', num2str(M_sat),' T']))


%% fin des paramètres variables %%

% vecteur positions latérales
y_start = -largeur;
y_end = largeur;
Y = y_start:y_step:y_end;


% vecteur longeurs d'onde
WL = wavelength-wl_mesure*wl_sigma:wl_step:wavelength+wl_mesure*wl_sigma;
% ratio d'intensité (gaussien) en fonction de la longeur d'onde
WL_RATIO = gaussian(WL,wl_sigma,wavelength);
WL_RATIO = WL_RATIO / sum(WL_RATIO);
% nombre de pas en longueur d'onde
wl_total_index = numel(WL_RATIO);
% nombre de pas total
total_index = numel(Y)*wl_total_index;

% tableau de résultats
clearvars polarNeutron
polarNeutron(1:numel(Y),1:wl_total_index)=0;


% Calcul de la polarisation pour chaque position lateral et chaque longeur
% d'onde
for y = Y;  
    for wl_index = 1:numel(WL);
        % longeur d'onde
        wl = WL(wl_index);
        %out = travel2(M, H,r,y,wl,0);
        out = travel3(M, H,r, X, y, wl, magFunc)
        % index de la position lateral dans le vecteur résultat
        index = round((y-y_start)/y_step + 1);
        % affiche l'avancement
        disp(strcat(['step ', num2str((index-1)*wl_total_index + wl_index), '/', num2str(total_index)])) % affiche la progression

        [U,V,W] = wirePlanField_vector(r, M, H, out(:,1), y);
        B = [U, V, W];
        theta = vectorAngles(B, out(:,2:4));
        polarNeutron(index,wl_index) = theta(end);
    end
end



% moyenne la precession sur la largeur du faisceau et sur la longeur d'onde
meanPolarNeutron1 = polarNeutron;
for wl_index = 1:numel(WL);
    meanPolarNeutron1(:,wl_index) = meanPolarNeutron1(:,wl_index)*WL_RATIO(wl_index)
end   
meanPolarNeutron1 = sum(meanPolarNeutron1, 2);

i = round(beam_width/y_step);
disp(strcat(['Beam integration on ', num2str(i), ' points']))
if i > 0;
    meanPolarNeutron2 = meanPolarNeutron1(1:end - i+1);
    for j = 1:i-1;
        meanPolarNeutron2 = meanPolarNeutron2 + meanPolarNeutron1(i-j +1 :end - j+1);
    end;
    figure;
    plot(Y(round(i/2):end - round(1+ i/2)+1),cos(meanPolarNeutron2/i));
    ylim([-1.1,1.1])
end

% u = out(:,2);
% v = out(:,3);
% w = out(:,4);
% x = out(:,1);
% y = x; y(:) =0;
% z = y;
% 
% figure;
% a=find(x > -2e-4 & x < 10e-4);
% quiver3(x(a),y(a),z(a),u(a),v(a),w(a)); 
% daspect([1 1 1]);
% 
% 
% [U,V,W] = wirePlanField_vector(1e-4, [1,0,0], [0.01,0,0], out(:,1), y);
% 
% B = [U, V, W];
% theta = vectorAngles(B, out(:,2:4));
% figure;
% plot(x,theta*180/pi)
