function [magLattice,nucLattice] = Make_NaFeAs_Lattice(numCell, spinFe)
%function [magLattice,nucLattice] = Make_NaFeAs_Lattice(numCell, spinFe)
%
%Computes nuclear and magnetic lattice of NaFeAs for a supercell
%
%
%INPUTS
%
%numCell                        Crystal length = 2*numCell + 1 in 
%                                   each dimension
%
%spinFe                         Spin vector for each iron in unit cell
%                                   (2 x 3 matrix)
%
%
%OUTPUTS
%
%magLattice                     magnetic lattice
%                                   Cell 1: positions
%                                   cell 2: spins
%
%nucLattice                     nuclear lattice
%                                   Cell 1: positions
%                                   Cell 2: spins
%
%XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
%X Author:          Sky Cheung                                                X
%X Date Created:    20 Apr 2016                                               X
%X Last Modified:   15 Jun 2016                                               X
%XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
%
%Input Parsing
p = inputParser;

%Parsing conditions
addRequired(p, 'numCell', @(x) validateattributes( x, ...
    {'numeric'}, {'nonempty','scalar','positive'} ) );
addRequired(p, 'spinFe', @(x) validateattributes( x, ...
    {'numeric'}, {'nonempty','ncols',3,'nrows',2} ) );

%Parsing Input
parse(p, numCell, spinFe);

%Retrieving Proper Inputs
numCell = p.Results.numCell;
spinFe = p.Results.spinFe;

%Deleting Parser
clear p;

%MATERIAL PARAMETERS

%Import Unit Cell Information (Angstrom to cm)
a = 3.94940 * 10^-8;
b = 3.94940 * 10^-8;
c = 7.03960 * 10^-8;

%Electronic moments
bohrMagnetonFe = 1.1;
bohrMagneton = 9.274 * 10 ^ (-21); %erg/G

%Nuclear moments
nucMagnetonFe = 0.09024;
nucMagnetonNa = 2.2174;
nucMagnetonAs = 1.439;
%nucMagnetonCu = 2.2226;
nucMagneton = 5.051 * 10^(-24);

%CONSTRUCTION
%Unit cell (ignore Fe layer on top, so 8 - 2 = 6 atoms in unit cell)
posFe1 = [3/4 * a, 1/4 * b, 0];
posFe2 = [1/4 * a, 3/4 * b, 0];

posNa1 = [3/4 * a, 3/4 * b, 0.35398 * c];
posNa2 = [1/4 * a, 1/4 * b, 0.64602 * c];

posAs1 = [3/4 * a, 3/4 * b, 0.79722 * c];
posAs2 = [1/4 * a, 1/4 * b, 0.20278 * c];

%FESC Stripe pattern

momentFe1 = spinFe(1,:);
momentFe2 = spinFe(2,:);

%Make lattice
lengthCrystal = 2 * numCell + 1;

latticeVec = cell(1,3);
latticeFe1 = cell(1,3);
latticeFe2 = cell(1,3);
latticeNa1 = cell(1,3);
latticeNa2 = cell(1,3);
latticeAs1 = cell(1,3);
latticeAs2 = cell(1,3);

%Lattice Vectors
[latticeVec{1},latticeVec{2},latticeVec{3}] = meshgrid( ...
    a * (-numCell:1:numCell) ,...
    b * (-numCell:1:numCell) ,...
    c * (-numCell:1:numCell) );

%Positions of all atoms in supercell
for j = 1:3
    latticeFe1{j} = posFe1(j) + latticeVec{j};
    latticeFe2{j} = posFe2(j) + latticeVec{j};
    latticeNa1{j} = posNa1(j) + latticeVec{j};
    latticeNa2{j} = posNa2(j) + latticeVec{j};
    latticeAs1{j} = posAs1(j) + latticeVec{j};
    latticeAs2{j} = posAs2(j) + latticeVec{j};
end

%Make AFM Tensor
M = (-1).^(1:lengthCrystal);
AFMPlane = M'*M;
AFMTensor = repmat(AFMPlane,[1,1,lengthCrystal]);
AFMTensor(:,:,2:2:end) = -1 * AFMTensor(:,:,2:2:end);

momFe1 = cell(1,3);
momFe2 = cell(1,3);

%Local Electronic Moments
for j = 1:3
    momFe1{j} = bohrMagnetonFe * bohrMagneton * momentFe1(j) * AFMTensor;
    momFe2{j} = bohrMagnetonFe * bohrMagneton * momentFe2(j) * AFMTensor;
end

numAtom = numel(momFe1{1});

%Nuclear Dipolar moments
nucFe1 = nucMagnetonFe * nucMagneton * normr(rand(numAtom,3) - 0.5);
nucFe2 = nucMagnetonFe * nucMagneton * normr(rand(numAtom,3) - 0.5);
nucNa1 = nucMagnetonNa * nucMagneton * normr(rand(numAtom,3) - 0.5);
nucNa2 = nucMagnetonNa * nucMagneton * normr(rand(numAtom,3) - 0.5);
nucAs1 = nucMagnetonAs * nucMagneton * normr(rand(numAtom,3) - 0.5);
nucAs2 = nucMagnetonAs * nucMagneton * normr(rand(numAtom,3) - 0.5);

%condensing moments and positions of all field-generating atoms
magLatticePosFe1 = [latticeFe1{1}(:),latticeFe1{2}(:),latticeFe1{3}(:)];
magLatticeMomFe1 = [momFe1{1}(:),momFe1{2}(:),momFe1{3}(:)];
magLatticePosFe2 = [latticeFe2{1}(:),latticeFe2{2}(:),latticeFe2{3}(:)];
magLatticeMomFe2 = [momFe2{1}(:),momFe2{2}(:),momFe2{3}(:)];

magLatticeNucFe1 = nucFe1;
magLatticeNucFe2 = nucFe2;
magLatticePosNa1 = [ latticeNa1{1}(:), latticeNa1{2}(:), latticeNa1{3}(:) ];
magLatticeNucNa1 = nucNa1;
magLatticePosNa2 = [ latticeNa2{1}(:), latticeNa2{2}(:), latticeNa2{3}(:) ];
magLatticeNucNa2 = nucNa2;
magLatticePosAs1 = [ latticeAs1{1}(:), latticeAs1{2}(:), latticeAs1{3}(:) ];
magLatticeNucAs1 = nucAs1;
magLatticePosAs2 = [ latticeAs2{1}(:), latticeAs2{2}(:), latticeAs2{3}(:) ];
magLatticeNucAs2 = nucAs2;

%Undoped Lattice
magLatticePos = [magLatticePosFe1;magLatticePosFe2];
magLatticeMom = [magLatticeMomFe1;magLatticeMomFe2];

nucLatticePos = [...
    magLatticePosFe1;magLatticePosFe2;...
    magLatticePosNa1;magLatticePosNa2;...
    magLatticePosAs1;magLatticePosAs2;...
    ];
nucLatticeMom = [...
    magLatticeNucFe1;magLatticeNucFe2;...
    magLatticeNucNa1;magLatticeNucNa2;...
    magLatticeNucAs1;magLatticeNucAs2;...
    ];

%Storing as cell
magLattice =  {magLatticePos, magLatticeMom};
nucLattice =  {nucLatticePos, nucLatticeMom};
