function GyroidFieldProfiles(Field,Idx,x,y,z,xcut,ycut,zcut)
%Use to create a 3D model of the gyroid metamaterial with electric field
%intensity enhancement cross sections at designated planes
%--------------------------------------------------------------------------
%Inputs
%Field - Electric field intensity enhancement. Should be a real-valued 3D 
% matrix with dimensions [length(x),length(y),length(z)].
% Use: Field = log10((E2_LCP + E2_RCP)./2) for log scale, unpolarized
% intensity enhancement.
%Idx - Spatially resolved complex refractive index with dimensions
% [length(x),length(y),length(z)].
%x, y, z - 1D spatial coordinate vectors.
%xcut, ycut, zcut - 1D vectors denoting the intensity enhancement cross
% section planes.
%--------------------------------------------------------------------------
%Example: GyroidFieldProfiles(log10((E2_LCP(:,:,:,2)+E2_RCP(:,:,:,2))./2),idx,x,y,z,max(y)*2,max(z)*2,[-1 1].*(65.*2.*sqrt(2)-2));
%--------------------------------------------------------------------------
%Reorginze the intensity and refractive index matrices to swap the x and z axes
Idx = permute(Idx,[3 2 1]);
Field = permute(Field,[3 2 1]);

%Expand the matrices to plot 2 unit cells by 2 unit cells
Field = [Field,Field;Field,Field];
Idx = [Idx,Idx;Idx,Idx];

%Update the y and z coordinates to account for the expansion
y = linspace(min(y)*2,max(y)*2,length(y)*2);
z = linspace(min(z)*2,max(z)*2,length(z)*2);

%Create a mesh grid with the x and z axis swapped
[X,Y,Z] = meshgrid(y,z,-x);

%Consider only the real part of the refractive index
Idx = real(Idx);

%Find the minimum and maximum values of the refractive index
Idx_max = max(max(max(Idx)));
Idx_min = min(min(min(Idx)));

if Idx_max == 1 %If the gyroid's refractive index is less than air
    %Calculate the faces and vectors needed to create the gyroid 3D model
    [f,v] = isosurface(X,Y,Z,Idx,Idx_min+0.001);
    [f2,v2] = isocaps(X,Y,Z,Idx,Idx_min+0.001,'below');
    f3 = [f ; f2+length(v(:,1))];
    v3 = [v ; v2];
else %Else the gyroid's refractive index is greater than air
    [f,v] = isosurface(X,Y,Z,Idx,Idx_max-0.001);
    [f2,v2] = isocaps(X,Y,Z,Idx,Idx_max-0.001,'above');
    f3 = [f ; f2+length(v(:,1))];
    v3 = [v ; v2];
end
%Plot the gyroid 3D model
figure;
p = patch('Faces',f3,'Vertices',v3);
p.FaceColor = [0.8 0.8 0.8];
p.EdgeColor = 'none';
daspect([1 1 1])
xlabel('x')
ylabel('y')
zlabel('z')
view(3)
axis off
camlight
lighting gouraud

hold on

%Add the intensity enhancement cross sections
h = slice(X,Y,Z,real(Field),xcut,ycut,zcut);
set(h,'edgecolor','none');
set(h,'FaceLighting','none');
colorbar;
end