function f_3PEF_PM(nn,xsection,pol,ratexs)
%   f_3PEF_PM - Calculating and visulaizing 3PEF-PM imaging of a given mid-plane xsection 
%   of a 3D director 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 - 1, 2, or 3, corrsponding to xy, xz, and yz cross-sections, 
%   respectively. This input pecifies the cross-section in which the 3PEF-PM 
%   iamge is calculated.
%   pol - 1, 2, or 3, corresponding to linear polarization along the 1st
%   dimension (x), linear polarization along the 2nd dimension (2), and
%   circular polarization in xy, respectively.
%   ratexs - a positive integer specifying the rate of interpolation.
%
%   Output:
%   A simulated 3PEF-PM image
%
%   Author: Jung-Shen Benny Tai. 11/08/2023
%
%   Example:
%   f_3PEF_PM(nn,xsection,pol,ratexs)

% Obtain midplane cross-section 
[m,n,p,~] = size(nn);

switch xsection
case 1 % xy-midplane
if mod(p,2) == 1
    nnxymid = nn(:,:,round((p+1)/2),:);
else 
    nnxymid = (nn(:,:,round((p+1)/2)-1,:)+nn(:,:,round((p+1)/2),:))/2;
    modmid = sqrt(nnxymid(:,:,:,1).^2+nnxymid(:,:,:,2).^2+nnxymid(:,:,:,3).^2);
    nnxymid(:,:,:,1) = nnxymid(:,:,:,1)./modmid;
    nnxymid(:,:,:,2) = nnxymid(:,:,:,2)./modmid;
    nnxymid(:,:,:,3) = nnxymid(:,:,:,3)./modmid;
end
        
case 2 % xz-midplane
if mod(n,2) == 1
    nnxzmid = nn(:,round((n+1)/2),:,:);
else 
    nnxzmid = (nn(:,round((n+1)/2)-1,:,:)+nn(:,round((n+1)/2),:,:))/2;
    modmid = sqrt(nnxzmid(:,:,:,1).^2+nnxzmid(:,:,:,2).^2+nnxzmid(:,:,:,3).^2);
    nnxzmid(:,:,:,1) = nnxzmid(:,:,:,1)./modmid;
    nnxzmid(:,:,:,2) = nnxzmid(:,:,:,2)./modmid;
    nnxzmid(:,:,:,3) = nnxzmid(:,:,:,3)./modmid;
end

case 3 % yz-midplane
if mod(m,2) == 1
    nnyzmid = nn(round((m+1)/2),:,:,:);
else 
    nnyzmid = (nn(round((m+1)/2)-1,:,:,:)+nn(round((m+1)/2),:,:,:))/2;
    modmid = sqrt(nnyzmid(:,:,:,1).^2+nnyzmid(:,:,:,2).^2+nnyzmid(:,:,:,3).^2);
    nnyzmid(:,:,:,1) = nnyzmid(:,:,:,1)./modmid;
    nnyzmid(:,:,:,2) = nnyzmid(:,:,:,2)./modmid;
    nnyzmid(:,:,:,3) = nnyzmid(:,:,:,3)./modmid;
end

end

% Interpolation 
switch xsection
    case 1
[mxs,nxs,~,~] = size(nnxymid);
x = linspace(1,mxs,mxs*ratexs);
y = linspace(1,nxs,nxs*ratexs);

[yy,xx] = meshgrid(y,x);
nn_int2(:,:,1) = interp2(squeeze(nnxymid(:,:,:,1)),yy,xx,'cubic');
nn_int2(:,:,2) = interp2(squeeze(nnxymid(:,:,:,2)),yy,xx,'cubic');
nn_int2(:,:,3) = interp2(squeeze(nnxymid(:,:,:,3)),yy,xx,'cubic');

modnn = sqrt(nn_int2(:,:,1).^2+nn_int2(:,:,2).^2+nn_int2(:,:,3).^2);
nn_int2(:,:,1) = nn_int2(:,:,1)./modnn;
nn_int2(:,:,2) = nn_int2(:,:,2)./modnn;
nn_int2(:,:,3) = nn_int2(:,:,3)./modnn;

nn_slice(:,:,1,:) = nn_int2;

figure
if pol <= 2
    imagesc(squeeze(nn_slice(:,:,1,pol))'.^6)
else
    imagesc((squeeze(nn_slice(:,:,1,1))'.^2+squeeze(nn_slice(:,:,1,2))'.^2).^3)
end
    case 2
[mxs,~,pxs,~] = size(nnxzmid);
x = linspace(1,mxs,mxs*ratexs);
z = linspace(1,pxs,pxs*ratexs);

[zz,xx] = meshgrid(z,x);
nn_int2(:,:,1) = interp2(squeeze(nnxzmid(:,:,:,1)),zz,xx,'cubic');
nn_int2(:,:,2) = interp2(squeeze(nnxzmid(:,:,:,2)),zz,xx,'cubic');
nn_int2(:,:,3) = interp2(squeeze(nnxzmid(:,:,:,3)),zz,xx,'cubic');

modnn = sqrt(nn_int2(:,:,1).^2+nn_int2(:,:,2).^2+nn_int2(:,:,3).^2);
nn_int2(:,:,1) = nn_int2(:,:,1)./modnn;
nn_int2(:,:,2) = nn_int2(:,:,2)./modnn;
nn_int2(:,:,3) = nn_int2(:,:,3)./modnn;

nn_slice(:,1,:,:) = nn_int2;

figure
if pol <= 2
    imagesc(squeeze(nn_slice(:,1,:,pol))'.^6)
else
    imagesc((squeeze(nn_slice(:,1,:,1))'.^2+squeeze(nn_slice(:,1,:,2))'.^2).^3)
end

case 3
[~,nxs,pxs,~] = size(nnyzmid);
y = linspace(1,nxs,nxs*ratexs);
z = linspace(1,pxs,pxs*ratexs);

[zz,yy] = meshgrid(z,y);
nn_int2(:,:,1) = interp2(squeeze(nnyzmid(:,:,:,1)),zz,yy,'cubic');
nn_int2(:,:,2) = interp2(squeeze(nnyzmid(:,:,:,2)),zz,yy,'cubic');
nn_int2(:,:,3) = interp2(squeeze(nnyzmid(:,:,:,3)),zz,yy,'cubic');

modnn = sqrt(nn_int2(:,:,1).^2+nn_int2(:,:,2).^2+nn_int2(:,:,3).^2);
nn_int2(:,:,1) = nn_int2(:,:,1)./modnn;
nn_int2(:,:,2) = nn_int2(:,:,2)./modnn;
nn_int2(:,:,3) = nn_int2(:,:,3)./modnn;

nn_slice(1,:,:,:) = nn_int2;

% plot 
figure
if pol <= 2
    imagesc(squeeze(nn_slice(1,:,:,pol))'.^6)
else
    imagesc((squeeze(nn_slice(1,:,:,1))'.^2+squeeze(nn_slice(1,:,:,2))'.^2).^3)
end

end

set(gca,'YDir','normal')
axis equal off

% colormap
k=15;
cmap = zeros(2^k,3);
t = 0:2^k-1;
cmap(:,2)=t/(2^k-1);
colormap(cmap)
set(gcf,'Color',[1 1 1])


end