%init
clc; clear all 
close all

% z-axis control parameters
z2_under_image = -300e-6;  %[um]
z2_upper_image = 300e-6;  %[um]
z2_series_step = 1e-6;   %[um]
z2_series_number = (z2_upper_image - z2_under_image)/z2_series_step + 1;

%constants
% integrating variabes in h : dx, dy ==> h = f(u,v)
% coordinate indices in image planes : u,v
syms x y u v;

syms M; %z2 = (-1)* z_under_image : z_series_number : z_upper_image;
wavelength_0 = 365e-9; %[m]
n = 1.47; % refractive index
wavelength = wavelength_0 / n;
k = 2 * pi / wavelength;

%-----lens type ------
f = 18e-3; %for 10x lens, [mm]
z1 = 198e-3; %[mm], 10x
z2_image = 19.8e-3; %[mm] , 10x
NA=0.3;
%pupil_size_real = 5.4e-3; %[mm]   10x lens D=2fNA= 2x18mmx0.3=10.8mm
pupil_size_real = 3.82e-3; %[mm]   inner rectangular, 10x

%f = 9e-3; %for 20x lens, [mm]
%z1 = 189e-3; %[mm] , 20x
%z2_image = 9.45e-3; %[mm]  , 20x
%NA=0.45;
%pupil_size_real = 2.86e-3;  %[mm]   inner rectangular, 20x

%---------------------------------------------------------
%object and image coordinates
image_size  = 151;                 % matrix size,,,
origin = (image_size+1)/2;
plane_side_size = 1e-6 * (image_size - 1)/2;     % area -75~75um
plane_resolution = 1e-6;                         % resolution 1um
% u_range & v_range ~ -75um~75um 肮荐嘎苗具窃,,,
u_range = (-1)*plane_side_size: plane_resolution : plane_side_size; 
v_range = (-1)*plane_side_size: plane_resolution : plane_side_size;
U_o = imread('source_ring_150.bmp');
U_o = double(U_o(:,:,1))/255;  % U_o (ox, oy) ~ not image size

%Zstack_I_MU_g = zeros(image_size, image_size, z2_series_number);
%Zstack_I_MH = zeros(image_size, image_size, z2_series_number);
Zstack_I_MU_i = zeros(image_size, image_size, z2_series_number);
Zvertical_I_MU_i = zeros(z2_series_number,image_size);


%% Use Gaussian beam

I_z=I_o(1/sqrt())^2*exp(-2(r^2)/w(z)^2);












%%
%1D test
%{
a=zeros(1,size(u_range));
index_a=1;
for u=u_range
    
    a(index_a)=exp(j*pi*(u^2));
    index_a=index_a+1;
end
a=horzcat(zeros(1,1000),a,zeros(1,1000));
b=fftshift(fft(a));
plot(abs(b))
%}
%%

%{
2D test

clear all
close all

func=zeros(151);
index_f1=1;
index_f2=1;
    for u=u_range
        
        for v=v_range
            func(index_f1,index_f2)=exp(j*pi*wavelength*(z2^2)*(1/z1+1/z2-1/f)*(u^2+v^2));
            index_f2=index_f2+1;
        end
        index_f2=1;
        index_f1=index_f1+1;
    end


%}