目录
本文篇幅较长，分为上中下三篇
定日镜场的输出功率优化（附获奖论文及MATLAB代码实现）（中）
定日镜场的输出功率优化（附获奖论文及MATLAB代码实现）（下）
代码实现
问题二代码（SecQues.m）
global W_a;
global W_b;
global h2;
global h1;
global L0;
WAAA=5.5:0.1:6.5;
WBBB=5.5:0.1:6.5;
R_year=[];
E_year=[];
for W_a=WAAA
for W_b=WBBB
if W_a<W_b
continue;
else
h2_range=linspace(W_b/2,6, 5);
for h2=h2_range
h1=80;
L0_range=linspace(0,350, 5);
for L0=L0_range
%L0=0;
%% 入射角
%% 计算太阳赤纬角 delta
% 输入：D 为以春分作为第 0 天起算的天数
%D = 31; % 4 月 21 日对应的天数
%%输入 D 得到 R
% R=Main(D);
D=[-59 -28 0 31 61 92 122 153 184 214 245 275];
R_total=[];
EEEE_total=[];
for i=1:1:size(D,1)
disp(i);
t = 1 - c;
x = axis(1);
y = axis(2);
z = axis(3);
RR = [
t*x*x + c, t*x*y - s*z, t*x*z + s*y;
t*x*y + s*z, t*y*y + c, t*y*z - s*x;
t*x*z - s*y, t*y*z + s*x, t*z*z + c
];
End
程序二 问题二代码（SecQues.m）
global W_a;
global W_b;
global h2;
global h1;
global L0;
WAAA=5.5:0.1:6.5;
WBBB=5.5:0.1:6.5;
R_year=[];
E_year=[];
for W_a=WAAA
for W_b=WBBB
if W_a<W_b
continue;
else
h2_range=linspace(W_b/2,6, 5);
for h2=h2_range
h1=80;
L0_range=linspace(0,350, 5);
for L0=L0_range
%L0=0;
%% 入射角
%% 计算太阳赤纬角 delta
% 输入：D 为以春分作为第 0 天起算的天数
%D = 31; % 4 月 21 日对应的天数
%%输入 D 得到 R
% R=Main(D);
D=[-59 -28 0 31 61 92 122 153 184 214 245 275];
R_total=[];
EEEE_total=[];
for i=1:1:size(D,1)
disp(i);
[R,EEEE]=Main(D(i));
R_total=[R_total;R];
EEEE_total=[EEEE_total;EEEE];
end
R_year=[R_year;mean(R_total)];
E_year=[E_year;mean(EEEE_total) h2 L0];
end
end
end
end
end
function [R,EEEE]=Main(Day)
global W_a;
global W_b;
global h2;
global h1;
global L0;
EEEE=[];
D=Day;
STT=[9 10.5 12 13.5 15]; %当地时间,假设为 15 点
result=[];%日结果详细矩阵 5*4 5 个时间点 4 个参数 R 才是最终值
eta_ref=0.92;%镜面反射率
O_R=[0 0 h1];%集热器中心坐标
%定日镜选取根据与中心的距离进行排序
% 计算 data 的程序
data=data_func(); %经过排序后的定日镜坐标 data---1745*3 x,y,r
% 计算太阳赤纬角 delta
sin_delta = sin(2*pi*D/365) * sin(2*pi/360 * 23.45);
delta = asin(sin_delta); % 得到赤纬角的弧度值
phi = deg2rad(39.4); % 北纬 39.4 度，转换为弧度
for XH_ST=1:size(STT,2)
ST=STT(XH_ST);
omega =pi/12*(ST-12); % 声明 omega 的值
% 计算 sin(alpha_s)
sin_alpha_s = cos(delta) * cos(phi) * cos(omega) + sin(delta) * sin(phi);
G_0=1.366; %太阳常数
Hh=3;%3km
ha=0.4237 - 0.00821*(6 - Hh)^2;%大气压
hb=0.5055 + 0.00595*(6.5 - Hh)^2;%大气压
hc=0.2711 + 0.01858*(2.5 - Hh)^2;%大气压
DNI = G_0*(ha + hb*exp(-hc/sin_alpha_s)); %直接法辐照度
alpha_s = asin(sin_alpha_s); % 得到 alpha_s 的弧度值
% 根据公式计算 cos(gamma_s)
cos_gamma_s = (sin(delta) - sin_alpha_s * sin(phi)) / (cos(alpha_s) * cos(phi));
[R,EEEE]=Main(D(i));
R_total=[R_total;R];
EEEE_total=[EEEE_total;EEEE];
end
R_year=[R_year;mean(R_total)];
E_year=[E_year;mean(EEEE_total) h2 L0];
end
end
end
end
end
function [R,EEEE]=Main(Day)
global W_a;
global W_b;
global h2;
global h1;
global L0;
EEEE=[];
D=Day;
STT=[9 10.5 12 13.5 15]; %当地时间,假设为 15 点
result=[];%日结果详细矩阵 5*4 5 个时间点 4 个参数 R 才是最终值
eta_ref=0.92;%镜面反射率
O_R=[0 0 h1];%集热器中心坐标
%定日镜选取根据与中心的距离进行排序
% 计算 data 的程序
data=data_func(); %经过排序后的定日镜坐标 data---1745*3 x,y,r
% 计算太阳赤纬角 delta
sin_delta = sin(2*pi*D/365) * sin(2*pi/360 * 23.45);
delta = asin(sin_delta); % 得到赤纬角的弧度值
phi = deg2rad(39.4); % 北纬 39.4 度，转换为弧度
for XH_ST=1:size(STT,2)
ST=STT(XH_ST);
omega =pi/12*(ST-12); % 声明 omega 的值
% 计算 sin(alpha_s)
sin_alpha_s = cos(delta) * cos(phi) * cos(omega) + sin(delta) * sin(phi);
G_0=1.366; %太阳常数
Hh=3;%3km
ha=0.4237 - 0.00821*(6 - Hh)^2;%大气压
hb=0.5055 + 0.00595*(6.5 - Hh)^2;%大气压
hc=0.2711 + 0.01858*(2.5 - Hh)^2;%大气压
DNI = G_0*(ha + hb*exp(-hc/sin_alpha_s)); %直接法辐照度
alpha_s = asin(sin_alpha_s); % 得到 alpha_s 的弧度值
% 根据公式计算 cos(gamma_s)
cos_gamma_s = (sin(delta) - sin_alpha_s * sin(phi)) / (cos(alpha_s) * cos(phi));
if cos_gamma_s > 1
cos_gamma_s = 1;
elseif cos_gamma_s < -1
cos_gamma_s = -1;
end
gamma_s = acos(cos_gamma_s); % 得到 gamma_s 的弧度值
x=cos(alpha_s)*sin(gamma_s);
y=cos(alpha_s)*cos(gamma_s);
z=sin(alpha_s);
S_i=[x y z];%太阳入射方向单位向量的相反
Store=zeros(size(data,1),4);%光学效率 阴影遮挡效率 余弦效率 截断效率
E_field=0;%定日镜场的输出热功率
for num=1:1:size(data,1)%%定日镜编号
disp(num);
O_A=[data(num,1) data(num,2) h2];%定日镜中心坐标
%计算集热器截断效率 eta_trunc
d_HR=sqrt((O_A(1)-O_R(1))^2+(O_A(2)-O_R(2))^2+(O_A(3)-O_R(3))^2);%定日镜镜面中心与集
热器中心的距离
eta_at = 0.99321 - 0.0001176*d_HR + 1.97e-8 * (d_HR^2); % 计算吸收塔的吸收率
S_AR=(O_R-O_A)/norm(O_R-O_A);%定日镜中心到集热器中心的单位向量 %% 反射角
%% 计算法向量
S_n=(S_i+S_AR)/norm(S_i+S_AR);%法向量单位向量
% 计算方位角 Azimuth (θ)
A_H = atan2(S_n(1), S_n(2));
% 计算俯仰角 Elevation (ϕ)
E_H = acos(S_n(3));
%
T=[-sin(E_H),-sin(A_H)*cos(E_H),cos(A_H)*cos(E_H);cos(E_H),-sin(A_H)*sin(E_H),cos(A_H)*sin(E_H
);0,cos(A_H),sin(A_H)];
M1 = [1, 0, 0; 0, cos(E_H), sin(E_H); 0, -sin(E_H), cos(E_H)];
M2 = [cos(pi-A_H), sin(pi-A_H), 0; -sin(pi-A_H), cos(pi-A_H), 0; 0, 0, 1];
% 计算矩阵乘积
T = M1 * M2; %镜坐标系变为地坐标系的转换矩阵
S_An=[0 0 1];%镜面 A 坐标系上的的法向量
SSS=S_An*T;
tolerance = 1e-10;
difference = abs(SSS - S_n);
if all(difference(:) >tolerance)
error('法向量转换错误');
end
%%
%计算余弦效率 eta_cos
eta_cos = dot(S_i, S_n) / (norm(S_i) * norm(S_n)); % 计算入射向量相反和法向量之间的夹
角的余弦值
if cos_gamma_s > 1
cos_gamma_s = 1;
elseif cos_gamma_s < -1
cos_gamma_s = -1;
end
gamma_s = acos(cos_gamma_s); % 得到 gamma_s 的弧度值
x=cos(alpha_s)*sin(gamma_s);
y=cos(alpha_s)*cos(gamma_s);
z=sin(alpha_s);
S_i=[x y z];%太阳入射方向单位向量的相反
Store=zeros(size(data,1),4);%光学效率 阴影遮挡效率 余弦效率 截断效率
E_field=0;%定日镜场的输出热功率
for num=1:1:size(data,1)%%定日镜编号
disp(num);
O_A=[data(num,1) data(num,2) h2];%定日镜中心坐标
%计算集热器截断效率 eta_trunc
d_HR=sqrt((O_A(1)-O_R(1))^2+(O_A(2)-O_R(2))^2+(O_A(3)-O_R(3))^2);%定日镜镜面中心与集
热器中心的距离
eta_at = 0.99321 - 0.0001176*d_HR + 1.97e-8 * (d_HR^2); % 计算吸收塔的吸收率
S_AR=(O_R-O_A)/norm(O_R-O_A);%定日镜中心到集热器中心的单位向量 %% 反射角
%% 计算法向量
S_n=(S_i+S_AR)/norm(S_i+S_AR);%法向量单位向量
% 计算方位角 Azimuth (θ)
A_H = atan2(S_n(1), S_n(2));
% 计算俯仰角 Elevation (ϕ)
E_H = acos(S_n(3));
%
T=[-sin(E_H),-sin(A_H)*cos(E_H),cos(A_H)*cos(E_H);cos(E_H),-sin(A_H)*sin(E_H),cos(A_H)*sin(E_H
);0,cos(A_H),sin(A_H)];
M1 = [1, 0, 0; 0, cos(E_H), sin(E_H); 0, -sin(E_H), cos(E_H)];
M2 = [cos(pi-A_H), sin(pi-A_H), 0; -sin(pi-A_H), cos(pi-A_H), 0; 0, 0, 1];
% 计算矩阵乘积
T = M1 * M2; %镜坐标系变为地坐标系的转换矩阵
S_An=[0 0 1];%镜面 A 坐标系上的的法向量
SSS=S_An*T;
tolerance = 1e-10;
difference = abs(SSS - S_n);
if all(difference(:) >tolerance)
error('法向量转换错误');
end
%%
%计算余弦效率 eta_cos
eta_cos = dot(S_i, S_n) / (norm(S_i) * norm(S_n)); % 计算入射向量相反和法向量之间的夹
角的余弦值
%%%如果在塔挡阴影范围内 eta_sb=0；
%计算阴影遮挡效率 eta_sb
%step1：塔挡损失 定日镜镜面中心是否在阴影中,如果是 eta_sb=0 %i 是-S_i
[P_top_1_intersect,P_top_2_intersect,P_bot_1_intersect,P_bot_2_intersect] =
shadow_range(-S_i, h2) ;%i 是-S_i
isInside = pointInRectangle(O_A(1), O_A(2), P_top_1_intersect, P_top_2_intersect,
P_bot_1_intersect, P_bot_2_intersect);
if isInside %在里面
eta_sb=0;
eta_trunc=0;
Store(num,1)=0;
Store(num,2)=0;
Store(num,3)=eta_cos;
Store(num,4)=0;
else %不在里面
R = 100;
X = linspace(-1/2 * W_a,1/2 * W_a, 5);
Y = linspace(-1/2 * W_b,1/2 *W_b, 5);
% Assuming S_i, S_AR, and T are defined elsewhere
T_inv = inv(T);
% 在指定的半径内寻找定日镜
distances = sqrt((data(1:num-1, 1) - O_A(1)).^2 + (data(1:num-1, 2) - O_A(2)).^2);
valid_indices = find(distances <= R);
% Main loops
valid_point_count = 0;%有效点
valid_points = 0;%每个点的有效光线
obstructed_mirrors = [];
for x = X
for y = Y
dot_A = [x, y, 0];
dot_AG = dot_A * T + O_A;
is_valid = true;
for idx = valid_indices'
O_B = [data(idx, 1), data(idx, 2), h2];
dot_B = (dot_AG - O_B)*T_inv;
S_iB = S_i*T_inv;
H2 = [(S_iB(3)*dot_B(1) - S_iB(1)*dot_B(3))/S_iB(3),
(S_iB(3)*dot_B(2) - S_iB(2)*dot_B(3))/S_iB(3), 0];
if ~(-1/2*W_a <= H2(1) && H2(1) <= 1/2*W_a && -1/2*W_b <= H2(2) && H2(2)
<= 1/2*W_b)%不在该平面内为真
S_ARB = S_AR*T_inv;
H1 = [(S_ARB(3)*dot_B(1) - S_ARB(1)*dot_B(3))/S_ARB(3),
(S_ARB(3)*dot_B(2) - S_ARB(2)*dot_B(3))/S_ARB(3), 0];
%%%如果在塔挡阴影范围内 eta_sb=0；
%计算阴影遮挡效率 eta_sb
%step1：塔挡损失 定日镜镜面中心是否在阴影中,如果是 eta_sb=0 %i 是-S_i
[P_top_1_intersect,P_top_2_intersect,P_bot_1_intersect,P_bot_2_intersect] =
shadow_range(-S_i, h2) ;%i 是-S_i
isInside = pointInRectangle(O_A(1), O_A(2), P_top_1_intersect, P_top_2_intersect,
P_bot_1_intersect, P_bot_2_intersect);
if isInside %在里面
eta_sb=0;
eta_trunc=0;
Store(num,1)=0;
Store(num,2)=0;
Store(num,3)=eta_cos;
Store(num,4)=0;
else %不在里面
R = 100;
X = linspace(-1/2 * W_a,1/2 * W_a, 5);
Y = linspace(-1/2 * W_b,1/2 *W_b, 5);
% Assuming S_i, S_AR, and T are defined elsewhere
T_inv = inv(T);
% 在指定的半径内寻找定日镜
distances = sqrt((data(1:num-1, 1) - O_A(1)).^2 + (data(1:num-1, 2) - O_A(2)).^2);
valid_indices = find(distances <= R);
% Main loops
valid_point_count = 0;%有效点
valid_points = 0;%每个点的有效光线
obstructed_mirrors = [];
for x = X
for y = Y
dot_A = [x, y, 0];
dot_AG = dot_A * T + O_A;
is_valid = true;
for idx = valid_indices'
O_B = [data(idx, 1), data(idx, 2), h2];
dot_B = (dot_AG - O_B)*T_inv;
S_iB = S_i*T_inv;
H2 = [(S_iB(3)*dot_B(1) - S_iB(1)*dot_B(3))/S_iB(3),
(S_iB(3)*dot_B(2) - S_iB(2)*dot_B(3))/S_iB(3), 0];
if ~(-1/2*W_a <= H2(1) && H2(1) <= 1/2*W_a && -1/2*W_b <= H2(2) && H2(2)
<= 1/2*W_b)%不在该平面内为真
S_ARB = S_AR*T_inv;
H1 = [(S_ARB(3)*dot_B(1) - S_ARB(1)*dot_B(3))/S_ARB(3),
(S_ARB(3)*dot_B(2) - S_ARB(2)*dot_B(3))/S_ARB(3), 0];
if -1/2*W_a <= H1(1) && H1(1) <= 1/2*W_a && -1/2*W_b <= H1(2) &&
H1(2) <= 1/2*W_b %在该平面内为真
obstructed_mirrors = [obstructed_mirrors; O_B,
distances(idx)];
is_valid = false;
break;
end
else%在该平面内为真
obstructed_mirrors = [obstructed_mirrors; O_B, distances(idx)];
is_valid = false;
break;
end
end
if is_valid
valid_point_count = valid_point_count + 1;
%这个点是有效的计算集热器接收到的光
% 计算入射向量和法向量之间的夹角
cosTheta_i = dot(S_i, S_n) / (norm(S_i) * norm(S_n));
S_nt = S_n / cosTheta_i;
d_ki = S_nt - S_AR;
d_ki = d_ki / norm(d_ki);
theta_range = linspace(-4.65e-3,4.65e-3, 5);
theta2_range = linspace(0, 2*pi, 12);
% 循环 theta 值
for theta = theta_range
d_ki2 = tan(theta) * d_ki;
% 循环不同的 theta2 值
for theta2 = theta2_range
RR = rotationMatrix(S_AR, -theta2); %顺时针
d_tki = RR * d_ki2';
d_tki=d_tki';
r = d_tki + S_AR;
% 解出 t 的范围 判断 r 是否与集热器相交 地系坐标
a = r(1)^2 + r(2)^2;
b = 2 * (dot_AG(1) * r(1) + dot_AG(2) * r(2));
c = dot_AG(1)^2 + dot_AG(2)^2 - 49/4;
discriminant = b^2 - 4*a*c;
if discriminant >= 0
t1 = (-b + sqrt(discriminant)) / (2*a);
t2 = (-b - sqrt(discriminant)) / (2*a);
% 使用 t 的范围得到 z 的范围 地系坐标
z1 = dot_AG(3) + t1 * r(3);
z2 = dot_AG(3) + t2 * r(3);
% 检查 z 的范围是否与(h1-4,h1+4]有交集
if -1/2*W_a <= H1(1) && H1(1) <= 1/2*W_a && -1/2*W_b <= H1(2) &&
H1(2) <= 1/2*W_b %在该平面内为真
obstructed_mirrors = [obstructed_mirrors; O_B,
distances(idx)];
is_valid = false;
break;
end
else%在该平面内为真
obstructed_mirrors = [obstructed_mirrors; O_B, distances(idx)];
is_valid = false;
break;
end
end
if is_valid
valid_point_count = valid_point_count + 1;
%这个点是有效的计算集热器接收到的光
% 计算入射向量和法向量之间的夹角
cosTheta_i = dot(S_i, S_n) / (norm(S_i) * norm(S_n));
S_nt = S_n / cosTheta_i;
d_ki = S_nt - S_AR;
d_ki = d_ki / norm(d_ki);
theta_range = linspace(-4.65e-3,4.65e-3, 5);
theta2_range = linspace(0, 2*pi, 12);
% 循环 theta 值
for theta = theta_range
d_ki2 = tan(theta) * d_ki;
% 循环不同的 theta2 值
for theta2 = theta2_range
RR = rotationMatrix(S_AR, -theta2); %顺时针
d_tki = RR * d_ki2';
d_tki=d_tki';
r = d_tki + S_AR;
% 解出 t 的范围 判断 r 是否与集热器相交 地系坐标
a = r(1)^2 + r(2)^2;
b = 2 * (dot_AG(1) * r(1) + dot_AG(2) * r(2));
c = dot_AG(1)^2 + dot_AG(2)^2 - 49/4;
discriminant = b^2 - 4*a*c;
if discriminant >= 0
t1 = (-b + sqrt(discriminant)) / (2*a);
t2 = (-b - sqrt(discriminant)) / (2*a);
% 使用 t 的范围得到 z 的范围 地系坐标
z1 = dot_AG(3) + t1 * r(3);
z2 = dot_AG(3) + t2 * r(3);
% 检查 z 的范围是否与(h1-4,h1+4]有交集
z_min = min(z1, z2);
z_max = max(z1, z2);
if (z_min <= (h1+4) && z_max > (h1-4))%有交集
valid_points = valid_points + 1;
end
end
end
end
end
end
end
eta_sb = valid_point_count / (length(X) * length(Y));
Dot_Sum_light=length(theta_range)*length(theta2_range);%每个区域总光线数量
eta_trunc=valid_points/(valid_point_count*Dot_Sum_light);
eta=eta_cos*eta_sb*eta_at*eta_trunc*eta_ref;%计算光学效率 eta
Store(num,1)=eta;
Store(num,2)=eta_sb;
Store(num,3)=eta_cos;
Store(num,4)=eta_trunc;
end
%计算定日镜场的输出热功率 E_field
Ai=W_a*W_b;%定日镜镜面积
E_field=E_field+Ai*Store(num,1);%E_field=E_field+Ai*eta;
if num== 2836
disp(num);
end
end
E_field=E_field*DNI;
result=[result;mean(Store) E_field/(num*W_a*W_b)];
end
R=mean(result);
EEEE=[W_a W_b R(5)*W_a*W_b*num];
%检验 入射的相反与法向量的夹角 反射与法向量夹角
% % 计算入射向量和法向量之间的夹角
% cosTheta_i = dot(S_i, S_n) / (norm(S_i) * norm(S_n));
% angle_i = acos(cosTheta_i);
%
% % 计算反射向量
% S_r =S_AR;
%
% % 计算反射向量和法向量之间的夹角
% cosTheta_r = dot(S_r, S_n) / (norm(S_r) * norm(S_n));
% angle_r = acos(cosTheta_r);
%
z_min = min(z1, z2);
z_max = max(z1, z2);
if (z_min <= (h1+4) && z_max > (h1-4))%有交集
valid_points = valid_points + 1;
end
end
end
end
end
end
end
eta_sb = valid_point_count / (length(X) * length(Y));
Dot_Sum_light=length(theta_range)*length(theta2_range);%每个区域总光线数量
eta_trunc=valid_points/(valid_point_count*Dot_Sum_light);
eta=eta_cos*eta_sb*eta_at*eta_trunc*eta_ref;%计算光学效率 eta
Store(num,1)=eta;
Store(num,2)=eta_sb;
Store(num,3)=eta_cos;
Store(num,4)=eta_trunc;
end
%计算定日镜场的输出热功率 E_field
Ai=W_a*W_b;%定日镜镜面积
E_field=E_field+Ai*Store(num,1);%E_field=E_field+Ai*eta;
if num== 2836
disp(num);
end
end
E_field=E_field*DNI;
result=[result;mean(Store) E_field/(num*W_a*W_b)];
end
R=mean(result);
EEEE=[W_a W_b R(5)*W_a*W_b*num];
%检验 入射的相反与法向量的夹角 反射与法向量夹角
% % 计算入射向量和法向量之间的夹角
% cosTheta_i = dot(S_i, S_n) / (norm(S_i) * norm(S_n));
% angle_i = acos(cosTheta_i);
%
% % 计算反射向量
% S_r =S_AR;
%
% % 计算反射向量和法向量之间的夹角
% cosTheta_r = dot(S_r, S_n) / (norm(S_r) * norm(S_n));
% angle_r = acos(cosTheta_r);
%
% % 转换夹角为度
% angle_i_deg = rad2deg(angle_i)
% angle_r_deg = rad2deg(angle_r)
end
function [P_top_1_intersect,P_top_2_intersect,P_bot_1_intersect,P_bot_2_intersect] =
shadow_range(i, h2) %i 是-S_i
% 参数:
% i: 入射向量, 例如 [1,1,1]
% h2: 圆柱阴影所在的 z 坐标
% 入射向量规范化
i = i / norm(i);
% 圆柱参数
global h1;
R = 7/2; % 圆柱半径
O_top= [0 0 (h1 + R)]; % 圆柱上底面中心坐标
% 计算与入射光 i 垂直的直径方向
diameter_direction = cross(i, [0, 0, 1]);
diameter_direction = diameter_direction / norm(diameter_direction); % 单位化
% 上底面和下底面上与入射光垂直的直径的两个端点
P_top_1 = O_top + R * diameter_direction;
P_top_2 = O_top - R * diameter_direction;
P_bot_1 = P_top_1 - [0, 0, 2*R];
P_bot_2 = P_top_2 - [0, 0, 2*R];
% 计算在 z=h 时的交点
t_top_1 = (h2 - P_top_1(3)) / i(3);
P_top_1_intersect = P_top_1 + t_top_1 * i;
t_top_2 = (h2 - P_top_2(3)) / i(3);
P_top_2_intersect = P_top_2 + t_top_2 * i;
t_bot_1 = (h2 - P_bot_1(3)) / i(3);
P_bot_1_intersect = P_bot_1 + t_bot_1 * i;
t_bot_2 = (h2 - P_bot_2(3)) / i(3);
P_bot_2_intersect = P_bot_2 + t_bot_2 * i;
end
function isInside = pointInRectangle(x, y, P_top_1_intersect, P_top_2_intersect,
P_bot_1_intersect, P_bot_2_intersect)
% 定义矩形的顶点
polygon = [P_top_1_intersect; P_top_2_intersect; P_bot_2_intersect; P_bot_1_intersect;
P_top_1_intersect];
% 计算交点数
crossings = 0;
for i = 1:length(polygon)-1
if ((polygon(i, 2) > y) ~= (polygon(i+1, 2) > y)) && ...
(x < (polygon(i+1, 1) - polygon(i, 1)) * (y - polygon(i, 2)) / ...
(polygon(i+1, 2) - polygon(i, 2)) + polygon(i, 1))
% % 转换夹角为度
% angle_i_deg = rad2deg(angle_i)
% angle_r_deg = rad2deg(angle_r)
end
function [P_top_1_intersect,P_top_2_intersect,P_bot_1_intersect,P_bot_2_intersect] =
shadow_range(i, h2) %i 是-S_i
% 参数:
% i: 入射向量, 例如 [1,1,1]
% h2: 圆柱阴影所在的 z 坐标
% 入射向量规范化
i = i / norm(i);
% 圆柱参数
global h1;
R = 7/2; % 圆柱半径
O_top= [0 0 (h1 + R)]; % 圆柱上底面中心坐标
% 计算与入射光 i 垂直的直径方向
diameter_direction = cross(i, [0, 0, 1]);
diameter_direction = diameter_direction / norm(diameter_direction); % 单位化
% 上底面和下底面上与入射光垂直的直径的两个端点
P_top_1 = O_top + R * diameter_direction;
P_top_2 = O_top - R * diameter_direction;
P_bot_1 = P_top_1 - [0, 0, 2*R];
P_bot_2 = P_top_2 - [0, 0, 2*R];
% 计算在 z=h 时的交点
t_top_1 = (h2 - P_top_1(3)) / i(3);
P_top_1_intersect = P_top_1 + t_top_1 * i;
t_top_2 = (h2 - P_top_2(3)) / i(3);
P_top_2_intersect = P_top_2 + t_top_2 * i;
t_bot_1 = (h2 - P_bot_1(3)) / i(3);
P_bot_1_intersect = P_bot_1 + t_bot_1 * i;
t_bot_2 = (h2 - P_bot_2(3)) / i(3);
P_bot_2_intersect = P_bot_2 + t_bot_2 * i;
end
function isInside = pointInRectangle(x, y, P_top_1_intersect, P_top_2_intersect,
P_bot_1_intersect, P_bot_2_intersect)
% 定义矩形的顶点
polygon = [P_top_1_intersect; P_top_2_intersect; P_bot_2_intersect; P_bot_1_intersect;
P_top_1_intersect];
% 计算交点数
crossings = 0;
for i = 1:length(polygon)-1
if ((polygon(i, 2) > y) ~= (polygon(i+1, 2) > y)) && ...
(x < (polygon(i+1, 1) - polygon(i, 1)) * (y - polygon(i, 2)) / ...
(polygon(i+1, 2) - polygon(i, 2)) + polygon(i, 1))
crossings = crossings + 1;
end
end
% 根据交点数判断点是否在矩形内
if mod(crossings, 2) == 1
isInside = true;
else
isInside = false;
end
end
function RR = rotationMatrix(axis, angle)
c = cos(angle);
s = sin(angle);
t = 1 - c;
x = axis(1);
y = axis(2);
z = axis(3);
RR = [
t*x*x + c, t*x*y - s*z, t*x*z + s*y;
t*x*y + s*z, t*y*y + c, t*y*z - s*x;
t*x*z - s*y, t*y*z + s*x, t*z*z + c
];
end
function data=data_func()% 这里求出的定日镜坐标是相对坐标
global W_a;
global L0;
R0 = 100;
delta_D=5+W_a+0.05;%给点容差%% !
coordinates = []; % 保存所有坐标
maxNumMirrors = 10000; % 估计的最大定日镜数
coordinates = NaN(maxNumMirrors, 2); % 预先分配空间
idx = 1; % 当前填充到的坐标索引
% 第一区域的五层
for numj = 1:5
R = R0 + (numj-1)*delta_D;
% 使用余弦定理计算 theta_D
cos_theta_D = (2*R^2 - delta_D^2) / (2*R^2);
theta_D = acos(cos_theta_D);
ND = floor(2*pi/theta_D);
% 重新计算 theta_D 的值
theta_D = 2*pi/ND;
% 计算 cos 和 sin 值
cos_values = cos((0:ND-1) * theta_D);
sin_values = sin((0:ND-1) * theta_D);
crossings = crossings + 1;
end
end
% 根据交点数判断点是否在矩形内
if mod(crossings, 2) == 1
isInside = true;
else
isInside = false;
end
end
function RR = rotationMatrix(axis, angle)
c = cos(angle);
s = sin(angle);
t = 1 - c;
x = axis(1);
y = axis(2);
z = axis(3);
RR = [
t*x*x + c, t*x*y - s*z, t*x*z + s*y;
t*x*y + s*z, t*y*y + c, t*y*z - s*x;
t*x*z - s*y, t*y*z + s*x, t*z*z + c
];
end
function data=data_func()% 这里求出的定日镜坐标是相对坐标
global W_a;
global L0;
R0 = 100;
delta_D=5+W_a+0.05;%给点容差%% !
coordinates = []; % 保存所有坐标
maxNumMirrors = 10000; % 估计的最大定日镜数
coordinates = NaN(maxNumMirrors, 2); % 预先分配空间
idx = 1; % 当前填充到的坐标索引
% 第一区域的五层
for numj = 1:5
R = R0 + (numj-1)*delta_D;
% 使用余弦定理计算 theta_D
cos_theta_D = (2*R^2 - delta_D^2) / (2*R^2);
theta_D = acos(cos_theta_D);
ND = floor(2*pi/theta_D);
% 重新计算 theta_D 的值
theta_D = 2*pi/ND;
% 计算 cos 和 sin 值
cos_values = cos((0:ND-1) * theta_D);
sin_values = sin((0:ND-1) * theta_D);
% 计算所有定日镜的坐标
x_values = R * cos_values;
y_values = R * sin_values;
% 更新坐标数组
coordinates(idx:idx+ND-1, :) = [x_values', y_values'];
idx = idx + ND;
end
i = 2; % 从第二个区域开始
while true
R = R + delta_D;% 对于新区域的第一层，从上一个区域的最后一层的 R 增加一个 delta_D
% 检查 R 的条件
if R > (350+L0)
break;
end
% 使用余弦定理计算 theta_D
cos_theta_D = (2*R^2 - delta_D^2) / (2*R^2);
theta_D = acos(cos_theta_D);
ND = floor(2*pi/theta_D);
theta_D = 2*pi/ND; % 重新计算 theta_D 的值
% 计算第一层的坐标
cos_values = cos((0:ND-1) * theta_D);
sin_values = sin((0:ND-1) * theta_D);
x_values = R * cos_values;
y_values = R * sin_values;
coordinates(idx:idx+ND-1, :) = [x_values', y_values'];
idx = idx + ND;
% 保存第一层坐标的距离
distance_first_layer = sqrt((x_values(2) - x_values(1))^2 + (y_values(2) -
y_values(1))^2);
% 其他层的坐标
layer = 2;
while true
prev=sqrt((coordinates(idx-1,1) - coordinates(idx-2,1))^2 + (coordinates(idx-1,2)-
coordinates(idx-2,2))^2);
R = R + sqrt(delta_D^2 - (prev/2)^2);
% 检查 R 的条件
if R > (350+L0)%%%%
break;
end
theta_D0 = atan2(coordinates(idx-ND,2), coordinates(idx-ND,1)) + theta_D/2;
cos_values = cos((0:ND-1) * theta_D + theta_D0); % 角平分线
sin_values = sin((0:ND-1) * theta_D + theta_D0); % 角平分线
x_values = R * cos_values;
y_values = R * sin_values;
% 计算所有定日镜的坐标
x_values = R * cos_values;
y_values = R * sin_values;
% 更新坐标数组
coordinates(idx:idx+ND-1, :) = [x_values', y_values'];
idx = idx + ND;
end
i = 2; % 从第二个区域开始
while true
R = R + delta_D;% 对于新区域的第一层，从上一个区域的最后一层的 R 增加一个 delta_D
% 检查 R 的条件
if R > (350+L0)
break;
end
% 使用余弦定理计算 theta_D
cos_theta_D = (2*R^2 - delta_D^2) / (2*R^2);
theta_D = acos(cos_theta_D);
ND = floor(2*pi/theta_D);
theta_D = 2*pi/ND; % 重新计算 theta_D 的值
% 计算第一层的坐标
cos_values = cos((0:ND-1) * theta_D);
sin_values = sin((0:ND-1) * theta_D);
x_values = R * cos_values;
y_values = R * sin_values;
coordinates(idx:idx+ND-1, :) = [x_values', y_values'];
idx = idx + ND;
% 保存第一层坐标的距离
distance_first_layer = sqrt((x_values(2) - x_values(1))^2 + (y_values(2) -
y_values(1))^2);
% 其他层的坐标
layer = 2;
while true
prev=sqrt((coordinates(idx-1,1) - coordinates(idx-2,1))^2 + (coordinates(idx-1,2)-
coordinates(idx-2,2))^2);
R = R + sqrt(delta_D^2 - (prev/2)^2);
% 检查 R 的条件
if R > (350+L0)%%%%
break;
end
theta_D0 = atan2(coordinates(idx-ND,2), coordinates(idx-ND,1)) + theta_D/2;
cos_values = cos((0:ND-1) * theta_D + theta_D0); % 角平分线
sin_values = sin((0:ND-1) * theta_D + theta_D0); % 角平分线
x_values = R * cos_values;
y_values = R * sin_values;
coordinates(idx:idx+ND-1, :) = [x_values', y_values'];
idx = idx + ND;
% 检查距离条件
distance = sqrt((x_values(2) - x_values(1))^2 + (y_values(2) - y_values(1))^2);
if distance > 1.33 * distance_first_layer
break;
end
layer = layer + 1;
end
i = i + 1;
end
% 删除未使用的空间
coordinates = coordinates(~isnan(coordinates(:,1)), :);
%筛选定日镜坐标
% 检查约束条件
N = size(coordinates,1);
valid_indices = [];
for i = 1:N
x1 = coordinates(i,1);
y1 = coordinates(i,2);
if x1^2 + (y1-L0)^2 <= 350^2
valid_indices = [valid_indices; i];
end
end
data=[coordinates(valid_indices,1), coordinates(valid_indices,2)];
error2 = check_distance(data);
if ~error2
error('定日镜距离不符合');
end
data(:,3) = sqrt(data(:,1).^2 + data(:,2).^2);
data = sortrows(data, 3);
end
function result = check_distance(coordinates)% 检验点是否满足两个定日镜之间的距离要大于等于 delta_D
% 给定的参数
global W_a;
delta_D = 5 + W_a;
% 计算所有点之间的距离
distances = pdist(coordinates);
% 找到最小距离
min_distance = min(distances);
% 检查最小距离是否大于等于 delta_D
if min_distance >= delta_D
result = true;
else
coordinates(idx:idx+ND-1, :) = [x_values', y_values'];
idx = idx + ND;
% 检查距离条件
distance = sqrt((x_values(2) - x_values(1))^2 + (y_values(2) - y_values(1))^2);
if distance > 1.33 * distance_first_layer
break;
end
layer = layer + 1;
end
i = i + 1;
end
% 删除未使用的空间
coordinates = coordinates(~isnan(coordinates(:,1)), :);
%筛选定日镜坐标
% 检查约束条件
N = size(coordinates,1);
valid_indices = [];
for i = 1:N
x1 = coordinates(i,1);
y1 = coordinates(i,2);
if x1^2 + (y1-L0)^2 <= 350^2
valid_indices = [valid_indices; i];
end
end
data=[coordinates(valid_indices,1), coordinates(valid_indices,2)];
error2 = check_distance(data);
if ~error2
error('定日镜距离不符合');
end
data(:,3) = sqrt(data(:,1).^2 + data(:,2).^2);
data = sortrows(data, 3);
end
function result = check_distance(coordinates)% 检验点是否满足两个定日镜之间的距离要大于等于 delta_D
% 给定的参数
global W_a;
delta_D = 5 + W_a;
% 计算所有点之间的距离
distances = pdist(coordinates);
% 找到最小距离
min_distance = min(distances);
% 检查最小距离是否大于等于 delta_D
if min_distance >= delta_D
result = true;
else
result = false;
end
end程序三问题三代码（ThiQues.m）
global W_a;
global W_b;
global h2;
global h1;
global L0;
WAAA=5.5:0.1:6.5;
WBBB=5.5:0.1:6.5;
R_year=[];
E_year=[];
for W_a=6.4
for W_b=6.4
if W_a<W_b
continue;
else
h2_range=linspace(W_b/2,6, 5);
for h2=3.2
h1=80;
L0_range=linspace(78,85, 9);
for L0=79
%L0=0;
%% 入射角
%% 计算太阳赤纬角 delta
% 输入：D 为以春分作为第 0 天起算的天数
%D = 31; % 4 月 21 日对应的天数
%%输入 D 得到 R
% R=Main(D);
D=[-59 -28 0 31 61 92 122 153 184 214 245 275 ];
R_total=[];
EEEE_total=[];
for i=1:1:size(D,1)
disp(i);
[R,EEEE]=Main(D(i));
R_total=[R_total;R];
EEEE_total=[EEEE_total;EEEE];
end
R_year=[R_year;mean(R_total)];
E_year=[E_year;mean(EEEE_total)];
result = false;
end
end
程序三 问题三代码（ThiQues.m）
global W_a;
global W_b;
global h2;
global h1;
global L0;
WAAA=5.5:0.1:6.5;
WBBB=5.5:0.1:6.5;
R_year=[];
E_year=[];
for W_a=6.4
for W_b=6.4
if W_a<W_b
continue;
else
h2_range=linspace(W_b/2,6, 5);
for h2=3.2
h1=80;
L0_range=linspace(78,85, 9);
for L0=79
%L0=0;
%% 入射角
%% 计算太阳赤纬角 delta
% 输入：D 为以春分作为第 0 天起算的天数
%D = 31; % 4 月 21 日对应的天数
%%输入 D 得到 R
% R=Main(D);
D=[-59 -28 0 31 61 92 122 153 184 214 245 275 ];
R_total=[];
EEEE_total=[];
for i=1:1:size(D,1)
disp(i);
[R,EEEE]=Main(D(i));
R_total=[R_total;R];
EEEE_total=[EEEE_total;EEEE];
end
R_year=[R_year;mean(R_total)];
E_year=[E_year;mean(EEEE_total)];
end
end
end
end
end
function [R,EEEE]=Main(Day)
global W_a;
global W_b;
global h2;
global h1;
global L0;
EEEE=[];
D=Day;
STT=[9 10.5 12 13.5 15]; %当地时间,假设为 15 点
result=[];%日结果详细矩阵 5*4 5 个时间点 4 个参数 R 才是最终值
eta_ref=0.92;%镜面反射率
O_R=[0 0 h1];%集热器中心坐标
%定日镜选取根据与中心的距离进行排序
% 计算 data 的程序
data=data_func(); %经过排序后的定日镜坐标 data---1745*3 x,y,r
% 计算太阳赤纬角 delta
sin_delta = sin(2*pi*D/365) * sin(2*pi/360 * 23.45);
delta = asin(sin_delta); % 得到赤纬角的弧度值
phi = deg2rad(39.4); % 北纬 39.4 度，转换为弧度
for XH_ST=1:size(STT,2)
ST=STT(XH_ST);
omega =pi/12*(ST-12); % 声明 omega 的值
% 计算 sin(alpha_s)
sin_alpha_s = cos(delta) * cos(phi) * cos(omega) + sin(delta) * sin(phi);
G_0=1.366; %太阳常数
Hh=3;%3km
ha=0.4237 - 0.00821*(6 - Hh)^2;%大气压
hb=0.5055 + 0.00595*(6.5 - Hh)^2;%大气压
hc=0.2711 + 0.01858*(2.5 - Hh)^2;%大气压
DNI = G_0*(ha + hb*exp(-hc/sin_alpha_s)); %直接法辐照度
alpha_s = asin(sin_alpha_s); % 得到 alpha_s 的弧度值
% 根据公式计算 cos(gamma_s)
cos_gamma_s = (sin(delta) - sin_alpha_s * sin(phi)) / (cos(alpha_s) * cos(phi));
if cos_gamma_s > 1
cos_gamma_s = 1;
elseif cos_gamma_s < -1
cos_gamma_s = -1;
end
end
end
end
end
end
function [R,EEEE]=Main(Day)
global W_a;
global W_b;
global h2;
global h1;
global L0;
EEEE=[];
D=Day;
STT=[9 10.5 12 13.5 15]; %当地时间,假设为 15 点
result=[];%日结果详细矩阵 5*4 5 个时间点 4 个参数 R 才是最终值
eta_ref=0.92;%镜面反射率
O_R=[0 0 h1];%集热器中心坐标
%定日镜选取根据与中心的距离进行排序
% 计算 data 的程序
data=data_func(); %经过排序后的定日镜坐标 data---1745*3 x,y,r
% 计算太阳赤纬角 delta
sin_delta = sin(2*pi*D/365) * sin(2*pi/360 * 23.45);
delta = asin(sin_delta); % 得到赤纬角的弧度值
phi = deg2rad(39.4); % 北纬 39.4 度，转换为弧度
for XH_ST=1:size(STT,2)
ST=STT(XH_ST);
omega =pi/12*(ST-12); % 声明 omega 的值
% 计算 sin(alpha_s)
sin_alpha_s = cos(delta) * cos(phi) * cos(omega) + sin(delta) * sin(phi);
G_0=1.366; %太阳常数
Hh=3;%3km
ha=0.4237 - 0.00821*(6 - Hh)^2;%大气压
hb=0.5055 + 0.00595*(6.5 - Hh)^2;%大气压
hc=0.2711 + 0.01858*(2.5 - Hh)^2;%大气压
DNI = G_0*(ha + hb*exp(-hc/sin_alpha_s)); %直接法辐照度
alpha_s = asin(sin_alpha_s); % 得到 alpha_s 的弧度值
% 根据公式计算 cos(gamma_s)
cos_gamma_s = (sin(delta) - sin_alpha_s * sin(phi)) / (cos(alpha_s) * cos(phi));
if cos_gamma_s > 1
cos_gamma_s = 1;
elseif cos_gamma_s < -1
cos_gamma_s = -1;
end
gamma_s = acos(cos_gamma_s); % 得到 gamma_s 的弧度值
x=cos(alpha_s)*sin(gamma_s);
y=cos(alpha_s)*cos(gamma_s);
z=sin(alpha_s);
S_i=[x y z];%太阳入射方向单位向量的相反
Store=zeros(size(data,1),4);%光学效率 阴影遮挡效率 余弦效率 截断效率
E_field=0;%定日镜场的输出热功率
for num=1:1:size(data,1)%%定日镜编号
disp(num);
h2=data(num,4);%定日镜高度
O_A=[data(num,1) data(num,2) h2];%定日镜中心坐标
%计算集热器截断效率 eta_trunc
d_HR=sqrt((O_A(1)-O_R(1))^2+(O_A(2)-O_R(2))^2+(O_A(3)-O_R(3))^2);%定日镜镜面中心与集
热器中心的距离
eta_at = 0.99321 - 0.0001176*d_HR + 1.97e-8 * (d_HR^2); % 计算吸收塔的吸收率
S_AR=(O_R-O_A)/norm(O_R-O_A);%定日镜中心到集热器中心的单位向量 %% 反射角
%% 计算法向量
S_n=(S_i+S_AR)/norm(S_i+S_AR);%法向量单位向量
% 计算方位角 Azimuth (θ)
A_H = atan2(S_n(1), S_n(2));
% 计算俯仰角 Elevation (ϕ)
E_H = acos(S_n(3));
%
T=[-sin(E_H),-sin(A_H)*cos(E_H),cos(A_H)*cos(E_H);cos(E_H),-sin(A_H)*sin(E_H),cos(A_H)*sin(E_H
);0,cos(A_H),sin(A_H)];
M1 = [1, 0, 0; 0, cos(E_H), sin(E_H); 0, -sin(E_H), cos(E_H)];
M2 = [cos(pi-A_H), sin(pi-A_H), 0; -sin(pi-A_H), cos(pi-A_H), 0; 0, 0, 1];
% 计算矩阵乘积
T = M1 * M2; %镜坐标系变为地坐标系的转换矩阵
S_An=[0 0 1];%镜面 A 坐标系上的的法向量
SSS=S_An*T;
tolerance = 1e-10;
difference = abs(SSS - S_n);
if all(difference(:) >tolerance)
error('法向量转换错误');
end
%%
%计算余弦效率 eta_cos
eta_cos = dot(S_i, S_n) / (norm(S_i) * norm(S_n)); % 计算入射向量相反和法向量之间的夹
角的余弦值
%计算阴影遮挡效率 eta_sb
%step1：塔挡损失 定日镜镜面中心是否在阴影中,如果是 eta_sb=0 %i 是-S_i
[P_top_1_intersect,P_top_2_intersect,P_bot_1_intersect,P_bot_2_intersect] =
shadow_range(-S_i, h2) ;%i 是-S_i
gamma_s = acos(cos_gamma_s); % 得到 gamma_s 的弧度值
x=cos(alpha_s)*sin(gamma_s);
y=cos(alpha_s)*cos(gamma_s);
z=sin(alpha_s);
S_i=[x y z];%太阳入射方向单位向量的相反
Store=zeros(size(data,1),4);%光学效率 阴影遮挡效率 余弦效率 截断效率
E_field=0;%定日镜场的输出热功率
for num=1:1:size(data,1)%%定日镜编号
disp(num);
h2=data(num,4);%定日镜高度
O_A=[data(num,1) data(num,2) h2];%定日镜中心坐标
%计算集热器截断效率 eta_trunc
d_HR=sqrt((O_A(1)-O_R(1))^2+(O_A(2)-O_R(2))^2+(O_A(3)-O_R(3))^2);%定日镜镜面中心与集
热器中心的距离
eta_at = 0.99321 - 0.0001176*d_HR + 1.97e-8 * (d_HR^2); % 计算吸收塔的吸收率
S_AR=(O_R-O_A)/norm(O_R-O_A);%定日镜中心到集热器中心的单位向量 %% 反射角
%% 计算法向量
S_n=(S_i+S_AR)/norm(S_i+S_AR);%法向量单位向量
% 计算方位角 Azimuth (θ)
A_H = atan2(S_n(1), S_n(2));
% 计算俯仰角 Elevation (ϕ)
E_H = acos(S_n(3));
%
T=[-sin(E_H),-sin(A_H)*cos(E_H),cos(A_H)*cos(E_H);cos(E_H),-sin(A_H)*sin(E_H),cos(A_H)*sin(E_H
);0,cos(A_H),sin(A_H)];
M1 = [1, 0, 0; 0, cos(E_H), sin(E_H); 0, -sin(E_H), cos(E_H)];
M2 = [cos(pi-A_H), sin(pi-A_H), 0; -sin(pi-A_H), cos(pi-A_H), 0; 0, 0, 1];
% 计算矩阵乘积
T = M1 * M2; %镜坐标系变为地坐标系的转换矩阵
S_An=[0 0 1];%镜面 A 坐标系上的的法向量
SSS=S_An*T;
tolerance = 1e-10;
difference = abs(SSS - S_n);
if all(difference(:) >tolerance)
error('法向量转换错误');
end
%%
%计算余弦效率 eta_cos
eta_cos = dot(S_i, S_n) / (norm(S_i) * norm(S_n)); % 计算入射向量相反和法向量之间的夹
角的余弦值
%计算阴影遮挡效率 eta_sb
%step1：塔挡损失 定日镜镜面中心是否在阴影中,如果是 eta_sb=0 %i 是-S_i
[P_top_1_intersect,P_top_2_intersect,P_bot_1_intersect,P_bot_2_intersect] =
shadow_range(-S_i, h2) ;%i 是-S_i
isInside = pointInRectangle(O_A(1), O_A(2), P_top_1_intersect, P_top_2_intersect,
P_bot_1_intersect, P_bot_2_intersect);
if isInside %在里面
eta_sb=0;
eta_trunc=0;
Store(num,1)=0;
Store(num,2)=0;
Store(num,3)=eta_cos;
Store(num,4)=0;
else %不在里面
R = 100;
X = linspace(-1/2 * W_a,1/2 * W_a, 20);
Y = linspace(-1/2 * W_b,1/2 *W_b, 20);
% Assuming S_i, S_AR, and T are defined elsewhere
T_inv = inv(T);
% 在指定的半径内寻找定日镜
distances = sqrt((data(1:num-1, 1) - O_A(1)).^2 + (data(1:num-1, 2) - O_A(2)).^2);
valid_indices = find(distances <= R);
% Main loops
valid_point_count = 0;%有效点
valid_points = 0;%每个点的有效光线
obstructed_mirrors = [];
for x = X
for y = Y
dot_A = [x, y, 0];
dot_AG = dot_A * T + O_A;
is_valid = true;
for idx = valid_indices'
O_B = [data(idx, 1), data(idx, 2), h2];
dot_B = (dot_AG - O_B)*T_inv;
S_iB = S_i*T_inv;
H2 = [(S_iB(3)*dot_B(1) - S_iB(1)*dot_B(3))/S_iB(3),
(S_iB(3)*dot_B(2) - S_iB(2)*dot_B(3))/S_iB(3), 0];
if ~(-1/2*W_a <= H2(1) && H2(1) <= 1/2*W_a && -1/2*W_b <= H2(2) && H2(2)
<= 1/2*W_b)%不在该平面内为真
S_ARB = S_AR*T_inv;
H1 = [(S_ARB(3)*dot_B(1) - S_ARB(1)*dot_B(3))/S_ARB(3),
(S_ARB(3)*dot_B(2) - S_ARB(2)*dot_B(3))/S_ARB(3), 0];
if -1/2*W_a <= H1(1) && H1(1) <= 1/2*W_a && -1/2*W_b <= H1(2) &&
H1(2) <= 1/2*W_b %在该平面内为真
obstructed_mirrors = [obstructed_mirrors; O_B,
distances(idx)];
is_valid = false;
break;
isInside = pointInRectangle(O_A(1), O_A(2), P_top_1_intersect, P_top_2_intersect,
P_bot_1_intersect, P_bot_2_intersect);
if isInside %在里面
eta_sb=0;
eta_trunc=0;
Store(num,1)=0;
Store(num,2)=0;
Store(num,3)=eta_cos;
Store(num,4)=0;
else %不在里面
R = 100;
X = linspace(-1/2 * W_a,1/2 * W_a, 20);
Y = linspace(-1/2 * W_b,1/2 *W_b, 20);
% Assuming S_i, S_AR, and T are defined elsewhere
T_inv = inv(T);
% 在指定的半径内寻找定日镜
distances = sqrt((data(1:num-1, 1) - O_A(1)).^2 + (data(1:num-1, 2) - O_A(2)).^2);
valid_indices = find(distances <= R);
% Main loops
valid_point_count = 0;%有效点
valid_points = 0;%每个点的有效光线
obstructed_mirrors = [];
for x = X
for y = Y
dot_A = [x, y, 0];
dot_AG = dot_A * T + O_A;
is_valid = true;
for idx = valid_indices'
O_B = [data(idx, 1), data(idx, 2), h2];
dot_B = (dot_AG - O_B)*T_inv;
S_iB = S_i*T_inv;
H2 = [(S_iB(3)*dot_B(1) - S_iB(1)*dot_B(3))/S_iB(3),
(S_iB(3)*dot_B(2) - S_iB(2)*dot_B(3))/S_iB(3), 0];
if ~(-1/2*W_a <= H2(1) && H2(1) <= 1/2*W_a && -1/2*W_b <= H2(2) && H2(2)
<= 1/2*W_b)%不在该平面内为真
S_ARB = S_AR*T_inv;
H1 = [(S_ARB(3)*dot_B(1) - S_ARB(1)*dot_B(3))/S_ARB(3),
(S_ARB(3)*dot_B(2) - S_ARB(2)*dot_B(3))/S_ARB(3), 0];
if -1/2*W_a <= H1(1) && H1(1) <= 1/2*W_a && -1/2*W_b <= H1(2) &&
H1(2) <= 1/2*W_b %在该平面内为真
obstructed_mirrors = [obstructed_mirrors; O_B,
distances(idx)];
is_valid = false;
break;
end
else%在该平面内为真
obstructed_mirrors = [obstructed_mirrors; O_B, distances(idx)];
is_valid = false;
break;
end
end
if is_valid
valid_point_count = valid_point_count + 1;
%这个点是有效的计算集热器接收到的光
% 计算入射向量和法向量之间的夹角
cosTheta_i = dot(S_i, S_n) / (norm(S_i) * norm(S_n));
S_nt = S_n / cosTheta_i;
d_ki = S_nt - S_AR;
d_ki = d_ki / norm(d_ki);
theta_range = linspace(-4.65e-3,4.65e-3, 5);
theta2_range = linspace(0, 2*pi, 12);
% 循环 theta 值
for theta = theta_range
d_ki2 = tan(theta) * d_ki;
% 循环不同的 theta2 值
for theta2 = theta2_range
RR = rotationMatrix(S_AR, -theta2); %顺时针
d_tki = RR * d_ki2';
d_tki=d_tki';
r = d_tki + S_AR;
% 解出 t 的范围 判断 r 是否与集热器相交 地系坐标
a = r(1)^2 + r(2)^2;
b = 2 * (dot_AG(1) * r(1) + dot_AG(2) * r(2));
c = dot_AG(1)^2 + dot_AG(2)^2 - 49/4;
discriminant = b^2 - 4*a*c;
if discriminant >= 0
t1 = (-b + sqrt(discriminant)) / (2*a);
t2 = (-b - sqrt(discriminant)) / (2*a);
% 使用 t 的范围得到 z 的范围 地系坐标
z1 = dot_AG(3) + t1 * r(3);
z2 = dot_AG(3) + t2 * r(3);
% 检查 z 的范围是否与(h1-4,h1+4]有交集
z_min = min(z1, z2);
z_max = max(z1, z2);
if (z_min <= (h1+4) && z_max > (h1-4))%有交集
valid_points = valid_points + 1;
end
end
end
else%在该平面内为真
obstructed_mirrors = [obstructed_mirrors; O_B, distances(idx)];
is_valid = false;
break;
end
end
if is_valid
valid_point_count = valid_point_count + 1;
%这个点是有效的计算集热器接收到的光
% 计算入射向量和法向量之间的夹角
cosTheta_i = dot(S_i, S_n) / (norm(S_i) * norm(S_n));
S_nt = S_n / cosTheta_i;
d_ki = S_nt - S_AR;
d_ki = d_ki / norm(d_ki);
theta_range = linspace(-4.65e-3,4.65e-3, 5);
theta2_range = linspace(0, 2*pi, 12);
% 循环 theta 值
for theta = theta_range
d_ki2 = tan(theta) * d_ki;
% 循环不同的 theta2 值
for theta2 = theta2_range
RR = rotationMatrix(S_AR, -theta2); %顺时针
d_tki = RR * d_ki2';
d_tki=d_tki';
r = d_tki + S_AR;
% 解出 t 的范围 判断 r 是否与集热器相交 地系坐标
a = r(1)^2 + r(2)^2;
b = 2 * (dot_AG(1) * r(1) + dot_AG(2) * r(2));
c = dot_AG(1)^2 + dot_AG(2)^2 - 49/4;
discriminant = b^2 - 4*a*c;
if discriminant >= 0
t1 = (-b + sqrt(discriminant)) / (2*a);
t2 = (-b - sqrt(discriminant)) / (2*a);
% 使用 t 的范围得到 z 的范围 地系坐标
z1 = dot_AG(3) + t1 * r(3);
z2 = dot_AG(3) + t2 * r(3);
% 检查 z 的范围是否与(h1-4,h1+4]有交集
z_min = min(z1, z2);
z_max = max(z1, z2);
if (z_min <= (h1+4) && z_max > (h1-4))%有交集
valid_points = valid_points + 1;
end
end
end
end
end
end
end
eta_sb = valid_point_count / (length(X) * length(Y));
Dot_Sum_light=length(theta_range)*length(theta2_range);%每个区域总光线数量
eta_trunc=valid_points/(valid_point_count*Dot_Sum_light);
eta=eta_cos*eta_sb*eta_at*eta_trunc*eta_ref;%计算光学效率 eta
Store(num,1)=eta;
Store(num,2)=eta_sb;
Store(num,3)=eta_cos;
Store(num,4)=eta_trunc;
end
%计算定日镜场的输出热功率 E_field
Ai=W_a*W_b;%定日镜镜面积
E_field=E_field+Ai*Store(num,1);%E_field=E_field+Ai*eta;
end
E_field=E_field*DNI;
result=[result;mean(Store) E_field/(num*W_a*W_b)];
end
R=mean(result);
EEEE=[W_a W_b R(5)*W_a*W_b*num];
end
function [P_top_1_intersect,P_top_2_intersect,P_bot_1_intersect,P_bot_2_intersect] =
shadow_range(i, h2) %i 是-S_i
% 参数:
% i: 入射向量, 例如 [1,1,1]
% h2: 圆柱阴影所在的 z 坐标
% 入射向量规范化
i = i / norm(i);
% 圆柱参数
global h1;
R = 7/2; % 圆柱半径
O_top= [0 0 (h1 + R)]; % 圆柱上底面中心坐标
% 计算与入射光 i 垂直的直径方向
diameter_direction = cross(i, [0, 0, 1]);
diameter_direction = diameter_direction / norm(diameter_direction); % 单位化
% 上底面和下底面上与入射光垂直的直径的两个端点
P_top_1 = O_top + R * diameter_direction;
P_top_2 = O_top - R * diameter_direction;
P_bot_1 = P_top_1 - [0, 0, 2*R];
end
end
end
end
end
eta_sb = valid_point_count / (length(X) * length(Y));
Dot_Sum_light=length(theta_range)*length(theta2_range);%每个区域总光线数量
eta_trunc=valid_points/(valid_point_count*Dot_Sum_light);
eta=eta_cos*eta_sb*eta_at*eta_trunc*eta_ref;%计算光学效率 eta
Store(num,1)=eta;
Store(num,2)=eta_sb;
Store(num,3)=eta_cos;
Store(num,4)=eta_trunc;
end
%计算定日镜场的输出热功率 E_field
Ai=W_a*W_b;%定日镜镜面积
E_field=E_field+Ai*Store(num,1);%E_field=E_field+Ai*eta;
end
E_field=E_field*DNI;
result=[result;mean(Store) E_field/(num*W_a*W_b)];
end
R=mean(result);
EEEE=[W_a W_b R(5)*W_a*W_b*num];
end
function [P_top_1_intersect,P_top_2_intersect,P_bot_1_intersect,P_bot_2_intersect] =
shadow_range(i, h2) %i 是-S_i
% 参数:
% i: 入射向量, 例如 [1,1,1]
% h2: 圆柱阴影所在的 z 坐标
% 入射向量规范化
i = i / norm(i);
% 圆柱参数
global h1;
R = 7/2; % 圆柱半径
O_top= [0 0 (h1 + R)]; % 圆柱上底面中心坐标
% 计算与入射光 i 垂直的直径方向
diameter_direction = cross(i, [0, 0, 1]);
diameter_direction = diameter_direction / norm(diameter_direction); % 单位化
% 上底面和下底面上与入射光垂直的直径的两个端点
P_top_1 = O_top + R * diameter_direction;
P_top_2 = O_top - R * diameter_direction;
P_bot_1 = P_top_1 - [0, 0, 2*R];
P_bot_2 = P_top_2 - [0, 0, 2*R];
% 计算在 z=h 时的交点
t_top_1 = (h2 - P_top_1(3)) / i(3);
P_top_1_intersect = P_top_1 + t_top_1 * i;
t_top_2 = (h2 - P_top_2(3)) / i(3);
P_top_2_intersect = P_top_2 + t_top_2 * i;
t_bot_1 = (h2 - P_bot_1(3)) / i(3);
P_bot_1_intersect = P_bot_1 + t_bot_1 * i;
t_bot_2 = (h2 - P_bot_2(3)) / i(3);
P_bot_2_intersect = P_bot_2 + t_bot_2 * i;
end
function isInside = pointInRectangle(x, y, P_top_1_intersect, P_top_2_intersect,
P_bot_1_intersect, P_bot_2_intersect)
% 定义矩形的顶点
polygon = [P_top_1_intersect; P_top_2_intersect; P_bot_2_intersect; P_bot_1_intersect;
P_top_1_intersect];
% 计算交点数
crossings = 0;
for i = 1:length(polygon)-1
if ((polygon(i, 2) > y) ~= (polygon(i+1, 2) > y)) && ...
(x < (polygon(i+1, 1) - polygon(i, 1)) * (y - polygon(i, 2)) / ...
(polygon(i+1, 2) - polygon(i, 2)) + polygon(i, 1))
crossings = crossings + 1;
end
end
% 根据交点数判断点是否在矩形内
if mod(crossings, 2) == 1
isInside = true;
else
isInside = false;
end
end
function RR = rotationMatrix(axis, angle)
c = cos(angle);
s = sin(angle);
t = 1 - c;
x = axis(1);
y = axis(2);
z = axis(3);
RR = [
t*x*x + c, t*x*y - s*z, t*x*z + s*y;
t*x*y + s*z, t*y*y + c, t*y*z - s*x;
t*x*z - s*y, t*y*z + s*x, t*z*z + c
];
P_bot_2 = P_top_2 - [0, 0, 2*R];
% 计算在 z=h 时的交点
t_top_1 = (h2 - P_top_1(3)) / i(3);
P_top_1_intersect = P_top_1 + t_top_1 * i;
t_top_2 = (h2 - P_top_2(3)) / i(3);
P_top_2_intersect = P_top_2 + t_top_2 * i;
t_bot_1 = (h2 - P_bot_1(3)) / i(3);
P_bot_1_intersect = P_bot_1 + t_bot_1 * i;
t_bot_2 = (h2 - P_bot_2(3)) / i(3);
P_bot_2_intersect = P_bot_2 + t_bot_2 * i;
end
function isInside = pointInRectangle(x, y, P_top_1_intersect, P_top_2_intersect,
P_bot_1_intersect, P_bot_2_intersect)
% 定义矩形的顶点
polygon = [P_top_1_intersect; P_top_2_intersect; P_bot_2_intersect; P_bot_1_intersect;
P_top_1_intersect];
% 计算交点数
crossings = 0;
for i = 1:length(polygon)-1
if ((polygon(i, 2) > y) ~= (polygon(i+1, 2) > y)) && ...
(x < (polygon(i+1, 1) - polygon(i, 1)) * (y - polygon(i, 2)) / ...
(polygon(i+1, 2) - polygon(i, 2)) + polygon(i, 1))
crossings = crossings + 1;
end
end
% 根据交点数判断点是否在矩形内
if mod(crossings, 2) == 1
isInside = true;
else
isInside = false;
end
end
function RR = rotationMatrix(axis, angle)
c = cos(angle);
s = sin(angle);
t = 1 - c;
x = axis(1);
y = axis(2);
z = axis(3);
RR = [
t*x*x + c, t*x*y - s*z, t*x*z + s*y;
t*x*y + s*z, t*y*y + c, t*y*z - s*x;
t*x*z - s*y, t*y*z + s*x, t*z*z + c
];
end
function data=data_func()% 这里求出的定日镜坐标是相对坐标
global W_a;
global L0;
global h2;
h0=h2;
R0 = 100;
delta_D=5+W_a+0.05;%给点容差
coordinates = []; % 保存所有坐标
maxNumMirrors = 10000; % 估计的最大定日镜数
coordinates = NaN(maxNumMirrors, 3); % 预先分配空间
idx = 1; % 当前填充到的坐标索引
for numj = 1
R = R0 + (numj-1)*delta_D;
% 使用余弦定理计算 theta_D
cos_theta_D = (2*R^2 - delta_D^2) / (2*R^2);
theta_D = acos(cos_theta_D);
ND = floor(2*pi/theta_D);
% 重新计算 theta_D 的值
theta_D = 2*pi/ND;
% 计算 cos 和 sin 值
cos_values = cos((0:ND-1) * theta_D);
sin_values = sin((0:ND-1) * theta_D);
% 计算所有定日镜的坐标
x_values = R * cos_values;
y_values = R * sin_values;
h_values=h0*ones(ND,1);
% 更新坐标数组
coordinates(idx:idx+ND-1, 1:2) = [x_values', y_values'];
coordinates(idx:idx+ND-1, 3) = h_values;
idx = idx + ND;
end
i = 2; % 从第二个区域开始
while true
R = R + delta_D;% 对于新区域的第一层，从上一个区域的最后一层的 R 增加一个 delta_D
h=h0+(i-1)*0.5;
% 检查 R 的条件
if R > (350+L0)
break;
end
% 使用余弦定理计算 theta_D
cos_theta_D = (2*R^2 - delta_D^2) / (2*R^2);
theta_D = acos(cos_theta_D);
ND = floor(2*pi/theta_D);
end
function data=data_func()% 这里求出的定日镜坐标是相对坐标
global W_a;
global L0;
global h2;
h0=h2;
R0 = 100;
delta_D=5+W_a+0.05;%给点容差
coordinates = []; % 保存所有坐标
maxNumMirrors = 10000; % 估计的最大定日镜数
coordinates = NaN(maxNumMirrors, 3); % 预先分配空间
idx = 1; % 当前填充到的坐标索引
for numj = 1
R = R0 + (numj-1)*delta_D;
% 使用余弦定理计算 theta_D
cos_theta_D = (2*R^2 - delta_D^2) / (2*R^2);
theta_D = acos(cos_theta_D);
ND = floor(2*pi/theta_D);
% 重新计算 theta_D 的值
theta_D = 2*pi/ND;
% 计算 cos 和 sin 值
cos_values = cos((0:ND-1) * theta_D);
sin_values = sin((0:ND-1) * theta_D);
% 计算所有定日镜的坐标
x_values = R * cos_values;
y_values = R * sin_values;
h_values=h0*ones(ND,1);
% 更新坐标数组
coordinates(idx:idx+ND-1, 1:2) = [x_values', y_values'];
coordinates(idx:idx+ND-1, 3) = h_values;
idx = idx + ND;
end
i = 2; % 从第二个区域开始
while true
R = R + delta_D;% 对于新区域的第一层，从上一个区域的最后一层的 R 增加一个 delta_D
h=h0+(i-1)*0.5;
% 检查 R 的条件
if R > (350+L0)
break;
end
% 使用余弦定理计算 theta_D
cos_theta_D = (2*R^2 - delta_D^2) / (2*R^2);
theta_D = acos(cos_theta_D);
ND = floor(2*pi/theta_D);
theta_D = 2*pi/ND; % 重新计算 theta_D 的值
% 计算第一层的坐标
cos_values = cos((0:ND-1) * theta_D);
sin_values = sin((0:ND-1) * theta_D);
x_values = R * cos_values;
y_values = R * sin_values;
h_values=h*ones(ND,1);
coordinates(idx:idx+ND-1, 1:2) = [x_values', y_values'];
coordinates(idx:idx+ND-1, 3) = h_values;
idx = idx + ND;
% 保存第一层坐标的距离
distance_first_layer = sqrt((x_values(2) - x_values(1))^2 + (y_values(2) -
y_values(1))^2);
% 其他层的坐标
layer = 2;
while true
prev=sqrt((coordinates(idx-1,1) - coordinates(idx-2,1))^2 + (coordinates(idx-1,2)-
coordinates(idx-2,2))^2);
R = R + sqrt(delta_D^2 - (prev/2)^2);
% 检查 R 的条件
if R > (350+L0)%%%%
break;
end
theta_D0 = atan2(coordinates(idx-ND,2), coordinates(idx-ND,1)) + theta_D/2;
cos_values = cos((0:ND-1) * theta_D + theta_D0); % 角平分线
sin_values = sin((0:ND-1) * theta_D + theta_D0); % 角平分线
x_values = R * cos_values;
y_values = R * sin_values;
coordinates(idx:idx+ND-1, 1:2) = [x_values', y_values'];
coordinates(idx:idx+ND-1, 3) = h_values;
idx = idx + ND;
% 检查距离条件
distance = sqrt((x_values(2) - x_values(1))^2 + (y_values(2) - y_values(1))^2);
if distance > 1.2* distance_first_layer
break;
end
layer = layer + 1;
end
i = i + 1;
end
% 删除未使用的空间
coordinates = coordinates(~isnan(coordinates(:,1)), :);
%筛选定日镜坐标
% 检查约束条件
theta_D = 2*pi/ND; % 重新计算 theta_D 的值
% 计算第一层的坐标
cos_values = cos((0:ND-1) * theta_D);
sin_values = sin((0:ND-1) * theta_D);
x_values = R * cos_values;
y_values = R * sin_values;
h_values=h*ones(ND,1);
coordinates(idx:idx+ND-1, 1:2) = [x_values', y_values'];
coordinates(idx:idx+ND-1, 3) = h_values;
idx = idx + ND;
% 保存第一层坐标的距离
distance_first_layer = sqrt((x_values(2) - x_values(1))^2 + (y_values(2) -
y_values(1))^2);
% 其他层的坐标
layer = 2;
while true
prev=sqrt((coordinates(idx-1,1) - coordinates(idx-2,1))^2 + (coordinates(idx-1,2)-
coordinates(idx-2,2))^2);
R = R + sqrt(delta_D^2 - (prev/2)^2);
% 检查 R 的条件
if R > (350+L0)%%%%
break;
end
theta_D0 = atan2(coordinates(idx-ND,2), coordinates(idx-ND,1)) + theta_D/2;
cos_values = cos((0:ND-1) * theta_D + theta_D0); % 角平分线
sin_values = sin((0:ND-1) * theta_D + theta_D0); % 角平分线
x_values = R * cos_values;
y_values = R * sin_values;
coordinates(idx:idx+ND-1, 1:2) = [x_values', y_values'];
coordinates(idx:idx+ND-1, 3) = h_values;
idx = idx + ND;
% 检查距离条件
distance = sqrt((x_values(2) - x_values(1))^2 + (y_values(2) - y_values(1))^2);
if distance > 1.2* distance_first_layer
break;
end
layer = layer + 1;
end
i = i + 1;
end
% 删除未使用的空间
coordinates = coordinates(~isnan(coordinates(:,1)), :);
%筛选定日镜坐标
% 检查约束条件
N = size(coordinates,1);
valid_indices = [];
for i = 1:N
x1 = coordinates(i,1);
y1 = coordinates(i,2);
if x1^2 + (y1-L0)^2 <= 350^2
valid_indices = [valid_indices; i];
end
end
data=[coordinates(valid_indices,1), coordinates(valid_indices,2)];
error2 = check_distance(data);
if ~error2
error('定日镜距离不符合');
end
data(:,3) = sqrt(data(:,1).^2 + data(:,2).^2);
data(:,4) = coordinates(valid_indices,3);
data = sortrows(data, 3);
end
function result = check_distance(coordinates)% 检验点是否满足两个定日镜之间的距离要大于等于 delta_D
% 给定的参数
global W_a;
delta_D = 5 + W_a;
% 计算所有点之间的距离
distances = pdist(coordinates);
% 找到最小距离
min_distance = min(distances);
% 检查最小距离是否大于等于 delta_D
if min_distance >= delta_D
result = true;
else
result = false;
end
end
N = size(coordinates,1);
valid_indices = [];
for i = 1:N
x1 = coordinates(i,1);
y1 = coordinates(i,2);
if x1^2 + (y1-L0)^2 <= 350^2
valid_indices = [valid_indices; i];
end
end
data=[coordinates(valid_indices,1), coordinates(valid_indices,2)];
error2 = check_distance(data);
if ~error2
error('定日镜距离不符合');
end
data(:,3) = sqrt(data(:,1).^2 + data(:,2).^2);
data(:,4) = coordinates(valid_indices,3);
data = sortrows(data, 3);
end
function result = check_distance(coordinates)% 检验点是否满足两个定日镜之间的距离要大于等于 delta_D
% 给定的参数
global W_a;
delta_D = 5 + W_a;
% 计算所有点之间的距离
distances = pdist(coordinates);
% 找到最小距离
min_distance = min(distances);
% 检查最小距离是否大于等于 delta_D
if min_distance >= delta_D
result = true;
else
result = false;
end
end
扫一扫
410
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
