function [E2_LCP_int,E2_RCP_int] = Integrate_E2(E2_LCP,E2_RCP,y,z)
%Calculate the integrated electric field intensity enhancement down the
%length of a gyroid metamaterial
%--------------------------------------------------------------------------
%Inputs
%E2_LCP, E2_RCP - Electric field intensity enhancements under LCP or RCP
% illumination, respectively. E2_LCP and E2_RCP should be a real-valued 4D 
% matrix with dimensions [length(x),length(y),length(z),length(f)]
%y, z - 1D spatial coordinate vectors
%--------------------------------------------------------------------------
%Outputs
%E2_LCP_int, E2_RCP_int - integrated electric field intensity enhancement 
% down the length of a gyroid metamaterial. Size = 
% [length(x),length(Wvlgth_EM]
%--------------------------------------------------------------------------

E2_LCP_int = squeeze(trapz(z,trapz(y,permute(E2_LCP,[2,3,1,4])))); %Integrate the LCP electric field intensity enhancement along the y and z directions
E2_LCP_int = E2_LCP_int./squeeze(trapz(z,trapz(y,permute(ones(size(E2_LCP)),[2,3,1,4])))); %Normalize to the integrated area

%Repeat for RCP illumination
E2_RCP_int = squeeze(trapz(z,trapz(y,permute(E2_RCP,[2,3,1,4]))));
E2_RCP_int = E2_RCP_int./squeeze(trapz(z,trapz(y,permute(ones(size(E2_RCP)),[2,3,1,4]))));

end