function ZeemanShift = NVZeemanShiftFromMagnetizedSampleEdge(...
    x, x0, Is, heightNV, phiNV, thetaNV)
% Calculates the Zeeman shift of an NV center (with NV axis defined by thetaNV, phiNV)
% that is scanned at a height z=heightNV along a line crossing the edge
% of a homogenously, out-of-plane magnetized 2D sample at z=0.

% Input Parameters
%   x (in um): position array of the NV center
%   x0 (in um): position of the sample edge
%   Is (in uA): 2D magnetisation of the sample
%   heightNV (in um): distance of the NV center above the sample
%   phiNV (in degree): azimuthal angle from x-axis to NV axis
%   thetaNV (in degree): polar angle from z-axis to NV axis

% Output Parameters
%   ZeemanShift (in MHz): Zeeman shift experienced by the NV center

    % Vacuum permeability
    mu0 = 4*pi*10^(-7); % Tum/uA
    % NV gyromagnetic ratio in MHz/T
    gammaNV = 28000;

    % Field of a Bloch wall
    Bx =  mu0 * Is/pi * heightNV./((x-x0).^2 + heightNV.^2);
    By = zeros(1,length(x));
    Bz = -mu0 * Is/pi * (x-x0)./((x-x0).^2 + heightNV.^2);
    BBloch = [Bx; By; Bz];

    % sample edge is half a Bloch wall
    Bedge = BBloch/2;

    % Field projection onto NV center
    NVx = sind(thetaNV) * cosd(phiNV);
    NVy = sind(thetaNV) * sind(phiNV);
    NVz = cosd(thetaNV);
    BNV = NVx*Bedge(1,:) + NVy*Bedge(2,:) + NVz*Bedge(3,:);

    % ZeemanShift in MHz
    ZeemanShift = gammaNV*BNV;
end



