%% Variables(change value if needed)
%number of pixels for each dimension
    N = 256;
%refractive index
    refractiveIndex = 1.333;
%wavelength[nm]
    wavelength_SI = 488e-9;%excitation 
    wavelength_sheet = 405e-9;%single photon activation 
    wavelength_sheet2p = 780e-9;%two photonactivation
    wavelength_det = 522e-9;% detection (Skylan-NS) 
%parameters for structured illumination
    ar = 0.6;%aperture ratio
    Nangle = 3;%number of strip angles
    Nphase = 5;%number of strip phases
%NA
    NA_SI_det =1.1;%detection objective
    NA_sheet =0.67;%effective NA of activation, high NA, standard lightsheet
    lowNA_sheet = 0.1;%effective NA of activation, low NA, standard lightsheet
    NA_besselsheet =0.67;%effective NA of activation, sheet by scanned Bessel beam

%Image size [INPUT image size as a number of waves in excitation wavelength]
    Nlambda = 14;

%% import pupil function and samples
%Pupils(choose one if needed.)
    load('pupils/pupils.mat')
    
    pupil_SI =SIM_AR06;% SIM, AR = 0.6. Choose one if needed.
    % pupil_SI =SIM_AR09;% SIM, AR = 0.9. Choose one if needed.

    pupil_sheet = sheet;%standard lightsheet
    pupil_bessel = bessel;%bessel beam
    pupil_psf = Ppsf;%detection 

%samples(choose one if needed.)
    load('sample/point.mat');%point object
    % load('sample/scatterpoints_d200nm_n5000.mat');%dispersed 5000 beads with a diameter of 200 nm 
    % load('sample/sphere_d10um.mat');%single bead with a diameter of 10 um
    
    xyz_sample = obj;
    %zero pudding for sample scanning
    xyz_sample_scan = cat(3,zeros(size(xyz_sample,1),size(xyz_sample,2),size(xyz_sample,3)/2-1),xyz_sample,zeros(size(xyz_sample,1),size(xyz_sample,2),size(xyz_sample,3)/2));

%% Calibration of image size for the other wavelength
    LL_SI =Nlambda*refractiveIndex;
    LL_sheet =Nlambda*refractiveIndex*wavelength_SI/wavelength_sheet;
    LL_sheet2p =Nlambda*refractiveIndex*wavelength_SI/wavelength_sheet2p;
    LL_psf =Nlambda*refractiveIndex*wavelength_SI/wavelength_det;

%% 3D distribution of sheet activation
xyz_sheet =permute(FP_LightDistribution_vol(N,NA_sheet/refractiveIndex,LL_sheet,pupil_sheet, "intensity"), [2 3 1]);%high NA, standard lightsheet
xyz_sheet_lowNA = permute(FP_LightDistribution_vol(N,lowNA_sheet/refractiveIndex,LL_sheet,pupil_sheet, "intensity"), [2 3 1]);%low NA, standard lightsheet
xyz_besselsheet_temp =FP_LightDistribution_vol(N,NA_besselsheet/refractiveIndex,LL_sheet,pupil_bessel, "intensity");%Bessel beam, single photon activation
xyz_besselsheet2p_temp =FP_LightDistribution_vol(N,NA_besselsheet/refractiveIndex,LL_sheet2p,pupil_bessel, "intensity").^2;%Bessel beam, two photon activation

%% detection PSF
xyz_psf =FP_LightDistribution_vol(N,NA_SI_det/refractiveIndex,LL_psf,pupil_psf, "intensity");

%% patterns of structured illumination
xyz_SI = simamaker3d(N,Nlambda*wavelength_SI,wavelength_SI,refractiveIndex, NA_SI_det, ar, Nangle, Nphase, 1);
%change strip direction
xyz_SI = permute(xyz_SI,[2 1 3 4]);

%% distribution of sheet by scanned Bessel beam
xyz_besselbeam = permute(xyz_besselsheet_temp, [2 3 1]);
base = zeros (N,2*N-1,N);
for i=1:N
    base(:,i:i+N-1,:) =base(:,i:i+N-1,:) + xyz_besselsheet_temp;
