非成像光学设计作业

5个算法的作业。

  1. CPC集光器设计
  2. CPC配光系统设计
  3. 菲涅尔透镜设计
  4. 菲涅尔透镜配光系统
  5. 剪裁法
  6. 椭流线法
  7. 同步多表面设计法

CPC集光器设计

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
%根据抛物线定义,令焦距f=p/2
%令焦点在原点,抛物线方程x^2=2*p*(y+p/2)
clear all;
z=zeros(1000,3);
syms x y
p=2;
y=x.^2/(2*p)-p/2;

%光线入射角度alpha
alpha=pi/6;
%旋转矩阵
rotate_matrix=[cos(alpha),-sin(alpha);
sin(alpha),cos(alpha)];
%旋转抛物线
Muxian_origin=rotate_matrix*[x;y]
%画出旋转后抛物线的图,母线是其中的一部分
hold on
fplot(Muxian_origin(1),Muxian_origin(2),[-5,5])
ax = gca;
ax.XAxisLocation = 'origin';
ax.YAxisLocation = 'origin';
xlim([-10,10]);
axis equal;

%------------------------------------------------------------
%解出选择的母线的原始坐标x的区间,然后带入到抛物线表达式中,标出有用线段
% -----------------------------------------------------------
%求出抛物线与x轴交点得到x1_min
x_jiaodian=double(solve(Muxian_origin(2)));
x_min_jie=x_jiaodian(x_jiaodian>0)+eps(1);
x=x_min_jie;
x1_min=eval(Muxian_origin(1));
%求出母线最大的x值
x_max_jie=double(solve(diff(Muxian_origin(1))));
x=x_max_jie;
x1_max=eval(Muxian_origin(1));
%标出有用线段
fplot(Muxian_origin(1),Muxian_origin(2),[x_min_jie,x_max_jie],'r');
fplot(@(x2) tan(alpha+pi/2)*x2,"k--");
xline(x1_max,"k--");
hold off
%取1000个点保存到z中
x=linspace(x_min_jie,x_max_jie,1000);
z(:,1:2)=[eval(Muxian_origin(1))',eval(Muxian_origin(2))'];

CPC配光系统设计

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
%%初始化参数
%------------------------
%发光LED直径5mm
%高度3m
%光斑直径2m
clear all;
d=5;
D=2000;
H=3000;
A = [-d/2,0];
B = [d/2,0];
E = [-D/2,H];
F = [D/2,H];
%发射角beta
beta=atand((D+d)/2/H);
% alpha为对称轴与x轴的夹角
alpha =90+beta;
%------------------------


%%定义句柄函数
%------------------------
%抛物线对称轴
u = [cosd(alpha),sind(alpha)];
%归一化
nrm = @(v) v/norm(v);
%两点之间距离
P1P2 = @(P2,P1) sqrt(dot(P2-P1,P2-P1));
%两矢量之间的夹角
ang = @(u,v) acosd(dot(nrm(u),nrm(v)));
%抛物线x,y的坐标
L1 = @(P,F,a,phi) (P1P2(P,F)-dot(P-F,[cosd(a),sind(a)]))./(1-cosd(phi))...
.*cosd(phi+a) + F(1);
L2 = @(P,F,a,phi) (P1P2(P,F)-dot(P-F,[cosd(a),sind(a)]))./(1-cosd(phi))...
.*sind(phi+a) + F(2);
%------------------------

%%确定phi的上下限
%------------------------
phi_begin = 360-ang(B-A,u);
phi_end = 360-ang(B-A,u)+ang(B-A,F-A);
phi = phi_begin:0.01:phi_end;
%------------------------

%%抛物线x,y的坐标
%------------------------
%抛物线x,坐标
L_1 =L1(B,A,alpha,phi);
%抛物线y坐标
L_2 =L2(B,A,alpha,phi);
%------------------------

%%画图
%------------------------
plot(L_1,L_2)
hold on
plot(A(1),A(2),'*',B(1),B(2),'*')
axis([-10,10,-5,50])
grid on
ax = gca;
ax.XAxisLocation = 'origin';
ax.YAxisLocation = 'origin';
axis equal;
text(A(1)-1,A(2),'A')
text(B(1)+1,B(2),'B')
plot([A(1) F(1)],[A(2) F(2)],'r--')
hold off
%------------------------

%%导出坐标
%------------------------
%抛物线在参考坐标系下的x,y坐标
z = zeros(length(L_1),3);
z(:,1)= L_1'; %x坐标
z(:,2)= L_2'; %y坐标
save('z.txt','z','-ascii')
%------------------------

菲涅尔透镜设计

