月平均余弦效率求解

import sympy as sp
import pandas as pd
import math
import numpy as np
import csv
from tqdm import tqdm

#D为给定日期从春分算起的天数
D = [306, 337, 0, 31, 61, 92, 122, 153, 184, 214, 245, 275]
#ST为当地时间
ST = [9, 10.5, 12, 13.5, 15]
#z为集热器中心距离任意一个定日镜的竖直距离，为定值
z = 76

delta = []
for d in D:
    delta.append(math.asin(math.sin(2*math.pi*d/365) * math.sin(2*math.pi*23.45/360)))

W = []
for st in ST:
    W.append(math.pi/12*(st-12))
latitude = math.radians(39.4)

Alpha = []
Gamma = []
for i in range(len(delta)):
    for j in range(len(W)):
        alpha = math.asin(math.cos(delta[i])*math.cos(latitude)*math.cos(W[j]) + math.sin(delta[i])*math.sin(latitude))
        Alpha.append(alpha)
        gamma = math.acos(
            (math.sin(delta[i])-math.sin(alpha)*math.sin(latitude)) / (math.cos(alpha)*math.cos(latitude)) + 1e-6
        )
        Gamma.append(gamma)

e_in = []   #入射光线的单位方向向量
for k in range(12 * 5):
    e_in.append(sp.Matrix([
        -math.cos(Alpha[k])*math.sin(Gamma[k]), -math.cos(Alpha[k])*math.cos(Gamma[k]), -math.sin(Alpha[k])
    ]))

excel = pd.read_excel("C:\\Users\\ChenCong\\Desktop\\test.xlsx")

e_out = []  #反射光线的单位方向向量
X = list(excel['x'])
Y = list(excel['y'])
for m in range(len(excel)):
    e_out.append(sp.Matrix([-X[m], -Y[m], z]) / math.sqrt((X[m])**2 + (Y[m])**2 + z**2))
#print(e_out)
#print(e_in[0].reshape(1, 3) * e_out[0])
sum = sp.Matrix([0])
monthly_average_cos_efficiency = []
for k in tqdm(range(12)):
    for i in range(5*k, 5 + 5*k):
        for j in range(1745):
            sum = sum + ((sp.Matrix([1]) - e_in[i].reshape(1, 3) * e_out[j])
                         / math.sqrt((((e_out[j] - e_in[i]).reshape(1, 3)) * (e_out[j] - e_in[i]))[0]))
    monthly_average_cos_efficiency.append(sum[0]/(5*1745))
    sum = sp.Matrix([0])
#print((((e_out[j] - e_in[i]).reshape(1, 3)) * (e_out[j] - e_in[i]))[0])
print(monthly_average_cos_efficiency)
print(np.mean(monthly_average_cos_efficiency))

f = open('table_1.csv', 'a', newline='', encoding='utf-8-sig')
global writer
writer = csv.writer(f)
writer.writerow(['时间', '月平均效率'])
for i in range(12):
    writer.writerow([f'{i+1}月21日', monthly_average_cos_efficiency[i]])