end
xyz_besselsheet = permute(base(:,N/2:N/2+N-1,:), [2 3 1]);%single photon activation

xyz_besselbeam2p = permute(xyz_besselsheet2p_temp, [2 3 1]);
base2p = zeros (N,2*N-1,N);
for i=1:N
    base2p(:,i:i+N-1,:) =base2p(:,i:i+N-1,:) + xyz_besselsheet2p_temp;
end
xyz_besselsheet2p = permute(base2p(:,N/2:N/2+N-1,:), [2 3 1]);%two photon activation

%% convolution for image.
xy_image_SI = zeros(N,N,N,Nangle*Nphase);
xy_image_sheetSI = zeros(N,N,N,Nangle*Nphase);
xy_image_lowNAsheetSI = zeros(N,N,N,Nangle*Nphase);
xy_image_besselsheetSI = zeros(N,N,N,Nangle*Nphase);
xy_image_besselsheet2pSI = zeros(N,N,N,Nangle*Nphase);
tic

for i= 1:Nangle*Nphase
    xyz_sheetSI = xyz_SI(:,:,:,i).*xyz_sheet;
    xyz_lowNAsheetSI = xyz_SI(:,:,:,i).*xyz_sheet_lowNA;
    xyz_besselsheetSI = xyz_SI(:,:,:,i).*xyz_besselsheet;
    xyz_besselsheet2pSI = xyz_SI(:,:,:,i).*xyz_besselsheet2p;

    for k = 1:size(xyz_sample,3)
        disp(['(pattern,z) = (',num2str(i),',',num2str(k),')']);
        toc
        %scan object
        xyz_sample_temp = xyz_sample_scan(:,:,k:k+size(xyz_sample,3)-1);
        
        %conventional
        image_SI = conv3_fft(xyz_SI(:,:,:,i).*xyz_sample_temp,xyz_psf);
        image_SI = image_SI(N/2:3*N/2-1,N/2:3*N/2-1,N/2:3*N/2-1);
        xy_image_SI(:,:,k,i) = image_SI(:,:,N/2);%image with conventional structured illumination
        
        %highNA sheet
        image_sheetSI= conv3_fft(xyz_sheetSI.*xyz_sample_temp,xyz_psf);
        image_sheetSI = image_sheetSI(N/2:3*N/2-1,N/2:3*N/2-1,N/2:3*N/2-1);
        xy_image_sheetSI(:,:,k,i) = image_sheetSI(:,:,N/2);%image with high NA standard lightsheet & structured illumination 
        
        %lowNAsheet
        image_lowNAsheetSI= conv3_fft(xyz_lowNAsheetSI.*xyz_sample_temp,xyz_psf);
        image_lowNAsheetSI = image_lowNAsheetSI(N/2:3*N/2-1,N/2:3*N/2-1,N/2:3*N/2-1);
        xy_image_lowNAsheetSI(:,:,k,i) = image_lowNAsheetSI(:,:,N/2);%image with low NA standard lightsheet & structured illumination 

        %bessel1photon
        image_besselsheetSI= conv3_fft(xyz_besselsheetSI.*xyz_sample_temp,xyz_psf);
        image_besselsheetSI = image_besselsheetSI(N/2:3*N/2-1,N/2:3*N/2-1,N/2:3*N/2-1);
        xy_image_besselsheetSI(:,:,k,i) = image_besselsheetSI(:,:,N/2);%image with single photon, scanned Bessel sheet & structured illumination 

        %bessel2photon
        image_besselsheet2pSI= conv3_fft(xyz_besselsheet2pSI.*xyz_sample_temp,xyz_psf);
        image_besselsheet2pSI = image_besselsheet2pSI(N/2:3*N/2-1,N/2:3*N/2-1,N/2:3*N/2-1);
        xy_image_besselsheet2pSI(:,:,k,i) = image_besselsheet2pSI(:,:,N/2);%image with two photon, scanned Bessel sheet & structured illumination
    end
end

toc