classdef C_BeamPolar<handle
%%% bp = C_BeamPolar(polarNeutron, Y, M, B, wavelength, wlValues, wlRatio, r, MagMod)
% 
% Object contenant la polarisation des neutrons en fonction de l'espace, et
% de la longeur d'onde. Les moyennes sur la largeur du fasceau et sur la
% longeur d'onde sont réalisées à postériori.
% La classe possède des méthodes pour tracer facilement la polarisation et
% le flipping ratio (voir doc des méthodes).
%
% Methode  de création de l'object N° 1 :
% ======================================
% En donnant tous les arguments
%
% Arguments obligatoires:
% -----------------------
%
% - polarNeutron: matrice de dimension 2 contenant la polarisation en
%   fonction de y selon une dimension et en fonction de la longueur d'onde
%   selon l'autre.
% - Y : vecteur positions en y
% - M : vecteur xyz d'aimantation
% - B : vecteur xyz de champ appliqué
% - wavelength : longeur d'onde centrale
% - wlValues : vecteur des longeurs d'onde utilisées
% - wlRatio : vecteur de pondération des longeurs d'onde utilisées
% - r rayon du fil
%
%
% Methode  de création de l'object N° 2 :
% ======================================
% Par chargement d'un fichier *.mat créé par la fonction computePolarization
% bp = C_BeamPolar(cheminDuFichier)
%
%
%
%
% Liste des Méthodes publiques : 
% ==============================
%
% Calcule de la polarisation ou du flipping ratio en fonction de y en moyennant
% sur les longeurs d'onde et sur la largeur du fasceau
% -> [Y_new, P_mean] = getPolar(self, beamWidth)
% -> [Y_new,R] = getFlippingRatio(self, beamWidth)
%
% Traçage de la polarisation ou du flipping ratio calculé par les fonctions 
% ci-dessus
% -> [h] = plotFlippingRatio(beamWidth, h, ST_legend)
% -> [h] = plotPolar(self, beamWidth, h, ST_legend)
%
% Changement de la dispestion en longeur d'onde
% -> resetGaussianDisp(sigma) 
%
%
% $Revision: 15 $
% $Author: lenertz $
% $Date: 2015-07-30 17:36:07 +0200 (jeu., 30 juil. 2015) $
% contact: lenertz@cea.fr
 
    % Properties definition
    % ---------------------
    properties (SetAccess = private, GetAccess = public)
        polarData = {};
        %Y = {};
        %y_step = {};
        %largeur = {};
        M = {};
        B = {};
        wavelength = {};
        wlValues =  {};
        r = {};
        wlRatio = {};
        LI_colors = lines(20);
        MagMod = {};
    end
    properties (SetAccess = public, GetAccess = public)
        Y = {};
        y_step = {};
        largeur = {};
    end

    % Class public methods
    % --------------------
    methods 
    
        % Constructor
        % -----------
        function self = C_BeamPolar(polarNeutron, Y,  M, B, wavelength, wlValues, wlRatio, r, MagMod)
          
          if nargin == 8
              MagMod = 0;
          elseif nargin == 2
              MagMod = Y;
          end
              
          if nargin < 3 
              load(polarNeutron)
              wlRatio = WL_RATIO;
              polarNeutron = cos(polarNeutron);
              B = H;
              wlValues = WL;
          
              
          end
            
          self.polarData = polarNeutron;
          self.Y = Y;
          self.largeur = (Y(end) - Y(1))/2;
          if numel(Y)==1;
              self.y_step = 0;
          else
              self.y_step = 2*self.largeur/(numel(Y)-1);
          end 
          
          self.M = M;
          self.B = B;
          self.wavelength = wavelength;

          if numel(wlRatio) == 1 && wlRatio ~= 1;
              self.wlRatio(1:numel(polarNeutron(1,:))) = 1/numel(polarNeutron(1,:));
          else
              self.wlRatio = wlRatio;
          end
          
          self.wlValues = wlValues;
          
          self.r = r;
          self.MagMod = MagMod;

          end

        % Methods
        % -------      
        function [Y_new, P_mean] = getPolar(self, beamWidth)
            %% Renvoie la polarisation
            % [Y,P] = getPolar(beamWidth)
            %
            % Calcule de ma polarisation P en sortie en fonction de Y et et
            % de la largeur du faisceau. Renvoie Y er P.

            [meanPolarNeutron] = self.wlMean();
            [Y_new, P_mean] = self.beamWidthMean(meanPolarNeutron, beamWidth);
        end
        
        function [h] = plotPolar(self, beamWidth, h, ST_legend)
            %% Trace la polarisation
            % [h] = plotPolar(beamWidth, h, ST_legend)
            % 
            % Trace la polarisation en sortie en fonction de y (largeur) et
            % de la largeur du faisceau et renvois la référence de la
            % figure.
            %
            % Arguments:
            % ----------
            % - beamWidth: largeur du faisceau (m) sur laquelle le signal
            %              est moyenné 
            % 
            % Arguments optionnels:
            % ---------------------
            % - h: Référence vers une figure déjà existante. Si non donné,
            %      une nouvelle figure est tracé
            % - ST_legend: chaîne de caractre remplaçant la valeur par
            %              defaut

            if nargin < 3; h = 0; end
            if nargin < 4; ST_legend = 0; end
            [meanPolarNeutron] = self.wlMean();
            [Y_new, P_mean] = self.beamWidthMean(meanPolarNeutron, beamWidth);
            [h] = self.plotData(Y_new, P_mean, 'P', beamWidth, h, ST_legend);
            ylim([-1.1,1.1])
        end
        
        function [Y_new,R] = getFlippingRatio(self, beamWidth)
            %% Renvoie le flipping ratio
            % [Y,R] = getFlippingRatio(beamWidth)
            %
            % Calcule le flipping ratio R en sortie en fonction de Y et et
            % de la largeur du faisceau. Renvoie Y er R.
            
            [meanPolarNeutron] = self.wlMean();
            [Y_new, P_mean] = self.beamWidthMean(meanPolarNeutron, beamWidth);
            R = (P_mean+1)./(1-P_mean);
            %Y_new = -Y_new;
        end   
        
        function [h] = plotFlippingRatio(self, beamWidth, h, ST_legend)
            %% Trace le flipping ratio
            % [h] = plotFlippingRatio(beamWidth, h, ST_legend)
            % 
            % Trace le flipping ratio en sortie en fonction de y (largeur) 
            % et de la largeur du faisceau et renvois la référence de la
            % figure.
            %
            % Arguments:
            % ----------
            % - beamWidth: largeur du faisceau (m) sur laquelle le signal
            %              est moyenné 
            % 
            % Arguments optionnels:
            % ---------------------
            % - h: Référence vers une figure déjà existante. Si non donné,
            %      une nouvelle figure est tracé
            % - ST_legend: chaîne de caractre remplaçant la valeur par
            %              defaut
            
            if nargin < 3; h = 0; end
            if nargin < 4; ST_legend = 0; end
            [Y_new,R] = getFlippingRatio(self, beamWidth);
            [h] = self.plotData(Y_new, R, 'flipping ratio', beamWidth, h, ST_legend);
        end
        
        function resetGaussianDisp(self, sigma)
            % resetGaussianDisp(sigma) 
            % 
            % Redéfinie les ratio de longeur d'onde en fonction de la
            % nouvelle dispertion de la gaussienne sigma donnée en argument
            
            self.wlRatio = gaussian(self.wlValues,sigma,self.wavelength);
            self.wlRatio = self.wlRatio / sum(self.wlRatio);
            
        end
    end   
    methods(Access = protected) % Méthode en cours de développement, non-utilisables en l'état
        
        
        function [finalMean, new_y] = getPolarVersusLambda(self, beamWidth, lambdaGWidth, lambdaSelectWidth)
            %% 
            %  TODO Finir fonction
            % 
            % 
            
            
            
            %% moyenne sur la largeur du faisceau
            % on détermine les index des positions en y a moyenner
            y_nbStep = round(beamWidth/self.y_step);
            if mod(numel(self.Y),2) == 1;
                y_middleIndex = numel(self.Y)/2 +0.5;
                if mod(y_nbStep, 2) == 1;
                   y = y_middleIndex -(y_nbStep-1)/2:y_middleIndex +(y_nbStep-1)/2;
                else
                    y = y_middleIndex - y_nbStep/2:y_middleIndex + y_nbStep/2;
                end
            else
                if mod(y_nbStep, 2) == 1;
                   y = numel(self.Y)/2 -(y_nbStep-1)/2:  numel(self.Y)/2 + 1 + (y_nbStep-1)/2; 
                else
                   y = numel(self.Y)/2 -y_nbStep/2 + 1:  numel(self.Y)/2 +y_nbStep/2; 
                end
            end
            if numel(y) - numel(self.Y) == -2 && numel(y) > 2; 
                y = y(2:end-1);
            end
            % moyenne
            meanPolarNeutron = sum(self.polarData(y,:), 1);
            
            
            
            
           %% Moyennes sur la dispertion en longeurs d'onde
           wl_nbStep = round(lambdaSelectWidth/abs(self.wlValues(2) - self.wlValues(1)));
           if wl_nbStep < 1
               finalMean = meanPolarNeutron;
               new_y = self.Y;
           else
               finalMean(1:numel(self.wlValues)- wl_nbStep +1) = 0;
               % Calcul de la gaussienne
               wlRatioLocal = gaussian(self.wlValues(1:wl_nbStep),lambdaGWidth,mean(self.wlValues(1:wl_nbStep)));
               for i = 1:numel(wlRatioLocal)
                   finalMean = finalMean + meanPolarNeutron(i:end-wl_nbStep+i) .* wlRatioLocal(i);
                   
               end
               finalMean = finalMean / wl_nbStep;
               new_y = (self.Y(1:end-wl_nbStep+1) + self.Y(wl_nbStep:end)) / 2; 
              
           end    

        end
    
        function [h] = plotPolarVersusLambda(self, beamWidth, lambdaGWidth, lambdaSelectWidth, h, ST_legend)
        %%
            if nargin < 3; h = 0; end
            if nargin < 4; ST_legend = 0; end
            [Y_new, P_mean] = self.getPolarVersusLambda(beamWidth, lambdaGWidth, lambdaSelectWidth);
            [h] = self.plotData(Y_new, P_mean, 'P', beamWidth, h, ST_legend);
            ylim([-1.1,1.1])
        
        end
    
    
    %get access methods for properties 
    %---------------------------------
