Skip to content

Commit cc81ace

Browse files
YiChenLaiYiChenLai
committed
new version
1 parent 06fb61e commit cc81ace

6 files changed

+456
-82
lines changed

‎Main_code.asv

+236
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,236 @@
1+
clc; clear; colordef black; format long; close all;
2+
%% Parameter Setting
3+
%------------------%
4+
lambda = 546.1; % Unit : nm
5+
sk16_schott = 1.62286;
6+
surface_num = 2;
7+
distance = [0.01, 0.01, 0.112384569991885]; % Unit : mm
8+
material = [1, sk16_schott, 1]; % Unit : mm
9+
y_radius = [inf, -0.07]; % Unit : mm
10+
aperture = 0.05; % Unit : mm
11+
%------------------%
12+
%------------------%
13+
ang_x = 0;
14+
ang_y = 0;
15+
cross_diameter_num = 501;
16+
%------------------%
17+
%------------------%
18+
Use_Paraxial_Solve = 0; % 0 = No, 1 = Yes
19+
%------------------%
20+
View_Lens = 1; % 0 = No, 1 = Yes
21+
display_line = 21;
22+
Spot_Diagram = 1; % 0 = No, 1 = Yes
23+
display_num__per_
24+
Transmission_Plane = 1; % 0 = No, 1 = Yes
25+
Point_Spread_Function = [1, 1]; % 0 = No, 1 = Yes, perspective:[yz, xy]
26+
Line_Spread_Function = 1; % 0 = No, 1 = Yes
27+
28+
29+
%% Source Setting
30+
lambda = lambda*1e-6; % nm -> mm
31+
[s_x, s_y, s_z, L, M, N] = light_source_setting(aperture,distance,cross_diameter_num,ang_x,ang_y);
32+
33+
%% Calculate Paraxial Focal Length
34+
[BFL, EFL] = paraxial_focal_length(surface_num,distance,material,y_radius);
35+
disp(['BFL = ',num2str(BFL),', EFL = ',num2str(EFL)])
36+
if Use_Paraxial_Solve == 1
37+
distance(end+1) = BFL-distance(end);
38+
material(end+1) = 1;
39+
y_radius(end+1) = inf;
40+
end
41+
42+
%% Ray Tracing
43+
curvature = 1./y_radius;
44+
s_x_all = cell(1,numel(distance)); s_y_all = cell(1,numel(distance)); s_z_all = cell(1,numel(distance));
45+
delta = zeros(size(s_x,1),size(s_x,2));
46+
Opti_Path_Diff = zeros(size(s_x,1),size(s_x,2));
47+
48+
for i = 1:numel(distance)
49+
if i == numel(distance)
50+
z0 = ones(size(z0,1),size(z0,2))*sum(distance);
51+
x0 = s_x+(L./N).*(z0-s_z);
52+
y0 = s_y+(M./N).*(z0-s_z);
53+
54+
x = [s_x;x0]; y = [s_y;y0]; z = [s_z;z0];
55+
s_x_all{i} = x; s_y_all{i} = y; s_z_all{i} = z;
56+
else
57+
z0 = s_z+distance(i)-delta;
58+
x0 = s_x+(L./N).*(z0-s_z);
59+
y0 = s_y+(M./N).*(z0-s_z);
60+
61+
B = N-curvature(i).*(L.*x0+M.*y0);
62+
C = curvature(i).*(x0.^2+y0.^2);
63+
delta = C./(B+sqrt(B.^2-curvature(i).*C));
64+
65+
if i <= surface_num
66+
Opti_Path_Diff = Opti_Path_Diff+abs(delta);
67+
end
68+
69+
x1 = x0+L.*delta; y1 = y0+M.*delta; z1 = z0+N.*delta;
70+
x = [s_x;x1]; y = [s_y;y1]; z = [s_z;z1];
71+
s_x_all{i} = x; s_y_all{i} = y; s_z_all{i} = z;
72+
73+
CosInc = sqrt(B.^2-curvature(i).*C);
74+
nTrans_CosTrans = sqrt((material(i+1).^2)-((material(i).^2).*(1-CosInc.^2)));
75+
k = curvature(i).*(nTrans_CosTrans-material(i).*CosInc);
76+
77+
L_Trans = (material(i).*L-k.*x1)./material(i+1); L = L_Trans;
78+
M_Trans = (material(i).*M-k.*y1)./material(i+1); M = M_Trans;
79+
N_Trans = sqrt(1-(L_Trans.^2+M_Trans.^2)); N = N_Trans;
80+
81+
s_x = x1; s_y = y1; s_z = z1;
82+
end
83+
end
84+
85+
Data = data_reshape(s_x_all,s_y_all,s_z_all,cross_diameter_num);
86+
%% View Lens
87+
if View_Lens == 1
88+
index_y_plane = find(Data.X_1{1}(1,:)==0);
89+
display_line_position = aperture/2*linspace(-1,1,display_line);
90+
index_display_lines = [];
91+
for n = 1:display_line
92+
inx = find(Data.Y_1{1}(1,:)==display_line_position(n));
93+
index_display_lines = [index_display_lines,inx];
94+
end
95+
index = intersect(index_y_plane,index_display_lines);
96+
97+
98+
figure('units','normalized','outerposition',[0 0 1 1])
99+
for n = 1:numel(distance)
100+
if material(n)==1
101+
line_color = 'g';
102+
lin_wid = 0.5;
103+
else
104+
line_color = 'w';
105+
lin_wid = 3;
106+
end
107+
plot(Data.Z_1{n}(:,index)-sum(distance(1:surface_num)),Data.Y_1{n}(:,index), ...
108+
'color',line_color,'linewidth',lin_wid)
109+
hold on
110+
end
111+
axis equal
112+
xlim([-sum(distance(1:surface_num)),sum(distance(surface_num+1:end))])
113+
ylim([-aperture*1.2/2,aperture*1.2/2])
114+
grid on
115+
xlabel('z (mm)')
116+
ylabel('y (mm)')
117+
title('View Lens')
118+
pause(0.01)
119+
end
120+
%% Spot Diagram
121+
if Spot_Diagram == 1
122+
123+
display_line_position = aperture/2*linspace(-1,1,display_line);
124+
index_display_lines = [];
125+
126+
127+
figure
128+
plot(Data.X_1{end}(2,:),Data.Y_1{end}(2,:),'.')
129+
axis equal
130+
title('Spot Diagram')
131+
pause(0.01)
132+
end
133+
%% Optical Path Difference (OPD) at Transmission Plane
134+
trans_plane_data = trans_plane_position_and_optical_path(surface_num, distance, material, Data, L, M, N);
135+
136+
if Transmission_Plane == 1
137+
figure
138+
pcolor(trans_plane_data.x,trans_plane_data.y,trans_plane_data.OP)
139+
axis equal; shading flat; colorbar; colormap('jet')
140+
pause(0.01)
141+
end
142+
143+
%%
144+
if Use_Paraxial_Solve == 1
145+
focal_plane_position = sum(distance(end-1:end));
146+
else
147+
focal_plane_position = distance(end);
148+
end
149+
diffra_limit = diffraction_limit(lambda,aperture,BFL,EFL,focal_plane_position);
150+
%% Line Spread Function
151+
if Line_Spread_Function == 1
152+
LSF_data = line_spread_function(lambda, aperture, trans_plane_data, focal_plane_position);
153+
%%
154+
figure
155+
plot(LSF_data.Monitor_y,LSF_data.Power_normalize,'linewidth',.5,'color',[0.93,0.69,0.13])
156+
hold on
157+
plot(diffra_limit.y,diffra_limit.I,':','linewidth',.5,'color','w')
158+
title('Line Spread Function',['Strehl Ratio = ',num2str(LSF_data.Strehl_ratio)])
159+
grid on
160+
ax = gca;
161+
ax.GridColor = [0.32 0.32 0.32];
162+
163+
end
164+
165+
%% Point Spread Function
166+
if sum(Point_Spread_Function) > 0
167+
PSF_data = point_spread_function(lambda, aperture, trans_plane_data, focal_plane_position, Point_Spread_Function);
168+
%%
169+
if Point_Spread_Function(1)==1
170+
plot_z = PSF_data.Monitor_z-trans_plane_data.dz;
171+
[~, index_z_1] = find(PSF_data.Intensity_normalize==1);
172+
[~, index_z_2] = min(abs(plot_z-focal_plane_position));
173+
disp(['Best Focus Plane = ',num2str(plot_z(index_z_1)),' mm'])
174+
175+
figure
176+
p1 = pcolor(plot_z,PSF_data.Monitor_y,PSF_data.Power_normalize);
177+
p2 = line([plot_z(index_z_1) plot_z(index_z_1)],[PSF_data.Monitor_y(1) PSF_data.Monitor_y(end)], ...
178+
'color','w','linewidth',0.5,'linestyle',':');
179+
p3 = line([plot_z(index_z_2) plot_z(index_z_2)],[PSF_data.Monitor_y(1) PSF_data.Monitor_y(end)], ...
180+
'color','y','linewidth',0.5,'linestyle','-');
181+
title('Point Spread Function YZ',['Best Focus Plane = ',num2str(plot_z(index_z_1)),' mm'])
182+
legend([p2, p3],'Best Focus Plane','Focus Plane','Location','northwest','NumColumns',1)
183+
colormap('jet')
184+
shading interp
185+
xlabel('z (mm)')
186+
ylabel('y (mm)')
187+
colorbar
188+
pause(0.01)
189+
% figure
190+
% plot(PSF_data.Monitor_y,PSF_data.Power_normalize(:,index_z_2),'linewidth',.5,'color',[0.93,0.69,0.13])
191+
% hold on
192+
% title(['Strehl Ratio = ',num2str(PSF_data.Strehl_ratio)])
193+
% grid on
194+
% ax = gca;
195+
% ax.GridColor = [0.32 0.32 0.32];
196+
% pause(0.01)
197+
end
198+
%%
199+
if Point_Spread_Function(2)==1
200+
figure
201+
surf(PSF_data.Monitor_x,PSF_data.Monitor_y,PSF_data.Power_normalize_xy)
202+
title('Point Spread Function XY',['Strehl Ratio = ',num2str(PSF_data.Strehl_ratio_xy)])
203+
colormap('jet')
204+
shading interp
205+
xlabel('x (mm)')
206+
ylabel('y (mm)')
207+
colorbar
208+
pause(0.01)
209+
end
210+
end
211+
212+
%% MTF
213+
LSF = LSF_data.Power;
214+
% LSF = PSF_data.Power(:,index_z_2);
215+
OTF = fftshift(fft(LSF));
216+
MTF = abs(OTF);
217+
MTF = MTF./max(MTF);
218+
size = length(MTF);
219+
D = aperture;
220+
a = linspace(-size/2,size/2,size);
221+
f = a*(1/D);
222+
223+
figure
224+
plot(diffra_limit.f,diffra_limit.MTF',':','linewidth',.5,'color','w')
225+
hold on
226+
plot(f,MTF,'linewidth',.5,'color',[0.93,0.69,0.13])
227+
xlim([0,diffra_limit.cutoff_freq])
228+
title('Diffraction MTF')
229+
xlabel('cycles / mm')
230+
ylabel('Modulation')
231+
grid on
232+
legend('Diffraction Limit',['Ang ',num2str(ang_y),' degree'])
233+
ax = gca;
234+
ax.GridColor = [0.32 0.32 0.32];
235+
236+

0 commit comments

Comments
 (0)