function Mz = MzFromBNV(...
    dpx, dpy, BNV, heightNV, thetaNV, phiNV, lambdaCrit, FFTPaddingFactor)
% Retrieves the 2D out-of-plane magnetization Mz(z=0) from a planar map of BNV 
% with NV axis (thetaNV, phiNV) at a height 'heightNV' above a given
% magnetization.

% Input Parameters
%   dpx (in m): real space pixel spacing in x direction
%   dpy (in m): real space pixel spacing in y direction
%   BNV (in Tesla): magnetic field matrix along NV axis
%   heightNV (in m): BNV(z=heightNV) --> Mz(z=0)
%   thetaNV (in degree): polar angle from z-axis to NV axis
%   phiNV (in degree): azimuthal angle from x-axis to NV axis
%   lamdbaCrit (in m): cut-off wavelength for filter
%   FFTPaddingFactor (scalar): padding factor of FFT

% Output Parameters
%   Mz (in Ampere = muB/nm^2) is 2D magnetization at z=0

    u0 = 4*pi*1e-7; % Tm/A permeability of free space

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % SETUP K-VECTOR
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    [nbY,nbX] = size(BNV); 
    nbKx = nbX*FFTPaddingFactor;
    nbKy = nbY*FFTPaddingFactor;
    [kx, ky, k] = kvalues(nbX, nbY, dpx, dpy, FFTPaddingFactor);

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % TRANSFORMATION MATRICES
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    expTerm = exp(-k*heightNV)/2;
    Tx = -1i*u0* kx .* expTerm;
    Ty = -1i*u0* ky .* expTerm;
    Tz =     u0* k  .* expTerm;

    % Projection onto NV axis
    TNV = sind(thetaNV)*cosd(phiNV) * Tx +...
          sind(thetaNV)*sind(phiNV) * Ty +...
          cosd(thetaNV)             * Tz;

    % Filter only used for back propagation (BNV --> magnetization)
    kkrit = 2*pi/lambdaCrit;
    W = 0.5*(1 + cos(k*lambdaCrit/2));
    W(abs(k)>kkrit) = 0;

    % Final Transformation Matrix
    Trans_BtoM = W./TNV;

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % CALCULATE MAGNETIZATION
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    bnv = fftshift(fft2(BNV, nbKy, nbKx));
    % get rid of DC component
    bnv(k<eps) = 0;

    mz = bnv.*Trans_BtoM;
    Mext = real(ifft2( ifftshift (mz)));
    Mz = Mext(1:nbY , 1:nbX); 
end