%     function test = get.test(self)
%       test = self.test + 1;
%     end
%     
%     function data = get.data(self)
%       data = self.data;
%     end
    
    end  % end of public methods
    
    %% Protected Methods
    %  -----------------
    methods(Access = protected)
    
        function [h] = plotData(self, X, Y, y_name, beamWidth, h, ST_legend)
            %% Trace la moyenne sur la largeur du faisceau
           
            if nargin < 6 || h == 0;
                h = figure;
                ST_title = strcat(['B=(',...
                              num2str(self.B(1)),',',...
                              num2str(self.B(2)),',',...
                              num2str(self.B(3)),') M=(',...
                              num2str(self.M(1)),',',...
                              num2str(self.M(2)),',',...
                              num2str(self.M(3)),') MagMod=',...
                              num2str(self.MagMod), ' wl=',...
                              num2str(mean(self.wavelength)),'m']);
                title(ST_title)
                
                xlabel('y (m)','FontSize',16); 
                ylabel(y_name,'FontSize',16);
                
                colorIndex = 1;
            else
                figure(h);
                ax = get(h,'CurrentAxes');
                colorIndex = numel(get(ax,'children')) +1;
                if colorIndex > numel(self.LI_colors(:,1))
                    colorIndex = 1;
                end
                
            end
              
            
            hold on
            plot(X,Y, 'color', self.LI_colors(colorIndex,:));
            legend('-DynamicLegend')
            legend('show')
            ax = get(h,'CurrentAxes');
            ax = get(ax,'children');
            if nargin < 7 || strcmp(num2str(ST_legend),'0');
                set(ax(1),'DisplayName',strcat(['beamWidth = ',num2str(beamWidth*1e6),' um']))
            else
                set(ax(1),'DisplayName',ST_legend)
            end
        end
        
        
        
        function [Y_new, P_mean] = beamWidthMean(self, meanPolarNeutron, beamWidth)
            %% Renvoie la polarisation moyennée par la largeur du fasceau
            % [Y,P] = beamWidthMean(beamWidth)
            %
            % Calcule de ma polarisation R en sortie en fonction de Y, de
            % la largeur du faisceau et de la longeur d'onde (voir 
            % self.wlMean). Renvoie Y er R.
            
   
            
            index = round(beamWidth/self.y_step);
             
            if index < 1
                index = 1;    
            elseif index > numel(meanPolarNeutron)
                index = numel(meanPolarNeutron);
            end
            disp(strcat(['Beam integration on ', num2str(index), ' points']))

                P_mean = meanPolarNeutron(1:end - index+1);
                for j = 1:(index-1);
                    P_mean = P_mean + meanPolarNeutron(index-j +1 :end - j+1);
                end;

                % plot(Y(round(i/2):end - round(1+ i/2)+1),cos(meanPolarNeutron2/i));
                Y_new = self.Y(round(index/2): round(index/2) + numel(P_mean) -1);
                P_mean = P_mean/index;
            
            
            %[Y_new, P_mean] = beamWidthMean(self.Y,meanPolarNeutron,beamWidth,BO_plot);        
        end
    
    
    
        function [meanPolarNeutron] = wlMean(self)
            %% Renvoie la polarisation moyennée sur la longeur d'onde
            % [meanPolarNeutron] = wlMean()
            %
            % Renvoie la moyenne de la polarisation par la longeur d'onde 
            % pondéré par les ratios définie en interne dans self.wlRatio
            % pour chaque longeur d'onde. Par contre la polarisation n'est
            % moyennée par la largeur du faisceau.
            
            
            meanPolarNeutron = self.polarData;
            for wl_index = 1:numel(self.wlRatio);
                meanPolarNeutron(:,wl_index) = meanPolarNeutron(:,wl_index)*self.wlRatio(wl_index);
            end   
            meanPolarNeutron = sum(meanPolarNeutron, 2);
        end  
	end
       

end  % end of class definition   
    