Fresnel_Lens01.m

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
clear all;
%%初始化参数
%---------------------------------
% 焦距f=200mm
% 菲涅尔透镜折射率n=1.5
% 球冠半径r=20mm,采取纵切方式
% 焦距公式f=n_1/((n-n1)/R1+(n2-n)R2)
% 空气n1=1,可算出平凸透镜f=2R
% 分成n段
f = 200;
n = 1.5;
r = 20;
R = f/2;
n = 5;
Point_per_sec=1000;
%---------------------------------
z = split_sections(f,R,r,n,Point_per_sec)

% 在solidworks中,竖直线段会被拟合成曲线,选择分段导出坐标
for i=1:n
filename = 'z'+string(i)+'.txt';
z_section = z(1+Point_per_sec*(i-1):Point_per_sec*(i),1:3);
save(filename,'z_section','-ascii');
end

% 导出坐标到z.txt
save('z.txt','z','-ascii');

split_sections.m

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
function z=split_sections(f,R,r,n,Point_per_sec)
% f:焦距
% R:平凸透镜曲率半径
% r:球冠半径
% n:分成n段
Delta_r = r/n;
% 圆心坐标
Cir_cen = [0;-sqrt(R^2-r^2)];
%% 圆心到各个分段点的连线与x轴夹角theta
%---------------------------------------
theta_i=zeros(1,n+1);
for i=1:n+1
theta_i(i) = acos((r-(i-1)*Delta_r)/R);
end
%---------------------------------------
%% 初始化分配内存
%---------------------------------------
z = zeros(n*Point_per_sec,3);
x = z(:,1);
y = z(:,2);
y_theta_i = Cir_cen(2)+R*sin(theta_i);
%---------------------------------------
%% 每个section的坐标
for j=1:n
range = linspace(theta_i(j),theta_i(j+1),Point_per_sec);
x(1+Point_per_sec*(j-1):Point_per_sec*j) = Cir_cen(1)+R*cos(range);
y(1+Point_per_sec*(j-1):Point_per_sec*j) = Cir_cen(2)+R*sin(range)-y_theta_i(j);
end
%% 存入矩阵z,画图
%---------------------------------------
z(:,1) = x;
z(:,2) = y;
hold on
plot(x,y)
yline(0,'--')
axis equal
%---------------------------------------
return

菲涅尔透镜配光系统

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
clear all;
%%初始化参数
%----------
% 在距离菲涅尔透镜d=200mm处放置接收面,设此处光斑大小为R_facula=10mm
% 菲涅尔透镜折射率n_fresnel=1.5,空气n1=1
% 球冠半径rho_N=50mm,采取纵切方式,纵切分成N段
d = 200; %接受面距离菲涅尔透镜底的距离
n_fresnel=1.5; %菲涅尔透镜折射率
rho_N = 50; %球冠半径
n = 3; %将光斑划分成n段
N = sum([1:n].^2)+n; %由划分的光斑数,确定纵切菲涅尔透镜为N段
Point_per_sec = 500; %每一段采样点的个数
R_facula = 10; %光斑大小
Delta_R_facula = R_facula/n; % 把要生成的光斑分成n段
Delta_rho_N = rho_N/N; % 菲涅尔透镜纵切分成N段
%----------

%%预分配内存,提高运行速度
%------------------------
k=1; %k是为每一段参数索引
x = zeros(1,N); %菲涅尔透镜上的第N段与x轴相交点的坐标x
x_2 = zeros(1,N); %光斑上细分的点的坐标x_2
Delta_x = zeros(1,N); %菲涅尔透镜上的第N段坐标x与光斑上细分的点的坐标x_2之差
temp_ang = zeros(1,N); %中间变量的角度
theta = zeros(1,N); %x轴上点和圆心的连线与x轴的夹角,由斯涅尔定律求得
Cir_cen_y = zeros(1,N); %每段圆弧对应的圆心的y坐标
R = zeros(1,N); %每段圆弧对应的半径
for i=1:n
for j=1:i^2+1
x_2(k) = (i-(i^2+1-j)/(i^2+1))*Delta_R_facula;
x(k) = Delta_rho_N*k;
Delta_x(k) = x(k)-x_2(k);
temp_ang(k) = atan(d/Delta_x(k));
theta(k)=atan((n_fresnel-sin(temp_ang(k)))/cos(temp_ang(k)));
Cir_cen_y(k) = -x(k)*tan(theta(k));
R(k) = x(k)/cos(theta(k));
k = k+1;
end
end
%------------------------


%% 圆心到各个分段点的连线与x轴夹角theta范围
%---------------------------------------
theta_i=zeros(N,2);
for i=1:N
theta_i(i,1) = acos(x(i)/R(i));
theta_i(i,2) = acos((x(i)-Delta_rho_N)/R(i));
end
%---------------------------------------

