%init
clc; clear 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 갯수맞춰야함,,,
x_range = (-1)*plane_side_size: plane_resolution : plane_side_size;
y_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);

%%

%{
% for specific z2----------------------z2z2z2z2z2z2z2z2z2z2z2z2z2z2z2-----
index_z2 = 1;
for z2 = z2_under_image + z2_image : z2_series_step : z2_upper_image + z2_image

%z2_variation = 0; %[um]
%z2 = z2_image + z2_variation;
M = (-1)*z2/z1;
pupil_size = pupil_size_real / (wavelength * z2); %실렌즈적분구간 for h
MU_g = (-1)/M *U_o; % MU_g (u,v) ~ image size
MU_i = zeros(image_size, image_size); % MU_i (u,v) ~ image size
MH = zeros(image_size, image_size); % MH(u,v) ~ image size

% lens plane resolution 1um level씩,,, in real lens plane,,,
% 계산값이므로 구간 제약없음
lens_plane_resolution_real = 10e-6;
lens_plane_resolution = lens_plane_resolution_real / (wavelength * z2);
x_range = (-1)*pupil_size : lens_plane_resolution : pupil_size;
y_range = (-1)*pupil_size : lens_plane_resolution : pupil_size;

%-----------------------------------------------------------
% Approach #1 : too slow...
% h = ss { h_1<f(x,y)> * h_2<f(u, v, x, y)> } dxdy
% h_1 = exp(j*pi*wavelength*(z2^2)*(1/z1+1/z2-1/f)*(x^2+y^2));
% h_2 = exp(-j*2*pi*(u*x + v*y));

% Approach #2 :
% h = s{h_x<f(u,x)>}dx * s{h_y<f(v,y)>}dy
% h_const=0;

syms x y u v;
h_const = wavelength*(z2^2)*(1/z1+1/z2-1/f);
h_x = exp(j*pi*(h_const*x^2 - 2*u*x));
h_y = exp(j*pi*(h_const*y^2 - 2*v*y));

% constructing h_component(N=image_size)
Mh_components = 0; % s(un)dx or s(vm)dy (they are identical)
index_u = 1;
for u=u_range %Here, u is in real dimension
h_x_temp = 0;
for x=x_range % integration in x, x is also in real dimension (transformed)
% 적분구간 변경 필요!!!! for circular aperture
h_x_temp = h_x_temp + eval(h_x);
end
Mh_components(index_u) = h_x_temp;
index_u = index_u + 1;
end % Now, we have Mh_components
%Skip v=v_range & y=y_range, (identical to u-x)

i_u=0; i_v=0;
% constructing MH(u,v)
for i_u = 1:1:image_size % u,v index scan
for i_v = 1:1:image_size
% assign the calculated h(u,v) = h_x(u) * h_y(v) value into MH matrix
MH(i_u,i_v) = Mh_components(i_u) * Mh_components(i_v);
end
end % end of construction MH(u,v)
%Now we have MH(u,v) (in complex numbers)

%MU_i = conv2(MH,MU_g,'same');

I_MU_g = abs(MU_g).^2;
I_MH = abs(MH).^2;
I_MU_i = conv2(I_MH,I_MU_g,'same');

Zvertical_I_MU_i(index_z2,:)=I_MU_i(origin,:);
%Zstack_I_MU_g(:,:,index_z2) = I_MU_g;
%Zstack_I_MH(:,:,index_z2) = I_MH;
Zstack_I_MU_i(:,:,index_z2) = I_MU_i;
index_z2 = index_z2 + 1;

end % step-wise z2 calculation complete
Zvertical_I_MU_i_N = Zvertical_I_MU_i/(max(max(Zvertical_I_MU_i)));
code_run = 1

%}
%%

%approach 3
%calculate h first
size_index_z2=(z2_upper_image-z2_under_image)/z2_series_step;
%final_h=zeros(image_size,image_size,size_index_z2);

index_z2=1;
%aa=zeros(151);
% lens plane resolution 1um level씩,,, in real lens plane,,,
% 계산값이므로 구간 제약없음

grid_size=383;%이부분이 렌즈 바뀔때 바뀜
fft_size=grid_size;
%length_xi=383
final_h=zeros(grid_size,1+(abs(z2_under_image)+abs(z2_upper_image))/z2_series_step);
final_image=zeros(grid_size,(abs(z2_under_image)+abs(z2_upper_image))/z2_series_step);
%383은 length(xi_range) 대신하여
%% define xi and eta range
syms x y u v;
P=U_o;
%half_pupil_size = (1/2)*pupil_size_real / (wavelength * z2); %실렌즈적분구간 for h
half_pupil_size = (1/2)*pupil_size_real / (wavelength * z2_image); %실렌즈적분구간 for h
lens_plane_resolution_real = 10e-6;
%lens_plane_resolution = lens_plane_resolution_real / (wavelength * z2);
% z2변하면서 xi_range의 길이가 달리지므로 그냥 z2=z2_image로 고정하고 쓰자
lens_plane_resolution = lens_plane_resolution_real / (wavelength * z2_image);
xi_range = (-1)*half_pupil_size : lens_plane_resolution : half_pupil_size;
eta_range = (-1)*half_pupil_size : lens_plane_resolution : half_pupil_size;

%%

for z2 = z2_under_image + z2_image : z2_series_step : z2_upper_image + z2_image

index_f1=1;
index_f2=1;
%final_h=zeros(length(xi_range));%여따 쓰니깐 자꾸 z2마다 초기화 되어버림
func=zeros(length(xi_range),1);%383=length(xi_range)

for xi=xi_range

%for eta=eta_range

