Files
wxchen ba13a3f053 Initial commit
new file:   Copy_of_Exz_with_hamamatsu.m
	new file:   Exz_mod.m
	new file:   Exz_with_hamamatsu.m
	new file:   GS.m
	new file:   MRAF.m
	new file:   MRAF_8bit.bmp
	new file:   PSF of spherical aberration/Exz.m
	new file:   PSF of spherical aberration/Thumbs.db
	new file:   PSF of spherical aberration/josaa-12-10-2136.pdf
	new file:   PSF of spherical aberration/josaa-12-2-325.pdf
	new file:   PSF of spherical aberration/main.m
	new file:   "PSF of spherical aberration/\351\207\215\350\246\201psf\346\226\207\347\214\256.pdf"
	new file:   Spherical_aberration_SiminCao.m
	new file:   gen_m.m
	new file:   gen_rectangle.m
	new file:   hamamatsu.m
	new file:   m.tif
	new file:   m2.tif
	new file:   photo/50-30.BMP
	new file:   photo/8bit_50-30.BMP
	new file:   photo/8bit_ellipse.BMP
	new file:   photo/convert_8bit.exe
	new file:   photo/ellipse.BMP
	new file:   rect_MRAF_SiminCao.m
	new file:   rectangle.tif
	new file:   size/.vscode/settings.json
	new file:   size/black_c_20THSize_4f_1.064lamda.bmp
	new file:   size/black_c_30THSize_4f_61.064lamda.bmp
	new file:   size/black_output.bmp
	new file:   size/black_rect_30THSize_4f_1.064lamda.bmp
	new file:   size/black_rect_30THSize_4f_6_1.064lamda.bmp
	new file:   size/c_20THSize_4f_1.064lamda.bmp
	new file:   size/c_20THSize_4f_1.064lamda_resize.bmp
	new file:   size/c_30THSize_4f_61.064lamda.bmp
	new file:   size/c_30THSize_4f_61.064lamda_resize.bmp
	new file:   size/noisy_output.bmp
	new file:   size/output.bmp
	new file:   size/rect_30THSize_4f_1.064lamda.bmp
	new file:   size/rect_30THSize_4f_1.064lamda_resize.bmp
	new file:   size/rect_30THSize_4f_6_1.064lamda.bmp
	new file:   size/rect_30THSize_4f_6_1.064lamda_resize.bmp
	new file:   size/resize_4.7z
	new file:   size/resize_black.7z
	new file:   size/size copy.py
	new file:   size/size.py
	new file:   size/wave.7z
	new file:   sp.m
	new file:   to8bit.m
	new file:   trans_8bit.zip
	new file:   wavef/A.bmp
	new file:   wavef/B_linear.bmp
	new file:   wavef/PHA SID230828-2003.csv
	new file:   wavef/PHA_bilinear_1280_1024.csv
	new file:   wavef/PHA_bilinear_reversal.csv
	new file:   wavef/PHA_output_1280_1024.csv
	new file:   wavef/Untitled-1.py
	new file:   wavef/filled.bmp
	new file:   wavef/from PIL import Image.py
	new file:   wavef/matrix_filled.csv
	new file:   wavef/output.bmp
	new file:   wavef/output.csv
	new file:   wavef/output2.bmp
	new file:   wavef/pha_wavef copy.py
	new file:   wavef/pha_wavef.py
	new file:   wavef/pha_wavef_step.py
	new file:   wavef/wavef.zip
	new file:   wavef/xy_values.csv
	new file:   "wavef/\346\226\260\345\273\272\346\226\207\344\273\266\345\244\271/matrix.csv"
	new file:   "wavef/\346\226\260\345\273\272\346\226\207\344\273\266\345\244\271/matrix_filled.csv"
2023-08-29 23:06:40 +08:00

88 lines
3.0 KiB
Matlab
Raw Permalink Blame History

This file contains ambiguous Unicode characters

This file contains Unicode characters that might be confused with other characters. If you think that this is intentional, you can safely ignore this warning. Use the Escape button to reveal them.

