function [N_sk,Bem] = f_Bem(nn,xsection,rate)
%   f_Bem - Calculating emergent magnetic field (Bem) and skyrmion number (Nsk) for a given
%   xsection of a 3D vector field nn
%
%   Input arguments:
%   nn - An mxnxpx3 array of 3D unit vector. m, n, p are the number of
%   points in x, y, and z, respectively. The last dimension specifies the
%   vecotr oientation.
%   xsection - 'xz', 'yz', or 'xy'. Specifies the cross-section in which
%   the emergent field and the skyrmion number density are calculated
%   rate - a positive integer specifying the rate of interpolation to
%   improve numerical accuracy
%
%   Output:
%   Nsk - Skymrion number calculated in the xsection of choice
%   Bem - Emergent magentic field in the xsection of choice. In units of
%   hbar/2e
%
%   Author: Jung-Shen Benny Tai. 11/08/2023
%
%   Example:
%   [N_sk,Bem] = f_Bem(nn,xsection,rate)

% Interpolate field
[m,n,p,~] = size(nn);

xxx = 0:m-1;
yyy = 0:n-1;
zzz = 0:p-1;

xxxq = 0:1/rate:m-1;
yyyq = 0:1/rate:n-1;
zzzq = 0:1/rate:p-1;

dxq = xxxq(2)-xxxq(1);
dyq = yyyq(2)-yyyq(1);
dzq = zzzq(2)-zzzq(1);

[yy,xx,zz] = meshgrid(yyy,xxx,zzz);
[yyq,xxq,zzq] = meshgrid(yyyq,xxxq,zzzq);

nnq(:,:,:,1) = interp3(yy,xx,zz,nn(:,:,:,1),yyq,xxq,zzq,'cubic');
nnq(:,:,:,2) = interp3(yy,xx,zz,nn(:,:,:,2),yyq,xxq,zzq,'cubic');
nnq(:,:,:,3) = interp3(yy,xx,zz,nn(:,:,:,3),yyq,xxq,zzq,'cubic');

% % normalize director field after interpolation
modnn = sqrt(nnq(:,:,:,1).^2+nnq(:,:,:,2).^2+nnq(:,:,:,3).^2);
nnq(:,:,:,1) = nnq(:,:,:,1)./modnn;
nnq(:,:,:,2) = nnq(:,:,:,2)./modnn;
nnq(:,:,:,3) = nnq(:,:,:,3)./modnn;

% Creating Skyrmion number density or emergent B field proportional to E(ijk)nn.(d_j nn cross d_k nn)
partial_nn = zeros(size(xxq,1),size(xxq,2),size(xxq,3),3,3);
% index: x,y,z coordinates, partial_i, n_j

[partial_nn(:,:,:,2,1),partial_nn(:,:,:,1,1),partial_nn(:,:,:,3,1)] = gradient(nnq(:,:,:,1),dxq,dyq,dzq);
[partial_nn(:,:,:,2,2),partial_nn(:,:,:,1,2),partial_nn(:,:,:,3,2)] = gradient(nnq(:,:,:,2),dxq,dyq,dzq);
[partial_nn(:,:,:,2,3),partial_nn(:,:,:,1,3),partial_nn(:,:,:,3,3)] = gradient(nnq(:,:,:,3),dxq,dyq,dzq); 

% emergent B field in units of hbar
% Bem = zeros(rate*m-1,rate*n-1,rate*p-1,3);
Bem(:,:,:,1) = dot(nnq(:,:,:,:),cross(squeeze(partial_nn(:,:,:,2,:)),squeeze(partial_nn(:,:,:,3,:)),4),4);
Bem(:,:,:,2) = dot(nnq(:,:,:,:),cross(squeeze(partial_nn(:,:,:,3,:)),squeeze(partial_nn(:,:,:,1,:)),4),4);
Bem(:,:,:,3) = dot(nnq(:,:,:,:),cross(squeeze(partial_nn(:,:,:,1,:)),squeeze(partial_nn(:,:,:,2,:)),4),4);

N_sk_density = Bem/(4*pi);

% Calculate skyrmion number
switch xsection
    case 'xy'
        N_sk_density2D = N_sk_density(:,:,round((1+size(N_sk_density,3)/2)),3);
        N_sk = trapz(trapz(N_sk_density2D))*dxq*dyq;
    case 'xz'
        N_sk_density2D = N_sk_density(:,round((1+size(N_sk_density,2)/2)),:,2);
        N_sk = trapz(trapz(N_sk_density2D))*dxq*dzq;
    case 'yz'
        N_sk_density2D = N_sk_density(round((1+size(N_sk_density,1)/2)),:,:,1);
        N_sk = trapz(trapz(N_sk_density2D))*dyq*dzq;
end
disp(['Skyrmion number in ' xsection '-midplane = ' num2str(N_sk)])


% Plot
imagesc(squeeze(N_sk_density2D)')
colorbar
axis equal off

end