%把点存入z中
z = zeros(N*Point_per_sec,3);
x = z(:,1);
y = z(:,2);
theta_range= zeros(Point_per_sec,N);
for i=1:N
theta_range(:,i) = linspace(theta_i(i,2),theta_i(i,1),Point_per_sec);
x(1+Point_per_sec*(i-1):Point_per_sec*i) = R(i)*cos(theta_range(:,i));
y(1+Point_per_sec*(i-1):Point_per_sec*i) = Cir_cen_y(i)+R(i)*sin(theta_range(:,i));
end
z(:,1) = x;
z(:,2) = y;

%画图
hold on
plot(x,y)
axis equal
yline(0,'--')

%分段导出
for i=1:N
filename = 'z'+string(i)+'.txt';
z_section = z(1+Point_per_sec*(i-1):Point_per_sec*(i),1:3);
save(filename,'z_section','-ascii');
end
save('z.txt','z','-ascii');

剪裁法

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
% 裁剪法的LED反射罩的光学设计
% 参数:d=20mm, H=3000mm, r_N=500mm
clear
format short
% 初始化参数
d = 20;
H = 3000;
r_N = 500;
N = 10000; %分成10000个section
r_i = zeros(1,N);
theta_i = zeros(1,N);

% 求接受面上的r_i
S_N = pi*r_N^2/N;
for i = 1:N
if i == 1
r_i(i) = sqrt(S_N/pi);
else
r_i(i) = sqrt(r_i(i-1)^2+S_N/pi);
end
end

% 求解光通量方程
syms theta I_0 theta_1 theta_2
fun = I_0*cos(theta)*sin(theta)
phi_total = int(fun,theta,[0 pi/2]) %phi_total = I_0/2
phi_N = int(fun,theta,[theta_1 theta_2])
relative = phi_N/phi_total % cos(theta_1)^2 - cos(theta_2)^2

% 求天顶角theta_i
% 由于relative = phi_N/phi_total = 1/N
% 可得theta_2 = acos(sqrt(cos(theta_1)^2-1/N))
for i = 1:N
if i == 1
theta_i(i) = acos(sqrt(1-1/N));
else
theta_i(i) = acos(sqrt(cos(theta_i(i-1))^2-1/N));
end
end

% 初始化参数
Px = zeros(1,N);
Py = zeros(1,N);
la = zeros(1,N);
lb = zeros(1,N);
alpha_i = zeros(1,N-1);
Px_0 = 0;
Py_0 = d;




% A,B,P,r均为1x2坐标数组
% dist为两点之间的距离
% alpha为两条线段夹角的一半
dist = @(A,B) norm(B-A);
vec_norm = @(v) v/norm(v);
alpha = @(P,r) acos(dot(-P,r-P)/(dist(P,[0,0])*dist(P,r)))/2;
% 旋转矩阵
rotate_matrix =@(alpha) [cos(alpha),-sin(alpha);
sin(alpha),cos(alpha)];


% Px,Py坐标
for i = 1:N
if i == 1
la(i) = d*tan(theta_i(i));
lb(i) = d/cos(theta_i(i));
Px(i) = Px_0+la(i);
Py(i) = Py_0;
else
% la(i)=lb(i-1)*sin(theta_i(i)-theta_i(i-1))/...
% sin(pi/2-alpha([Px(i-1),Py(i-1)],[r_i(i-1),-H]))
alpha_i(i-1) = alpha([Px(i-1),Py(i-1)],[r_i(i-1),-H]);
la(i) = lb(i-1)*sin(theta_i(i)-theta_i(i-1))/...
sin(pi/2-alpha_i(i-1)+theta_i(i-1)-theta_i(i));
lb(i) = la(i)*sin(pi/2+alpha_i(i-1))/sin(theta_i(i)-theta_i(i-1));
dP = la(i)*...
rotate_matrix(pi/2+alpha_i(i-1))*vec_norm(-[Px(i-1);Py(i-1)]);
Px(i) = Px(i-1)+dP(1);
Py(i) = Py(i-1)+dP(2);
end
end

%画图并导出P点坐标到z.txt
plot([Px_0 Px],[Py_0 Py])
axis image
length = N+1;
z = zeros(length,3);
z(:,1)=[Px_0 Px]';
z(:,2)=[Py_0 Py]';
save('z.txt','z','-ascii')

椭流线法

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
% 椭流线法的LED反射罩的光学设计
% 参数:d=15mm,H=3300mm,r_N=1000mm.
clear all
format short
% 初始化参数
d = 15;
H = 3300;
r_N = 1000;
N = 6000; %分成100个section
N_sub = 5;
r_i = zeros(1,N);
theta_i = zeros(1,N);

% 求接受面上的r_i
S_N = pi*r_N^2/N;
for i = 1:N
if i == 1
r_i(i) = sqrt(S_N/pi);
else
r_i(i) = sqrt(r_i(i-1)^2+S_N/pi);
end
end
r_i = [0,r_i];