func(index_f1)=exp(j*pi*wavelength*(z2^2)*(1/z1+1/z2-1/f)*(xi^2));
%func(index_f1,index_f2)=exp(1/(wavelength*z2^2)*j*pi*(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

%take circurlar pupil into account
asdf=ones(length(func),1);
% for as=1:length(func)
% for df=1:length(func)
% if abs((as-length(func)/2)^2+(df-length(func)/2)^2)>(length(func)/2)^2
% asdf(as,df)=0;
% else
% end
% end
% end

func=func.*asdf;
%aa=fft2(P.*func);

%% zeropadding part
padded_func=zeros(fft_size,1);

padded_func(1:length(func))=func;

aa=fft(padded_func);
%% non padded part
% aa=fft2(func);
%%

aa=fftshift(aa);

%h(x,y,z) at specific z point
final_h(:,index_z2)=aa;

%normalize h
final_h(:,index_z2)=abs(final_h(:,index_z2))/sum(abs(final_h(:,index_z2)));

%final_image(:,:,index_z2)=conv2(aa,U_o,'same');

index_z2
index_z2=index_z2+1;

%break

end

figure;%display sliced h image
sliced_h=zeros(grid_size,601);
sliced_h=final_h;
imagesc(abs(sliced_h));
colormap(gray);

figure;
imagesc(abs(aa));

figure;
bb=zeros(grid_size,1);
bb=final_image(:,1);
imagesc(bb);
colormap(gray);

%%
%-------------------ABSORPTION --------------------------------------
%NOW, we have all data with light diffraction only, (no absorption)
% I(z) = I(0) exp (-z/z_c)
PI=1014; %concentration of photoinitiator [mol/m^3], Darocur 1173 - from doyle (PEGDA1023, ETPTA1014)
molar_extinction_coeff = 1.6; % of photoinitiator [m^3/(mol m)], Darocur 1173 - from doyle
z_c = 1/(PI*molar_extinction_coeff); % = 0.0037 = 3.7mm
% in our case 1/1600 ~ 600um

% 위쪽부터 빛이 아래로 갈때
% z_ab_start = 301; % set z=zo point (no absorption at this point)
% Zvertical_I_MU_i_N_AB = Zvertical_I_MU_i_N;
% for ab_index = (z_ab_start + 1) : 1 : z2_series_number
% Zvertical_I_MU_i_N_AB(ab_index,:) = exp(-1e-6 * (ab_index - z_ab_start)/z_c) * Zvertical_I_MU_i_N_AB(ab_index,:);
% end

% 아래쪽부터 빛이 위로 갈때 (보기 편함)
z_ab_start = 300; % set z=zo point (no absorption at this point)
Zvertical_I_MU_i_N_AB = Zvertical_I_MU_i_N;
for ab_index = (z_ab_start - 1) : -1 : 1
Zvertical_I_MU_i_N_AB(ab_index,:) = exp(-1e-6 * (z_ab_start - ab_index)/z_c) * Zvertical_I_MU_i_N_AB(ab_index,:);
end
%---------------------------------------------------------------------
%--------------------------on/off operation -------------------------

onoff_threshold = 0.103; % 높을수록 1이 적어져,,,
%for onoff_threshold = 0.12:0.002:0.13

Zvertical_I_MU_i_N_AB_onoff = Zvertical_I_MU_i_N_AB;
onoff_index_x = 1; onoff_index_y = 1;
for onoff_index_x = 1 : z2_series_number
for onoff_index_y = 1 : image_size
if Zvertical_I_MU_i_N_AB_onoff(onoff_index_x,onoff_index_y) > onoff_threshold
Zvertical_I_MU_i_N_AB_onoff(onoff_index_x,onoff_index_y) = 255;
else Zvertical_I_MU_i_N_AB_onoff(onoff_index_x,onoff_index_y) = 0;
end
end
end % on-off operation
% cut the range of z > z_ab_start (빛이 아래서오니까, zindex가 큰거 죽여,,,)
Zvertical_I_MU_i_N_AB_onoff((z_ab_start+1):z2_series_number,:) = 0;
% show guideline to indicate center line
Zvertical_I_MU_i_N_AB_onoff((z2_series_number+1)/2,:) = 255;
imshow(Zvertical_I_MU_i_N_AB_onoff)

%figure;
%end

% ----------------------------display functions---------------------------
mesh(x_range,y_range,I_MH)
imshow(I_MH, [min(min(I_MH)) max(max(I_MH))])
imshow(I_MU_i, [min(min(I_MU_i)) max(max(I_MU_i))])
imshow(Zvertical_I_MU_i_N, [min(min(Zvertical_I_MU_i_N(:,:))) max(max(Zvertical_I_MU_i_N(:,:)))])
min_compensate_factor = 0.0;
max_compensate_factor = 0.0;
imshow(Zvertical_I_MU_i_N, [min(min(Zvertical_I_MU_i_N(:,:)))+min_compensate_factor max(max(Zvertical_I_MU_i_N(:,:)))+max_compensate_factor])

imshow(Zvertical_I_MU_i_N_AB, [min(min(Zvertical_I_MU_i_N_AB(:,:))) max(max(Zvertical_I_MU_i_N_AB(:,:)))])
imshow(Zvertical_I_MU_i_N_AB, [min(min(Zvertical_I_MU_i_N_AB(:,:)))+min_compensate_factor max(max(Zvertical_I_MU_i_N_AB(:,:)))+max_compensate_factor])
imshow(Zvertical_I_MU_i_N_AB_onoff)

for z_display_count = 300:25:425
figure;
imshow(1.5*Zstack_I_MU_i(:,:,z_display_count), [min(min(min(Zstack_I_MU_i(:,:,:))))+min_factor max(max(max(Zstack_I_MU_i(:,:,:))))+min_factor])
z_display_count
end
%--------------------------display functions---------------------------