import numpy as np
import pandas as pd
from math import cos,sin,acos,asin,pi,exp,sqrt
import time
from tqdm import tqdm
def find_mirror(x0,y0,position_data,D_max=29.16):
"""
计算目标定日镜区域内所有符合条件定日镜的坐标
:参数 x0，y0: 目标定日镜坐标
:return: 所有符合条件定日镜的坐标
(sqrt((x-x0)**2+(y-y0) ** 2)>0.5)
"""
#position_data的第一项不是坐标数据，因此从第二项开始
position_data_np = position_data[1:]
a=((position_data_np[0] - x0) ** 2 + (position_data_np[1] - y0) ** 2).astype(float)
distances = np.array(np.sqrt(a))
within_mirrors = position_data_np[(distances <= D_max) & (distances >0.1)]
return within_mirrors
def F(L,W,H,D_0,ST_0,delta_t,HO,H1,HR,p1,x_c,y_c):
"""参数：L,W,H:定日镜的长、宽、安装高度
D_0:日期数据
ST_0:每天的时间数据
delta_t：计算阴影的步长
H0:吸收塔高度
H1:集热器高度
HR:集热器半径
p1:定日镜场的坐标数据
x_c,y_c:吸收塔中心
"""
position_data = pd.read_excel('./附件.xlsx', header=None)
N=(W/delta_t+1)**2*4
yita_ref=0.92
S=L*W
new_column=np.array([H for i in range(len(p1))])#安装高度4m
print(new_column)
P=np.column_stack((p1,new_column))
# 5*1745
Yita =np.zeros([len(ST_0),len(P)])
Yita_trunc =np.zeros([len(ST_0),len(P)])
Yita_at=np.zeros([len(ST_0),len(P)])
Yita_cos = np.zeros([len(ST_0), len(P)])
Yita_sb = np.zeros([len(ST_0), len(P)])
# 12*5
E_0 = np.zeros([len(D_0), len(ST_0)])
Y1 = np.zeros([len(D_0), len(ST_0)])
Y2 = np.zeros([len(D_0), len(ST_0)])
Y3 =np.zeros([len(D_0), len(ST_0)])
Y4 = np.zeros([len(D_0), len(ST_0)])
# 12*5
DF = np.zeros([len(D_0), 5])
for di, D in enumerate(D_0):
for sti, ST in enumerate(ST_0): # 计算太阳位置以及相关参数
fai = 39.4 * pi / 180;
delta = asin(sin(2 * pi * D / 365) * sin(2 * pi * 23.45 / 360))
w = (pi / 12) * (ST - 12)
# 太阳高度角
sin_as = cos(delta) * cos(fai) * cos(w) + sin(delta) * sin(fai)
cos_as=sqrt(1-sin_as**2)
# 太阳方位角
cos_rs = (sin(delta) - sin_as * sin(fai))/(sqrt(1 - sin_as ** 2) * cos(fai))
if abs(cos_rs)>1:
cos_rs = cos_rs / abs(cos_rs)
if ST <= 12:
sin_rs = sqrt(1 - cos_rs ** 2)
else:
sin_rs = -sqrt(1 - cos_rs ** 2)
a = 0.4237 - 0.00821 * (6 - 3) ** 2
b = 0.5055 + 0.00595 * (6.5 - 3) ** 2
c = 0.2711 + 0.01858 * (2.5 - 3) ** 2
DNI = 1.366 * (a + b * exp(-c / sin_as))
A0 = np.array([x_c, y_c, HO]) # 集热器中心
Ls = np.array([-cos_as *sin_rs, -cos_as * cos_rs, -sin_as]) # 入射光线单位向量
print("Ls",Ls)
s=np.array([1,0,0])
for i in tqdm(range(len(P))):
within_mirrors = np.array(find_mirror(P[i][0], P[i][1], position_data))
# A镜中点 反射向量
Di=P[i]
Lr=A0-Di
nl = -Ls + Lr / np.linalg.norm(Lr)# 定日镜的的法向量
nl=nl/np.linalg.norm(nl)#法向量的单位向量
print("nl",nl)
print("np.linalg.norm(nl)",np.linalg.norm(nl))
beta1 = asin(nl.dot([0, 0, 1])/np.linalg.norm(nl))
print("alpha_A",beta1)# 俯仰角
n0 = np.array([nl[0], nl[1], 0])
if nl[1] >= 0:
alpha1 = acos(n0.dot(s) / np.linalg.norm(n0))# 方位角
else:
alpha1 = -acos(n0.dot(s) / np.linalg.norm(n0))
print("gamma_A", alpha1)
# A的旋转矩阵
Ta = np.array([
[cos(alpha1) * cos(pi / 2 - beta1), -sin(alpha1), cos(alpha1) * sin(pi / 2 - beta1)],
[sin(alpha1) * cos(pi / 2 - beta1), cos(alpha1), sin(alpha1) * sin(pi / 2 - beta1)],
[-sin(pi / 2 - beta1), 0, cos(pi / 2 - beta1)]
])
print("Ta",Ta)
light = 0
empty = 0
barr_tower = 0
barr_s = 0
barr_r = 0
# 遍历每一个点
for dx in np.arange(-W / 2, W / 2 + 0.1, delta_t):
for dy in np.arange(-L / 2, L / 2 + 0.1, delta_t):
Dxy = np.array([dx, dy,0]) # A镜上的某个点在A镜坐标系
# print("Dxy",Dxy)
Di_d=Ta.dot(Dxy)+Di #A镜上的点转置到地面坐标系
print("Di_d",Di_d)
################遍历每一个入射光线圆锥的光线##################
the1 = 0.002
for the2 in np.arange(0, 2 * pi, pi / 2):
#if_barr=0,假设不会造成阴影遮挡损失
if_barr = 0
# g是入射光线在主光线锥体系的坐标
g=np.array([sin(the1)*cos(the1),sin(the1)*sin(the2),cos(the1)])
# 根据入射主光线，计算主光线锥体系的旋转矩阵Ls:入射的太阳主光线
v = pi / 2 - acos(Ls.dot(np.array([0, 0, 1]))/np.linalg.norm(Ls))
nl_gO = np.array([Ls[0], Ls[1], 0])
if Ls[1] >= 0:
u = acos(nl_gO.dot(s) / np.linalg.norm(nl_gO))
else:
u = -acos(nl_gO.dot(s) / np.linalg.norm(nl_gO))
T_s = np.array([
[cos(u) * cos(pi / 2 - v), -sin(u), cos(u) * sin(pi / 2 - v)],
[sin(u) * cos(pi / 2 - v), cos(u), sin(u) * sin(pi / 2 - v)],
[-sin(pi / 2 - v), 0, cos(pi / 2 - v)]
])
g_d = T_s.dot(g) # 转置到地面坐标系g;入射光锥中的某条入射线
g_r=g_d-2*g_d.dot(nl)*(nl)#对应的反射向量#nl:A的法向量
##################一、 判断入射光线是否被塔遮挡##########H#######
a, b, c = g_d
x0, y0, z0 = Di_d
delta_tower=4 * (a * (x0 - x_c) + b * (y0 - y_c)) ** 2 - 4 * (a ** 2 + b ** 2) * ((x0 -
x_c) ** 2 + (y0 - y_c) ** 2 - HR ** 2)
if delta_tower >= 0:
t1=(-2 * (a * (x0-x_c)+b * (y0-y_c))+sqrt(delta_tower))/(2 * (a ** 2+b ** 2))
t2=(-2 * (a * (x0-x_c)+b * (y0-y_c))-sqrt(delta_tower))/(2 * (a ** 2+b ** 2))
if min(t1 * c + z0, t2 * c + z0) <= (HO + H1 / 2) and min(t1 * c + z0,t2 * c + z0) >= 0:
barr_tower += 1
continue
##################二、判断入射光线是否被其他光镜遮挡##################
for j in range(within_mirrors.shape[0]):
# B的中心坐标alpha,beta直接调取提前计算的
B=within_mirrors[j]
B = np.append(B,H)
Lrb=A0-B
nlb = -Ls + Lrb / np.linalg.norm(Lrb) # 法向量
beta2 = asin(nlb.dot([0, 0, 1]) / np.linalg.norm(nlb)) # 俯仰角
nOb = np.array([nlb[0], nlb[1], 0])
if nlb[1] >= 0:
alpha2 = acos(nOb.dot(s) / np.linalg.norm(nOb)) # 方位角
else:
alpha2 = -acos(nOb.dot(s) / np.linalg.norm(nOb))
Tb = np.array([
[cos(alpha2) * cos(pi / 2 - beta2), -sin(alpha2), cos(alpha2) * sin(pi / 2 - beta2)],
[sin(alpha2) * cos(pi / 2 - beta2), cos(alpha2), sin(alpha2) * sin(pi / 2 - beta2)],
[-sin(pi / 2 - beta2), 0, cos(pi / 2 - beta2)]
])
Di_b = Tb.T.dot(Di_d - B) # A镜上的点从地面坐标系→B镜坐标系
g_b=Tb.T.dot(g_d) #A入射光线从地面坐标系→B镜坐标系
t = -Di_b[2]/g_b[2]# 计算入线光线的在B镜坐标系的交点
x_b = g_b[0] * t + Di_b[0]
y_b = g_b[1] * t + Di_b[1]
D0 = np.array([x_b, y_b, 0])
D0 = Tb.dot(D0) + B# 交点转到地面
if abs(x_b) <= W / 2 and abs(y_b) <= L / 2 and D0[2] > Di_d[2]: # 如果被遮挡了，就直接算下一条入射的线
barr_s += 1
if_barr = 1
break
##################三、对应的反射##################
g_r = g_d - 2 * g_d.dot(nl) * (nl) # 对应的反射向量#nl:A的法向量
g_r_b=Tb.T.dot(g_r)#转置到b镜面坐标系
t = -Di_b[2] / g_r_b[2] # 计算反射光线的交点
x_b=g_r_b[0]*t+Di_b[0]
y_b=g_r_b[1]*t+Di_b[1]
D0 = np.array([x_b, y_b, 0])
D0 = Tb.dot(D0) + B # 交点转到地面
if abs(x_b) <= W / 2 and abs(y_b) <= L / 2 and D0[2] > Di_d[2]:
barr_r += 1
if_barr=1
break
##################四、是否吸收##################
if if_barr==0:
a, b, c = g_r
x0, y0, z0 = Di_d
delta_recieve = 4 * (a * (x0-x_c) + b * (y0-y_c)) ** 2 - 4 * (a ** 2 + b ** 2) * ((x0-
x_c) ** 2 + (y0 - y_c) ** 2 - HR ** 2)
if delta_recieve >= 0:
t1 = (-2 * (a * (x0 - x_c) + b * (y0 - y_c)) + sqrt(delta_recieve))/(2 * (a ** 2 + b ** 2))
t2=(-2 * (a * (x0 - x_c) + b * (y0 - y_c)) - sqrt(delta_recieve)) /(2 * (a ** 2 + b ** 2))
if min(t1 * c + z0, t2 * c + z0) <= (HO + H1 / 2) and min(t1 * c + z0, t2 * c + z0) >= (HO - H1 / 2):
light += 1
else:
empty += 1
# 这里是这个时间点，计算第i个面板的数值，并记录，每个列表是1745的长度
yita_sb=1-(barr_r+barr_s+barr_tower)/N
yita_cos=abs(Ls.dot(-nl)/np.linalg.norm(Ls));
print("yita_cos",yita_cos)
HR0=np.linalg.norm(Lr);
yita_at = 0.99321 - 0.0001176 * HR0 + 1.97e-8 * (HR0 ** 2);
if N - barr_s - barr_r - barr_tower == 0:
yita_trunc = 1
else:
yita_trunc = (light)/(N - barr_s - barr_r - barr_tower)
yita = yita_sb * yita_cos * yita_at * yita_trunc * yita_ref;
# print(empty)
# print(yita_sb,yita_cos,yita_at,yita_trunc,yita)
#5*1745
Yita_sb[sti, i] = yita_sb
Yita_cos[sti, i] = yita_cos
Yita_at[sti, i] = yita_at
Yita_trunc[sti, i] = yita_trunc
Yita[sti, i] = yita
# 这里是在D,ST的循环里，计算第sti个时间点
E_0[di, sti] = DNI*sum(S*Yita[sti, :])
Y1[di, sti] = np.mean(Yita[sti, :])
Y2[di, sti] = np.mean(Yita_cos[sti, :])
Y3[di, sti] = np.mean(Yita_sb[sti,:])
Y4[di, sti] = np.mean(Yita_trunc[sti, :])
# print(Y1[di, sti],Y2[di, sti],Y3[di, sti],Y4[di, sti])
# 已经计算完一个具体时间点
DF[di,0]=np.mean(Y1[di,:])
DF[di,1]=np.mean(Y2[di,:])
DF[di,2]=np.mean(Y3[di,:])
DF[di,3]=np.mean(Y4[di,:])
DF[di,4]=sum(E_0[di,:])/(len(P)*S)/len(ST_0)
return DF
#导入题目的附件并添加乙值
P=pd.read_excel(r'附件.xlsx').values
# Dis=np.load(distance.npy)
D_0=[306,337,0,31,61,92,122,153,184,214,245,275]
ST_0=[9,10.5,12,13.5,15]
# D_0=[0]
delta_t=1.5
x_c=0;y_c=0;L=W=6;HO=80;H1=8;H=4;
HR=3.5
start_time = time.time()
# L,W,H,D_0,ST_0,delta_t,HO,HI,HR,p1,x_c,y_c
Wp = F(L, W, H, D_0, ST_0,delta_t,HO, H1, HR, P, x_c, y_c)
print(Wp)
# 记录程序运行时间
end_time =time.time()
elapsed_time_seconds = end_time - start_time
elapsed_time_minutes = round(elapsed_time_seconds / 60, 2)
print(f"代码运行时间：{elapsed_time_minutes}分钟")
# [[0.48722307 0.71993576 0.89510487 0.88293641 0.42453269]
# [0.50997021 0.74044035 0.91598281 0.87011053 0.4808807 ]
# [0.52629131 0.76114002 0.9218384 0.86245742 0.52352441]
# [0.54095191 0.77933997 0.922251 0.86109865 0.55672885]
# [0.55104133 0.78931698 0.9218086 0.86361524 0.57581782]
# [0.55449635 0.79235882 0.92160344 0.86488058 0.58190626]
# [0.55096519 0.78921185 0.92187393 0.86359106 0.57565267]
# [0.54035807 0.77863695 0.92229914 0.86105514 0.55547768]
# [0.52557774 0.76009198 0.92172951 0.86284416 0.52162896]
# [0.50755579 0.73783541 0.91382808 0.87183096 0.4746393 ]
# [0.48481803 0.71819561 0.89271748 0.88383722 0.41888792]
# [0.47425281 0.71108228 0.87946017 0.88953137 0.39437817]]
# [[0.48722307 0.71993576 0.89510487 0.88293641 0.42453269]
# [0.50997021 0.74044035 0.91598281 0.87011053 0.4808807 ]
# [0.52629131 0.76114002 0.9218384 0.86245742 0.52352441]
# [0.54095191 0.77933997 0.922251 0.86109865 0.55672885]
# [0.55104133 0.78931698 0.9218086 0.86361524 0.57581782]
# [0.55449635 0.79235882 0.92160344 0.86488058 0.58190626]
# [0.55096519 0.78921185 0.92187393 0.86359106 0.57565267]
# [0.54035807 0.77863695 0.92229914 0.86105514 0.55547768]
# [0.52557774 0.76009198 0.92172951 0.86284416 0.52162896]
# [0.50755579 0.73783541 0.91382808 0.87183096 0.4746393 ]
# [0.48481803 0.71819561 0.89271748 0.88383722 0.41888792]
# [0.47425281 0.71108228 0.87946017 0.88953137 0.39437817]]
该篇文章详细介绍了使用Python进行定日镜系统的设计，包括计算目标定日镜区域内的镜像位置，以及根据太阳位置、镜子角度和塔结构等因素分析阴影遮挡和能量吸收情况。函数`find_mirror`用于寻找符合条件的定日镜，`F`函数则实现了整个过程的模拟计算。
扫一扫
1867
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