% 求天顶角theta_i
% 由于relative = phi_N/phi_total = 1/N
% 可得theta_2 = acos(sqrt(cos(theta_1)^2-1/N))
for i = 1:N
if i == 1
theta_i(i) = acos(sqrt(1-1/N));
else
theta_i(i) = acos(sqrt(cos(theta_i(i-1))^2-1/N));
end
end
theta_i = [0,theta_i];

Q = zeros(N+1,2);
Q(:,2) = -H;
Q(:,1) = r_i';

% theta = linspace(pi/4,pi/4+pi/6,1000);


%%函数
% 两矢量夹角
ang = @(v1,v2) acos(dot(v1/norm(v1),v2/norm(v2)))

% x轴与长轴的夹角,以x的逆时针方向为正方向
% 参数Q为Q点坐标
alpha = @(Q) -ang([1,0],Q)

% P点与长轴的夹角
varphi = @(theta,Q) pi/2-theta-alpha(Q)

% alpha + varphi
alpha_plus_varphi = @(theta) pi/2-theta


%%P坐标
k = zeros(1,N);
f = zeros(1,N+1);
alpha_i = zeros(1,N+1);
varphi_i = zeros(N,N_sub);
alpha_plus_varphi_i = zeros(N,N_sub);
for i=1:N
theta_ii = linspace(theta_i(i),theta_i(i+1),N_sub);
varphi_i(i,:) = varphi(theta_ii,Q(i,:));
alpha_plus_varphi_i(i,:) = alpha_plus_varphi(theta_ii);
end

for i=1:N+1
f(i) = norm(Q(i,:));
alpha_i(i) = alpha(Q(i,:));
end

k = zeros(1,N+1);
k(1) = 2*d+H;
Px = zeros(N*N_sub,1);
Py = zeros(N*N_sub,1);
for i=1:N
for j=1:N_sub
Px((i-1)*N_sub+j) = (k(i)^2-f(i)^2).*...
cos(alpha_plus_varphi_i(i,j))/...
(2*k(i)-2*f(i).*cos(varphi_i(i,j)));

Py((i-1)*N_sub+j) = (k(i)^2-f(i)^2).*...
sin(alpha_plus_varphi_i(i,j))/...
(2*k(i)-2*f(i).*cos(varphi_i(i,j)));
end
k(i+1) = norm([Px(i*N_sub),Py(i*N_sub)])+norm([Px(i*N_sub),Py(i*N_sub)]-Q(i+1,:));
end

%%画图
% 在solidworks中会有点过于临近错误,这里选择导出z的奇数行
plot(Px,Py)
axis image
z = zeros(N*N_sub,3);
z(:,1) = Px;
z(:,2) = Py;
z=z(1:2:end,:);
save('z.txt','z','-ascii');

同步多表面设计法

SMS.gif

(1)SMS_stage1.m

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
%%SMS设计要求
% 光源AB=0.5m
% 光斑EF=1m
% 光源和接收面的距离H=1m
% 透镜折射率n=1.5
% 人为选取点P_0,l_0
% z(P_0)约为2/3的H
% l_0约为H/5到H/4
clear all;
format short
%% 初始化参数
AB = 500;
EF = 1000;
H = 1000;
n_0 = 1;
n_1 = 1.5;
% 尝试选在EF的中点D,同时设AB中点C
A = [-AB/2;H];
B = [AB/2;H];
C = [0;H];
E = [-EF/2;0];
F = [EF/2;0];
D = [0;0];
P_0 = [0;H*76.52/100];
l_0 = H*0.20; %暂定l_0=H/5


%%旋转矩阵
rotate_matrix=@(alpha) [cos(alpha),-sin(alpha);
sin(alpha),cos(alpha)];

%%已知法线和入射光线,求折射光线方向的单位矢量
% 函数文件:refra.m
% 函数说明:
% 输入:1.入射介质折射率2.折射介质折射率...
% 3.入射起点4.入射光终点5.法矢量
% 输出:1.折射光线单位矢量2.折射光线的起点

%%已知光线的起点和单位矢量,求终点
% 函数文件:Point_end.m
% 函数说明:
% 输入:1.单位矢量2.起点3.长度
% 输出:终点坐标

%%已知入射光线和折射光线,求法线
% 函数文件:Normal_vec.m
% 函数说明:
% 输入:1.入射介质折射率2.折射介质折射率...
% 3.入射光起点4.界面交点5.折射光终点
% 输出:1.入射角2.折射角3.法线单位矢量


