function [pol] = Compute_Spin_Polarization(B,t,S0,varargin)
%function [pol] = Compute_Spin_Polarization(B,tRange,S0,<Parameters>  )
%
%Computes MuSR spin polarization as a function of time
%
%XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
%X Author:          Sky Cheung                                                X
%X Date Created:    19 Apr 2016                                               X
%X Last Modified:   15 Jun 2016                                               X
%XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
%

%Input Parsing
p = inputParser;

%Parsing conditions
addRequired(p, 'B', @(x) validateattributes( x, ...
    {'numeric'}, {'nonempty','2d','ncols',3} ) );
addRequired(p, 't', @(x) validateattributes( x, ...
    {'numeric'}, {'nonempty','vector'} ) );
addRequired(p, 'S0', @(x) validateattributes( x, ...
    {'numeric'}, {'nonempty','2d','ncols',3} ) );

%Parsing Input
parse(p, B, t, S0, varargin{:});

%Retrieving Proper Inputs
B = p.Results.B;
t = p.Results.t;
S0 = p.Results.S0;

%Deleting Parser
clear p;

%Preallocaitng
numField = size(B,1);
numTimes = length(t);
%forcing times to be a row
t = t(:)';

%S = zeros(numTimes,3,numField);
% %For loop over all fields
% parfor j = 1:numFields
%     [tempT,tempSolve] = ode23( ...
%         @(t,spin) SpinEOM(t, B(j,:), spin ), times, S0, ODEOption);
%     tempS = interp1(tempT,tempSolve,times);
%     S(:,:,j) = tempS;
% end
%Compute ZF
%pol = sum(S,3)/numField;

%Time Dependence

%Input check
if isequal(size(B), size(S0))
    S0Mat = S0;
elseif isvector(S0)
    %Making S0Mat with S0
    S0Mat = repmat(S0,numField,1);
else
    error('Check B and S0');
end

BUnitMat = normr(B);
BMag = sqrt(sum(B,2).^2);
S0DotBU = repmat( dot( S0Mat, BUnitMat, 2 ), 1, numTimes); %Dot product by row
S0CrossBU = cross( S0Mat, BUnitMat, 2); %Cross Product by row
muonGyro = 13500; %Hz/G

scaledTimeMat = BMag * muonGyro * t;
S = cell(1,3);
for j = 1:3
    S{j} = bsxfun(@times, BUnitMat(:,j), S0DotBU ) + ...
        bsxfun(@times, ...
            (bsxfun(@minus, S0Mat(:,j), ...
                bsxfun(@times, BUnitMat(:,j), S0DotBU ) ...
                ) ...
            ), cos(scaledTimeMat) ...
        ) + ...
        bsxfun(@times, S0CrossBU(:,j), sin(scaledTimeMat) );
end
%Compute ZF
pol = cell2mat( cellfun(@sum, S, 'UniformOutput',0)' )' / numField;

