function [totalField] = Compute_Field_NaFeCuAs(muonLoc,impCon,numMC, ...
    magLatticePos,magLatticeMom,nucLatticePos,nucLatticeMom)
%function [totalField] = Compute_Field_NaFeCuAs(muonLoc,impCon,numMC)
%
%Computes local magnetic field at the muon site in NaFeCuAs
%
%
%INPUTS
%
%muonLoc                        Coordinates of muon location 
%                                   fractions of a,b,c values
%
%impCon                         Impurity Concentration
%                                   Scalar between 0 and 1
%
%numMC                          Number of indpenedent trials
%
%magLatticePos                  Matrix of magnetic positions
%
%magLatticeMom                  Matrix of magnetic moments
%
%nucLatticePos                  Matrix of nuclear positions
%
%nucLatticeMom                  Matrix of nuclear moments
%
%
%OUTPUTS
%
%totalField                     N x 3 matrix of the local field
%
%
%XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
%X Author:          Sky Cheung                                                X
%X Date Created:    19 Apr 2016                                               X
%X Last Modified:   19 Apr 2016                                               X
%XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
%

%Input Parsing
p = inputParser;

%Parsing conditions
addRequired(p, 'muonLoc', @(x) validateattributes( x, ...
    {'numeric'}, {'nonempty','vector','numel',3} ) );
addRequired(p, 'impCon', @(x) validateattributes( x, ...
    {'numeric'}, {'nonempty','scalar','<=',1,'>=',0} ) );
addRequired(p, 'numMC', @(x) validateattributes( x, ...
    {'numeric'}, {'nonempty','positive','integer'} ) );
addRequired(p, 'magLatticePos', @(x) validateattributes( x, ...
    {'numeric'}, {'nonempty','2d','ncols',3} ) );
addRequired(p, 'magLatticeMom', @(x) validateattributes( x, ...
    {'numeric'}, {'nonempty','2d','ncols',3} ) );
addRequired(p, 'nucLatticePos', @(x) validateattributes( x, ...
    {'numeric'}, {'nonempty','2d','ncols',3} ) );
addRequired(p, 'nucLatticeMom', @(x) validateattributes( x, ...
    {'numeric'}, {'nonempty','2d','ncols',3} ) );

%Parsing Input
parse(p, muonLoc, impCon, numMC, magLatticePos, ...
    magLatticeMom, nucLatticePos, nucLatticeMom);

%Retrieving Proper Inputs
muonLoc = p.Results.muonLoc;
impCon = p.Results.impCon;
numMC = p.Results.numMC;
magLatticePos = p.Results.magLatticePos;
magLatticeMom = p.Results.magLatticeMom;
nucLatticePos = p.Results.nucLatticePos;
nucLatticeMom = p.Results.nucLatticeMom;

%Deleting Parser
clear p;

%Material constants
%Import Unit Cell Information (Angstrom to cm)
a = 3.94940 * 10^-8;
b = 3.94940 * 10^-8;
c = 7.03960 * 10^-8;
%Cu Nuclear
nucMagnetonCu = 2.2226;
nucMagneton = 5.051 * 10^(-24);

%Preallocate for loop over MC
totalField = zeros(numMC,3);

%Determining number of impuriites
numFeAtoms = size(magLatticePos,1);
numImp = floor( impCon * numFeAtoms );

%Muon site
muonSite = muonLoc .* [a,b,c];

%Impurity For loop
parfor j = 1:numMC
    %Non-magnetic impurity
    impurityIdx = randperm(numFeAtoms, numImp);
    
    %MC Moment Lattice
    magLatticePosTemp = magLatticePos;
    magLatticeMomTemp = magLatticeMom;
    
    nucLatticePosTemp = nucLatticePos;
    nucLatticeMomTemp = nucLatticeMom;
    
    %Inserting Impurities
    if ~isempty(impurityIdx)
        magLatticeMomTemp(impurityIdx,:) = zeros( numImp, 3);
         
        %Copper Impurity substituted into iron positions
        magLatticePosCu = magLatticePosTemp(impurityIdx,:);
        magLatticeNucCu = nucMagnetonCu * nucMagneton * normr(rand(numImp,3) - 0.5);
    else %no impuriites
        magLatticePosCu = [];
        magLatticeNucCu = [];
    end
    
    %Assembling complete lattice (with copper impurities)
    fullLatticePos = [magLatticePosTemp; nucLatticePosTemp; magLatticePosCu];
    fullLatticeMom = [magLatticeMomTemp; nucLatticeMomTemp; magLatticeNucCu];

    %Compute normalized deviations
    fullDeviationPos = fullLatticePos - repmat(muonSite,size(fullLatticePos,1),1);
    fullDeviationNorm = sqrt( nansum(fullDeviationPos.^2, 2) );
    fullDeviationNormPos = fullDeviationPos ./ repmat(fullDeviationNorm,1,3);
    
    %Compute Dipolar Field at muon position
    magField = ( 3 * repmat( dot( fullLatticeMom, fullDeviationNormPos, 2), 1, 3) ...
        .* fullDeviationNormPos - fullLatticeMom ) ./ repmat( fullDeviationNorm, 1, 3).^3;

    %Compute Magnetic Field at muon site
    totalField(j,:) = nansum(magField);
end