%%根据光程相等,求l的长度
% 函数文件:Calculate_l.m
% 函数说明:
% 输入:1.入射介质折射率2.折射介质折射率...
% 3.总光程4.整段光线的起点5.第一个界面交点...
% 6.整段光线的终点7.折射光单位矢量
% 输出:1.l的长度2.第二个界面的交点



N_vecP0 = [0;1]
[Refra_unit_vec_P0_to_Q0,Refra_P0_to_Q0_at_P0] = Refra(n_0,n_1,B,P_0,N_vecP0)
Q_0 = Point_end(Refra_unit_vec_P0_to_Q0,Refra_P0_to_Q0_at_P0,l_0)
[theta_1,theta_2,N_vecQ0] = Normal_vec(n_1,n_0,P_0,Q_0,E)

% 光程
Delta_all = n_0*norm(P_0-B)+n_1*norm(Q_0-P_0)+n_0*norm(E-Q_0)

%%以下是从F到Q_0到P_1到A
[Refra_unit_vec_Q0_to_P1,Refra_Q0_to_P1_at_Q0] = Refra(n_0,n_1,F,Q_0,N_vecQ0)
[l1,P_1] = Calculate_l(n_0,n_1,Delta_all,F,Refra_Q0_to_P1_at_Q0,A,Refra_unit_vec_Q0_to_P1)
[theta_3,theta_4,N_vecP1] = Normal_vec(n_1,n_0,Q_0,P_1,A)


%%第二次循环
%%以下是从B到P_1到Q_1到E
[Refra_unit_vec_P1_to_Q1,Refra_P1_to_Q1_at_P1] = Refra(n_0,n_1,B,P_1,N_vecP1)
[l2,Q_1] = Calculate_l(n_0,n_1,Delta_all,B,Refra_P1_to_Q1_at_P1,E,Refra_unit_vec_P1_to_Q1)
[theta_5,theta_6,N_vecQ1] = Normal_vec(n_1,n_0,P_1,Q_1,E)

%%以下是从F到Q_1到P_2到A
[Refra_unit_vec_Q1_to_P2,Refra_Q1_to_P2_at_Q1] = Refra(n_0,n_1,F,Q_1,N_vecQ1)
[l3,P_2] = Calculate_l(n_0,n_1,Delta_all,F,Refra_Q1_to_P2_at_Q1,A,Refra_unit_vec_Q1_to_P2)
[theta_7,theta_8,N_vecP2] = Normal_vec(n_1,n_0,Q_1,P_2,A)

%%第三次循环
%%以下是从B到P_1到Q_1到E
[Refra_unit_vec_P2_to_Q2,Refra_P2_to_Q2_at_P2] = Refra(n_0,n_1,B,P_2,N_vecP2)
[l4,Q_2] = Calculate_l(n_0,n_1,Delta_all,B,Refra_P2_to_Q2_at_P2,E,Refra_unit_vec_P2_to_Q2)
[theta_9,theta_10,N_vecQ2] = Normal_vec(n_1,n_0,P_2,Q_2,E)

%%以下是从F到Q_2到P_3到A
[Refra_unit_vec_Q2_to_P3,Refra_Q2_to_P3_at_Q2] = Refra(n_0,n_1,F,Q_2,N_vecQ2)
[l5,P_3] = Calculate_l(n_0,n_1,Delta_all,F,Refra_Q2_to_P3_at_Q2,A,Refra_unit_vec_Q2_to_P3)
[theta_11,theta_12,N_vecP3] = Normal_vec(n_1,n_0,Q_2,P_3,A)

%%按对称轴对称到另一边
Change=[-1,0;0,1];
Q_0dash = Change*Q_0
Q_1dash = Change*Q_1
Q_2dash = Change*Q_2
P_1dash = Change*P_1
P_2dash = Change*P_2
P_3dash = Change*P_3

N = 100;
hold on

plot(A(1),A(2),"*")
plot(B(1),B(2),"*")
plot(C(1),C(2),"*")
plot(E(1),E(2),"*")
plot(F(1),F(2),"*")
plot(D(1),D(2),"*")
xline(0,"--")
yline(H,'--')
plot(P_0(1),P_0(2),"r*")
plot(Q_0(1),Q_0(2),"b*")
plot(P_1(1),P_1(2),"r*")
plot(Q_1(1),Q_1(2),"b*")
plot(P_2(1),P_2(2),"r*")
plot(Q_2(1),Q_2(2),"b*")
plot(P_3(1),P_3(2),"r*")

plot(Q_0dash(1),Q_0dash(2),"b*")
plot(Q_1dash(1),Q_1dash(2),"b*")
plot(Q_2dash(1),Q_2dash(2),"b*")
plot(P_1dash(1),P_1dash(2),"r*")
plot(P_2dash(1),P_2dash(2),"r*")
plot(P_3dash(1),P_3dash(2),"r*")

axis image

