function [polarNeutronObject] = computePolarization(H,M,MagMod,Params, supMagParams, startPolar, angPhase, hp, hf)
%% Calcule la polarisation du faisecau de neutrons
% computePolarization(H,MagMod,Params)
%
% Evaluation de la polarisation du faisceau au travers d'un fil. Le faiseau
% de neutron se déplace en x et fil est tiré le long de z
%
% Paramètres:
% -------------
%
% H=[Hx,Hy,Hz] : Champ magnétique appliqué
%
% M: Aimantation à saturation du fil (en [T])
%   - Si M est un vecteur, sa norm sera M_sat.
%   - Si M est un scalaire, il représentera M_sat et l'aimantation sera 
%     dans la même direction que H.
%   - Si M = 0, l'aimantation prend la valeur par défaut M_sat = 0.6T et M
%     sera dans la même direction que H
%
% MagMod=1: modèle magnétique d'un aimantation uniforme dans le film
%       =2: Aimantation sous forme de vortex
%       = @magFunction: où magFunc est une fonction de r_wire, M, B_0, X, Y
%       et supMagParams (même s'il ne sont pas utilisés dans la fonction)
%                            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
%                            supMagParams -> liste te tous les autres paramètres
%                                            nécessaires
%                       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.
%

% Paramètres Optionnels:
% ----------------------
% Params = [longeur, x_step, largeur, y_step, r, wavelength, wl_sigma, wl_nbStep,...
%           wl_mesure, beam_width, BO_save, angPhase]
% (paramètre (valeur par defaut): explication)
% Ces paramètres doivent être données dans l'ordre ci dessus. Un paramètre
% déclaré avec -1 sera remplacé par sa valeur par defaut sauf pour les 
% derniers paramètres (startPolar, angPhase). Si la longeur de
% 'Params' est inférieur au nombre de paramètres, les paramètres non
% déclarés seront également remplacé par leur valeur par defaut.
% 
% - longueur (50*r): +- longueur de mesure [en m]
% - x_step (longueur*1e-4): pas de mesure en x [en m]
% - largeur (500e-6 m): +-largeur de mesure latéral, selon y [en m]
% - y_step (5e-6) : pas de mesure latéral, selon y [en m]
% - r (100e-6) : rayon du fil [en m]
% - wavelength (4e-10) : longeur d'onde [en m]
% - wl_sigma (wavelength*0.07): dispersion en longueur d'onde, ecart-type de
%                               la gaussienne (la largeur à mi-hauteur de la
%                               gaussienne vaut environ 2,35 sigma) [en m]
% - wl_nbStep (3): nombre de pas en longeur d'onde pour la convolution [int>0]
% - wl_mesure (0): +- largeur de mesure en longeur d'onde en unité de 
%                  divergence (largeur total = 2*wl_mesure). 
% - beam_width (100e-6) : largeur du faiseau [en m]. Ne sert à la fin qu'à
%                         tracer le résultat. L'object polarNeutronObject
%                         peremt de changer ce paramètre à postériori.
% - BO_save (false/0) : Booléen décriavant si on sauvegarde ou non les
%                       données à la fin dans une fichier par défaut
%                       (voir variable ST_folderName pour le dossier par
%                        defaut).
%
% supMagParams (0): liste de tous les autres paramètres nécessaires au
%                   calcule du champ et qui sera passé en dernier argument
%                   à la fonction de calcule du champ. 
% startPolar (1): Polarisation initiale du faisceau(comprise entre -1 et 1)
% 
% angPhase (0) : phase de départ [en rad] pour une polarisation de départ 
%                différente de 1 ou -1 (paramètre utilisé dans travel3). Si 
%                angPhase est une liste de nombre le resultats sera une 
%                moyenne des différentes phases. 
% 
% hp: référence à une figure pour y ajouter la coube de polarisation en
%     fonction de y. Si hp = -1, aucune figure n'est tracés.
%
% hp: référence à une figure pour y ajouter la coube de rapport de flipping
%     en fonction de y Si hp = -1, aucune figure n'est tracés.
%
%
% $Revision: 14 $
% $Author: lenertz $
% $Date: 2015-07-27 18:22:01 +0200 (lun., 27 juil. 2015) $

    


    % ajoute les sous dossiers
    addSubFolders();

    

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %% Initialisation des paramètres du système %%
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    fprintf('------------------\n|   Parameters   |\n------------------\n')
    
    % mode de calcul du champ magnétique
    if isa(MagMod, 'function_handle') || isa(MagMod, 'char');
        magFunc = MagMod;
        MagMod = 0;
    elseif MagMod == 1;
    	magFunc = @wirePlanField_vector;
    elseif MagMod == 2;
    	magFunc = @wireVortexField_vector;
    elseif MagMod == 3;
    	magFunc = @wire3Field;
    end
    disp(strcat(['Magnetic function: ', func2str(magFunc)]))
    
    

    
    % On complète la la longeur du paramètre 'Params' avec les 0
    if nargin < 4
        Params(1:12) = -1;
    elseif numel(Params) <  11;
        Params(numel(Params)+1:11) = -1;
    end
    
    
   
    % largeur de mesure, selon x (en m)
    if Params(3) == -1;
        largeur = 500e-6;
    else
        largeur = Params(3);
    end
    disp(strcat(['largeur = ', num2str(largeur),' m']))

    % pas d'intégration latéral, selon y(en m)
    if Params(4) == -1;
        y_step = 5e-6;
    else
        y_step = Params(4);
    end
    disp(strcat(['y_step = ', num2str(y_step),' m']))

    % rayon du fil
    if Params(5) == -1;
        r = 100e-6; % (m)
    else
        r = Params(5);
    end
    disp(strcat(['r = ', num2str(r),' m']))

    
    % longueur de mesure, selon x (en m)
    if Params(1) == -1;
        longueur = 50*r;
    else
        longueur = Params(1);
    end
    disp(strcat(['longeur = ', num2str(longueur),' m']))
    
    % pas d'intégration, selon x (en m)
    if Params(2) == -1;
        x_step = longueur*1e-4;
    else
        x_step = Params(2);
    end
    disp(strcat(['x_step = ', num2str(longueur),' m']))
    
    
    % longeur d'onde 
    if Params(6) == -1;
        wavelength = 4e-10; %(m)
    else
        wavelength = Params(6);
    end
    disp(strcat(['wavelength = ', num2str(wavelength),' m']))

    % dispertion en longueur d'onde (largeur à mi-hauteur de la gaussienne)
    if Params(7) == -1;
        wl_sigma = wavelength*0.07; %(m)
    else
        wl_sigma = Params(7);
    end
    disp(strcat(['wl_sigma = ', num2str(wl_sigma),' m']))

    % pas de mesure en longeur d'onde
    if Params(8) == -1;
        wl_nbStep = wl_sigma;
    else
        wl_nbStep = Params(8);
    end
    wl_nbStep = round(abs(wl_nbStep));
    disp(strcat(['wl_nbStep = ', num2str(wl_nbStep)]))

    % Sur conbien de sigma on mesure la divergence en longeur d'onde
    if Params(9) == -1;
        wl_mesure = 0;
    else
        wl_mesure = Params(9);
    end
    disp(strcat(['wl_mesure = ', num2str(wl_mesure),' sigma']))

    % largeur du faiseau
    if Params(10) == -1;
        beam_width = 100e-6; %(m)
    else
        beam_width = Params(10);
    end
    disp(strcat(['beam_width = ', num2str(beam_width),' m']))

    % sauvegarde à la fin
    if Params(11) == -1;
        BO_save = false;
    else
        BO_save = Params(11);
    end
    
    if nargin < 5;
       supMagParams = 0; 
    end
    disp(strcat(['Supplementary mag conditions: ', num2str(supMagParams)]))
    
    if nargin < 6;
       startPolar = 1; 
    end
    disp(strcat(['start polarization = ', num2str(startPolar)]))
    
    
    if nargin < 7 || startPolar==1 || startPolar==-1;
       angPhase = 0; 
    end
    if not(startPolar == 1 || startPolar == -1);
        disp(strcat(['Polarization phase = ', num2str(angPhase), ' rad']))
    end
    
    
    % Magnetic field
    %H = [0.00,0.010,0.000];
    if numel(M) == 1;
        if M == 0;
            M_sat = 0.6; %(Tesla)
        else 
            M_sat = M;
        M = M_sat*H/norm(H);
        end
    else
        M_sat = norm(M);
    end
    disp(strcat(['H = ', num2str(H),' T']))
    disp(strcat(['(H_sat = ',num2str(norm(H)),')' ]))
    disp(strcat(['M = ', num2str(M),' T']))
    disp(strcat(['(M_sat = ',num2str(M_sat),')' ]))
    fprintf('\n') 
    
     
    
    
    %% fin des paramètres variables %%

    % vecteur position longitudinale
    x_start = -longueur;
    x_end = longueur;
    X = x_start:x_step:x_end;
    
    % vecteur positions latérale
    y_start = -largeur;
    y_end = largeur;
    Y = y_start:y_step:y_end;


    % vecteur longeurs d'onde
    if wl_nbStep == 1;
        WL = wavelength;
    else
        WL = wavelength-wl_mesure*wl_sigma:2*wl_mesure*wl_sigma/(wl_nbStep-1):wavelength+wl_mesure*wl_sigma;
    end
    disp(strcat(['wavelength vector = ', num2str(WL),' m']))
    % 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);
    disp(strcat(['wavelength ratios = ', num2str(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;
    polarNeutronSum = polarNeutron;

     
    
    %% Calcul de la polarisation pour chaque position laterale et chaque 
    %  longeur d'onde
    
    for phaseIndex = 1:numel(angPhase);
        if numel(angPhase) > 1
            disp(strcat(['Phase computation ',...
                         num2str(phaseIndex),...
                         '/',...
                         num2str(numel(angPhase))]))
        end
        tic;
        phaseInLoop = angPhase(phaseIndex);
        progress = 0.05;
        for index = 1:numel(Y);
            y = Y(index);
            parfor wl_index = 1:numel(WL); % for wl_index = 1:numel(WL);
                % longeur d'onde
                wl = WL(wl_index);
                out = travel3(M, H, r, X, y, wl, magFunc, startPolar, phaseInLoop, supMagParams);

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

            if progress < index/numel(Y);
                % affiche la progression
                duration = toc;
                disp(strcat(['progress: ',...
                             num2str(progress*100),...
                             '% step ', ...
                             num2str(index),...
                             '/', num2str(numel(Y)),...
                             ' remaining time: ',...
                             num2str(floor(duration*(1-progress)/(60*progress))),...
                             'min',...
                             num2str(round(mod(duration*(1-progress)/progress,60))),...
                             's']))
                 progress = progress + 0.05;
            end
        end
        polarNeutronSum = polarNeutronSum + polarNeutron;
    end
    polarNeutron = polarNeutronSum / numel(angPhase);
    fprintf('Done !\n\n')

    %% Empaquetage des données et traçage des figures
     
    polarNeutronObject = C_BeamPolar(cos(polarNeutron), Y,  M, H, wavelength, WL, WL_RATIO, r, MagMod);
    
    % Plot polarization and flipping Ratio
    if nargin < 8
        polarNeutronObject.plotPolar(beam_width);
    else
        if hp ~= -1;
            polarNeutronObject.plotPolar(beam_width, hp);
        end
    end
    if nargin < 9
        polarNeutronObject.plotFlippingRatio(beam_width);
    else
        if hf ~= -1;
            polarNeutronObject.plotFlippingRatio(beam_width, hf);
        end
    end
    
    
    %% Sauvegarde des données
    if BO_save == 1;
        
        ST_folderName = 'D:\Data simul mag\magneticWire\';
        ST_subFolderName = strcat(['B=(',...
                                  num2str(H(1)),',',...
                                  num2str(H(2)),',',...
                                  num2str(H(3)),') M=(',...
                                  num2str(M(1)),',',...
                                  num2str(M(2)),',',...
                                  num2str(M(3)),') MagMod=',...
                                  num2str(MagMod)]);
        ST_saveFolder = strcat(ST_folderName,ST_subFolderName,'\');
        
        if not(exist(ST_subFolderName,'dir'))
            mkdir(ST_saveFolder)
        end

        fileName = strcat(ST_saveFolder,'Data','.mat')
        save(fileName,...
            'Y',...
            'polarNeutron',...
            'WL',...
            'WL_RATIO',...
            'largeur',...
            'y_step',...
            'r',...
            'wavelength',...
            'wl_sigma',...
            'wl_nbStep',...
            'wl_mesure',...
            'beam_width',...
            'M_sat',...
            'H',...
            'M',...
            'MagMod')
    
    end


