function SI=simamaker3d(N, imsize, wavelength, refractiveindex, NA, ar, Nangle, Nphase, zeroOrder)
%Copyright Ryosuke Oketani 2022, Ver1.1
% make structured illumination pattern in 3D with Planewave_pol.m function%
%
%----------
%Syntax
%    simamaker3d(N,imsize,wavelength,refractiveindex, NA, ar, Nangle, Nphase)
%    simamaker3d(N,imsize,wavelength,refractiveindex, NA, ar, Nangle, Nphase,zeroOrder)
%----------
%
%Variables
%   N: number of pixel.
%   imsize: image sizeof calculated area.
%   wavelength: wavelength of incident light.
%   refractiveindex: refractive index of surrounding media.
%   NA: NA of illumination objective lens. 
%   ar: aperture ratio of beam position at pupil. 
%   Nangle: number of angles
%   Nphase: number of phases
%   zeroOrder:0 or 1. When 1, 0th order planewave used. default: 0.
%
%Example: simamaker3d(128, 3e-6, 488e-9, 1.333, 1.2, 0.6, 3, 5, 1)
%% 
if nargin<9
    zeroOrder = 0;
end

%% parameters
Ns = [N N N];
imsizes = [imsize imsize imsize];
Kvalue = 2*pi*refractiveindex/wavelength;
theta1 = asin(NA*ar/refractiveindex);
theta2 = -asin(NA*ar/refractiveindex);
theta0th = 0;
phase1 = 0;
phase0th = 0;

%% make SI patterns
SI = zeros(Ns(1),Ns(2),Ns(3),Nangle*Nphase);
for i = 1:Nangle
    phi = (i-1)*2*pi/Nangle;
    for k = 1:Nphase
        phase2 = (k-1)*2*pi/Nphase;
        SI(:,:,:,k+(i-1)*Nphase) = abs((Planewave_pol(Ns,imsizes,[theta1 phi],Kvalue,[],phase1-phase2.*zeroOrder)+Planewave_pol(Ns,imsizes,[theta2 phi],Kvalue,[],phase2) + zeroOrder .* Planewave_pol(Ns, imsizes,[theta0th phi],Kvalue,[],phase0th)).^2);
    end
end

end