(2)SMS_stage.m

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
%%首先求一下每个Q点对应的与y轴的交点
k_0 = N_vecQ0(2)/N_vecQ0(1);
k_1 = N_vecQ1(2)/N_vecQ1(1);
k_2 = N_vecQ2(2)/N_vecQ2(1);

y_0 =@(x_0) k_0*(x_0-Q_0(1))+Q_0(2);
y_1 =@(x_1) k_1*(x_1-Q_1(1))+Q_1(2);
y_2 =@(x_2) k_2*(x_2-Q_2(1))+Q_2(2);
b_0 = y_0(0);
b_1 = y_1(0);
b_2 = y_2(0);

%%于是确定圆心
c_0 = [0;b_0];
c_1 = [0;b_1];
c_2 = [0;b_2];

r_0 = norm(c_0-Q_0)
r_1 = norm(c_1-Q_1)
r_2 = norm(c_2-Q_2)

N = 100; %插值100个点
ang = @(A,Com,B) acos(dot((A-Com)/norm(A-Com),(B-Com)/norm(B-Com)))

%%先对Q_0到Q_0dash进行插值
theta_Q0_Q0_dash = ang(Q_0,c_0,Q_0dash)
theta_0j = linspace(-pi/2-theta_Q0_Q0_dash/2,-pi/2+theta_Q0_Q0_dash/2,N);
hold on
Q_0j = zeros(2,N);
P_0j = zeros(2,N);
Refra_unit_vec_0j = zeros(2,N);
Refra_Q0j_to_P0j_at_Q0j = zeros(2,N);
N_vecQ0j = zeros(2,N);
l_0j = zeros(1,N);
theta_IN_Refra_0j = zeros(2,N);


for j = 1:N
Q_0j(1,j) = r_0*cos(theta_0j(j))+c_0(1);
Q_0j(2,j) = r_0*sin(theta_0j(j))+c_0(2);
N_vecQ0j(:,j) = (c_0-Q_0j(:,j))/norm(c_0-Q_0j(:,j));
[Refra_unit_vec_0j(:,j),Refra_Q0j_to_P0j_at_Q0j(:,j)] = Refra(n_0,n_1,F,Q_0j(:,j),N_vecQ0j(:,j));
[l_0j(j),P_0j(:,j)] = Calculate_l(n_0,n_1,Delta_all,F,Refra_Q0j_to_P0j_at_Q0j(:,j),A,Refra_unit_vec_0j(:,j));
% [theta_IN_Refra_0j(1,j),theta_IN_Refra_0j(2,j),~] = Normal_vec(n_1,n_0,Q_0j(:,j),P_0j(:,j),A)
plot(Q_0j(1,j),Q_0j(2,j),'k*')
plot(P_0j(1,j),P_0j(2,j),'k*')
end

%%对另外两个区段也进行插值
theta_Q0_Q1 = ang(Q_1,c_1,[0;0])-ang(Q_0,c_0,[0;0]);
theta_1j = linspace(-pi/2-ang(Q_1,c_1,[0;0]),-pi/2-ang(Q_0,c_0,[0;0]),N);
c_1j = zeros(2,N);
Q_1j = zeros(2,N);
P_1j = zeros(2,N);
Refra_unit_vec_1j = zeros(2,N);
Refra_Q1j_to_P1j_at_Q1j = zeros(2,N);
N_vecQ1j = zeros(2,N);
l_1j = zeros(1,N);
theta_IN_Refra_1j = zeros(2,N);
r_1j = zeros(1,N);
for j = 1:N
c_1j(2,j) = b_1-(j-1)*(b_1-b_0)/(N-1);
r_1j(j) = r_1-(j-1)*(r_1-r_0)/(N-1)
Q_1j(1,j) = r_1j(j)*cos(theta_1j(j));
Q_1j(2,j) = r_1j(j)*sin(theta_1j(j))+c_1j(2,j);
N_vecQ1j(:,j) = (c_1j(:,j)-Q_1j(:,j))/norm(c_1j(:,j)-Q_1j(:,j));
[Refra_unit_vec_1j(:,j),Refra_Q1j_to_P1j_at_Q1j(:,j)] = Refra(n_0,n_1,F,Q_1j(:,j),N_vecQ1j(:,j));
[l_1j(j),P_1j(:,j)] = Calculate_l(n_0,n_1,Delta_all,F,Refra_Q1j_to_P1j_at_Q1j(:,j),A,Refra_unit_vec_1j(:,j));
% [theta_IN_Refra_0j(1,j),theta_IN_Refra_0j(2,j),~] = Normal_vec(n_1,n_0,Q_0j(:,j),P_0j(:,j),A)
plot(Q_1j(1,j),Q_1j(2,j),'k*')
plot(P_1j(1,j),P_1j(2,j),'k*')
end

