2023年高教社杯全国大学生数学建模竞赛A题

参考代码

import sympy as sp
import pandas as pd
import math
from math import sin, cos, tan, asin, acos, atan, pi, sqrt
import numpy as np
import csv
from tqdm import tqdm

# D为给定日期从春分算起的天数
# days_from_spring_equinoxs: 1 * 12 列表
days_from_spring_equinoxs = [306, 337, 0, 31, 61, 92, 122, 153, 184, 214, 245, 275]
np_days_from_spring_equinox = np.array(days_from_spring_equinoxs)

# local_time为当地时间，后面加s表示列表
# local_times = 1 * 5 列表
local_times = [9, 10.5, 12, 13.5, 15]
np_local_times = np.array([9, 10.5, 12, 13.5, 15])
# z为集热器中心距离任意一个定日镜的竖直距离，为定值
z = 76

# 计算太阳赤纬角：solar_declination_angle，和days_from_spring_equinoxs维度一样: 1 * 12 列表,
# 注意使用numpy的arcsin, sin, pi，而非math
solar_declination_angles = np.arcsin(np.sin(2 * np.pi * np_days_from_spring_equinox / 365) * np.sin(2 * np.pi * 23.45 /360))
# 计算太阳时角，和local_times维度一样，1 * 5 列表
solar_hour_angles = np.pi / 12 * (np_local_times - 12)

# latitude为当地纬度
latitude = math.radians(39.4)

# 计算太阳高度角solar_altitude_angles
temp = np.outer(np.cos(solar_declination_angles) * np.cos(latitude), np.cos(solar_hour_angles)) + \
       np.outer(np.sin(solar_declination_angles) * np.sin(latitude), np.ones(solar_hour_angles.shape))
# 先通过clip限定范围在[-1.0, 1.0]之间，然后计算反正弦, shape = (12, 5)
solar_altitude_angles = np.arcsin(np.clip(temp, -1.0, 1.0))
print(f'{solar_altitude_angles.shape = }')

# 计算太阳方位角solar_azimuth_angle, shape = (12, 5)
temp = np.divide(
    np.outer(np.sin(solar_declination_angles), np.ones(solar_hour_angles.shape)) - np.sin(solar_altitude_angles) * np.sin(latitude),
    np.cos(solar_altitude_angles) * np.cos(latitude)
)
solar_azimuth_angles = np.arccos(np.clip(temp, -1.0, 1.0))
print(f'{solar_azimuth_angles.shape = }')