%% setting
clear;
depthin = zeros(200, 1);
zpsflong = zeros(200, 1);
intensity_record = [];
for jj = 0:600
lambda = 0.8; %
d = jj; %
n1 = 1.08;
n2 = 1.43;
NA = 0.75;
alpha = asin(NA / n1);
k0 = 2 * pi / lambda;
x = 0; %linspace(-6,6,200);
D_z = linspace(-10, 50, 200);
z = D_z * n2 / n1 + (n2 / n1 - 1) * d;
v = 2 * pi * n1 / lambda * sqrt(2 * x .^ 2) * (NA / n1); %
u = 2 * pi * n2 / lambda * z * (NA / n1) ^ 2;
[u, v] = meshgrid(u, v);
%% Fresnel transmission coefficients
Rs = @(phi1, phi2)((n1 * cos(phi1) - n2 * cos(phi2)) / (n1 * cos(phi1) + n2 * cos(phi2)));
Rp = @(phi1, phi2)((n1 * cos(phi2) - n2 * cos(phi1)) / (n1 * cos(phi2) + n2 * cos(phi1)));
Ts = @(phi1, phi2)(1 - Rs(phi1, phi2));
Tp = @(phi1, phi2)(1 - Rp(phi1, phi2));
%% Sphere Aberration
Phi_d = @(phi1, phi2, d)(-d * (n1 * cos(phi1) - n2 * cos(phi2)));
%% I transfored integrals
I0_in = @(phi1, phi2)cos(phi1) ^ 0.5 * sin(phi1) * (Ts(phi1, phi2) + Tp(phi1, phi2) * cos(phi2)) ...
* exp(1i * u .* cos(phi2) / sin(alpha) ^ 2) * exp(1i * k0 * Phi_d(phi1, phi2, d)) .* besselj(0, v * sin(phi1) / sin(alpha));
I1_in = @(phi1, phi2)cos(phi1) ^ 0.5 * sin(phi1) * (Tp(phi1, phi2) * sin(phi2)) ...
* exp(1i * u .* cos(phi2) / sin(alpha) ^ 2) * exp(1i * k0 * Phi_d(phi1, phi2, d)) .* besselj(1, v * sin(phi1) / sin(alpha));
I2_in = @(phi1, phi2)cos(phi1) ^ 0.5 * sin(phi1) * (Ts(phi1, phi2) - Tp(phi1, phi2) * cos(phi2)) ...
* exp(1i * u .* cos(phi2) / sin(alpha) ^ 2) * exp(1i * k0 * Phi_d(phi1, phi2, d)) .* besselj(2, v * sin(phi1) / sin(alpha));
%% integree
intnum = 100;
Dphi = alpha / intnum;
I0 = 0;
I1 = 0;
I2 = 0;
for ii = 0:intnum
disp(ii);
phi1 = ii * Dphi;
phi2 = asin(sin(phi1) * n1 / n2);
I0 = I0 + I0_in(phi1, phi2) * Dphi;
I1 = I1 + I1_in(phi1, phi2) * Dphi;
I2 = I2 + I2_in(phi1, phi2) * Dphi;
end
%%
l0 = 1; %l_0 is an amplitude factor
f = 250; %f<><66><EFBFBD><EFBFBD>
K = pi * n1 * f * l0 / lambda;
%<25><>xz<78><7A>һ<EFBFBD><D2BB><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD> theta_p=0
e_2x = -1i * K * (I2 + I0);
e_2y = 0;
e_2z = -2 * K * I1;
intensity = (e_2x .* conj(e_2x) + e_2y .* conj(e_2y) + e_2z .* conj(e_2z));
intensity_record = [intensity_record; intensity];
z_pix = (max(z) - min(z)) / size(z, 2);
[depthin(jj + 1), maxidx] = max(intensity(:));
intensity_mainlobe = intensity(maxidx - 33:maxidx + 33);
intensity_mainlobe(intensity_mainlobe < max(intensity_mainlobe) / 2) = [];
zpsflong(jj + 1) = size(intensity_mainlobe, 2) * z_pix;
end
figure1 = figure;
ribbon(z, depthin, 0.1);
colormap('autumn');
% Create ylabel
ylabel({'Longitudinal length/um'});
% Create zlabel
zlabel({'intensity'});
legend('50um', '100um', '150um', '200um', '250um', '300um', '350um', '400um', '450um');
% figure;
% plot(intensity(100,:));
% hold on;
% plot(intensity(:,101),'r-');
% figure;
% [c,h]=contour(z,x,intensity/max(intensity(:)),[0.025 0.05 0:0.1:1],'-');
%clabel(c,h);
depthin = depthin / max(depthin);