theta_Q1_Q2 = ang(Q_2,c_2,[0;0])-ang(Q_1,c_1,[0;0]);
theta_2j = linspace(-pi/2-ang(Q_2,c_2,[0;0]),-pi/2-ang(Q_1,c_1,[0;0]),N);
c_2j = zeros(2,N);
Q_2j = zeros(2,N);
P_2j = zeros(2,N);
Refra_unit_vec_2j = zeros(2,N);
Refra_Q2j_to_P2j_at_Q2j = zeros(2,N);
N_vecQ2j = zeros(2,N);
l_2j = zeros(1,N);
theta_IN_Refra_2j = zeros(2,N);
r_2j = zeros(1,N);
for j = 1:N
c_2j(2,j) = b_2-(j-1)*(b_2-b_1)/(N-1);
r_2j(j) = r_2-(j-1)*(r_2-r_1)/(N-1)
Q_2j(1,j) = r_2j(j)*cos(theta_2j(j));
Q_2j(2,j) = r_2j(j)*sin(theta_2j(j))+c_2j(2,j);
N_vecQ2j(:,j) = (c_2j(:,j)-Q_2j(:,j))/norm(c_2j(:,j)-Q_2j(:,j));
[Refra_unit_vec_2j(:,j),Refra_Q2j_to_P2j_at_Q2j(:,j)] = Refra(n_0,n_1,F,Q_2j(:,j),N_vecQ2j(:,j));
[l_2j(j),P_2j(:,j)] = Calculate_l(n_0,n_1,Delta_all,F,Refra_Q2j_to_P2j_at_Q2j(:,j),A,Refra_unit_vec_2j(:,j));
% [theta_IN_Refra_0j(1,j),theta_IN_Refra_0j(2,j),~] = Normal_vec(n_1,n_0,Q_0j(:,j),P_0j(:,j),A)
plot(Q_2j(1,j),Q_2j(2,j),'k*')
plot(P_2j(1,j),P_2j(2,j),'k*')
end

%%关于对称轴的对称点
P_0jdash = Change*P_0j;
Q_1jdash = Change*Q_1j;
P_1jdash = Change*P_1j;
Q_2jdash = Change*Q_2j;
P_2jdash = Change*P_2j;

for j = 1:N
plot(Q_1jdash(1,j),Q_1jdash(2,j),'k*')
plot(Q_2jdash(1,j),Q_2jdash(2,j),'k*')
plot(P_0jdash(1,j),P_0jdash(2,j),'k*')
plot(P_1jdash(1,j),P_1jdash(2,j),'k*')
plot(P_2jdash(1,j),P_2jdash(2,j),'k*')
end


xline(0,"--")
yline(H,'--')
axis([-500 500 0 1100])
axis equal

%% 导出点
z_P0 = zeros(N,3);
z_P1 = zeros(N,3);
z_P2 = zeros(N,3);
z_Q0 = zeros(N,3);
z_Q1 = zeros(N,3);
z_Q2 = zeros(N,3);

z_P0(:,1:2) = P_0j';
z_P1(:,1:2) = P_1j';
z_P2(:,1:2) = P_2j';
z_Q0(:,1:2) = Q_0j';
z_Q1(:,1:2) = Q_1j';
z_Q2(:,1:2) = Q_2j';

save('z_P0.txt','z_P0','-ascii')
save('z_P1.txt','z_P1','-ascii')
save('z_P2.txt','z_P2','-ascii')
save('z_Q0.txt','z_Q0','-ascii')
save('z_Q1.txt','z_Q1','-ascii')
save('z_Q2.txt','z_Q2','-ascii')

(3)Refra.m

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
function [unit_vec,In_Point_end] = Refra(In_n,Refra_n,In_Point_start,In_Point_end,Normal_line)
%本函数为已知法线和入射光线,求折射光线方向的单位矢量

% 归一化入射光线和法线
Incident = In_Point_end - In_Point_start;
Incident = Incident/norm(Incident);
Normal_line = Normal_line/norm(Normal_line);
theta = acos(dot(Incident,Normal_line));
% 确定入射角,并确保法矢量方向与折射光线夹角为锐角
disp('未修正的法线')
disp(Normal_line)
disp('未修正入射角')
disp(theta*180/pi)
if theta > pi/2
theta = pi-theta;
Normal_line = -Normal_line;
end
disp('入射角')
disp(theta*180/pi)
disp('修正后与折射方向成锐角的法线')
disp(Normal_line)
% 确定切线方向与折射、入射方向成锐角
tan_line = [0,-1;1,0]*Normal_line;
if dot(tan_line,Incident) < 0
tan_line = -tan_line;
end
% 由snell定律求出折射角
theta_1 = asin(In_n/Refra_n*sin(theta));
disp('折射角')
disp(theta_1*180/pi)

if isreal(theta_1) == 0
disp("total reflection")
unit_vec = 0;
return
end
% 求出折射光线单位矢量
unit_vec = [cos(theta_1),-sin(theta_1);sin(theta_1),cos(theta_1)]*Normal_line;
if dot(unit_vec,tan_line) < 0
unit_vec = [cos(-theta_1),-sin(-theta_1);sin(-theta_1),cos(-theta_1)]*Normal_line;
end
% 检验
disp('入射n_0sin(theta)')
disp(In_n*sin(theta))
disp('折射n1sin(theta_1)')
disp(Refra_n*sin(theta_1))


return
end

(4)Point_end.m

1
2
3
4
5
6
function Point_end = Point_end(unit_vec,Point_start,length)
%本函数已知光线的起点和单位矢量,求终点

Point_end = Point_start+length*unit_vec;
return
end

(5)Normal_vec.m

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
function [Sol_theta_in,Sol_theta_refra,Normal_line] = Normal_vec(In_n,Refra_n,Point_start,Point_intermediate,Point_end)
%已知入射光线和折射光线,求法线
% 入射光线和折射光线归一化
Incident = Point_intermediate - Point_start;
Refraction = Point_end - Point_intermediate;
Incident = Incident/norm(Incident);
Refraction = Refraction/norm(Refraction);
% 入射光线和折射光线的夹角
Ang = acos(dot(Incident,Refraction));
% 求入射角和折射角
syms theta_in theta_refra
if In_n > Refra_n
eqn1 = theta_in - theta_refra == -Ang;
eqn2 = In_n * sin(theta_in) == Refra_n * sin(theta_refra);
[Sol_theta_in,Sol_theta_refra] = vpasolve(eqn1,eqn2,[theta_in,theta_refra],[0 pi]);
else
eqn1 = theta_in - theta_refra == Ang;
eqn2 = In_n * sin(theta_in) == Refra_n * sin(theta_refra);
[Sol_theta_in,Sol_theta_refra] = vpasolve(eqn1,eqn2,[theta_in,theta_refra],[0 pi]);
end
Sol_theta_in = eval(Sol_theta_in)
Sol_theta_refra = eval(Sol_theta_refra);
disp('入射角')
disp(Sol_theta_in)
disp('折射角')
disp(Sol_theta_refra)
disp('全反射临界角')
disp(asin(Refra_n/In_n))
% 求法矢量,以同方向旋转入射光线和折射光线,如果入射光线旋转入射角和折射光线旋转折射角之后,两条光线内积为1,则此旋转后的位置就是法线
% 若不为0,则反向旋转。
Temp1 = [cos(Sol_theta_in),-sin(Sol_theta_in);sin(Sol_theta_in),cos(Sol_theta_in)]*Incident;
Temp2 = [cos(Sol_theta_refra),-sin(Sol_theta_refra);sin(Sol_theta_refra),cos(Sol_theta_refra)]*Refraction;
disp('看看有没有转反,旋转之后两个矢量相等就没问题')
disp(Temp1)
disp(Temp2)
% disp(1-dot(Temp1,Temp2)>=eps(1))
if 1-dot(Temp1,Temp2)>eps(1)
Temp1 = [cos(-Sol_theta_in),-sin(-Sol_theta_in);sin(-Sol_theta_in),cos(-Sol_theta_in)]*Incident;
Temp2 = [cos(-Sol_theta_refra),-sin(-Sol_theta_refra);sin(-Sol_theta_refra),cos(-Sol_theta_refra)]*Refraction;
Ang2 = acos(dot(Temp1,Temp2));
disp('似乎转反了,反向转一次再看看两个矢量是否相等')
disp(Temp1)
disp(Temp2)
end
Normal_line = Temp1;

end

(6)Calculate_l.m

1
2
3
4
5
6
7
8
9
10
11
12
function [Sol_l,P_point] = Calculate_l(In_n,Refra_n,Delta_all,Known_Point1,Known_Point2,Known_Point3,Refra_unit_vec)
% 光程相等
Delta_remain = Delta_all - In_n*norm(Known_Point2-Known_Point1);
syms l
assume(l>0)
P = [Known_Point2(1)+l*Refra_unit_vec(1);Known_Point2(2)+l*Refra_unit_vec(2)];
eqn = Refra_n*norm(l*Refra_unit_vec) + In_n*norm(Known_Point3-P) == Delta_remain;
Sol_l = eval(vpasolve(eqn,[0,Inf]));
l = Sol_l;
P_point = eval(P);
return
end
打赏
  • Copyrights © 2020-2021 haoyu fang
  • 访问人数: | 浏览次数:

请我喝杯啤酒吧~

支付宝
微信