import os
import numpy as np
import matplotlib.pyplot as plt

np.random.seed(42)

# -----------------------------
# 参数
# -----------------------------
A = 1.0                 # m^2
N = 100                 # 灯管数量
lamp_length = 0.5       # m
P_nominal = 9.0         # kW
rows = 2
lamps_per_row = 50
H = 0.12                # m，安装高度（工程假设）
grid_n = 160
rays_per_lamp = 4000
n_mc = 200_000

# -----------------------------
# 蒙特卡洛：总有效热流密度
# -----------------------------
P_lamp = np.random.triangular(left=8.5, mode=9.0, right=9.0, size=n_mc)  # kW
eta_r = np.random.uniform(0.50, 0.60, size=n_mc)
eta_a = np.random.uniform(0.80, 0.95, size=n_mc)

q_theoretical = N * P_nominal / 1000.0
q_incident_samples = N * P_lamp * eta_r / 1000.0
q_effective_samples = N * P_lamp * eta_r * eta_a / 1000.0

# -----------------------------
# 射线追踪：空间分布
# -----------------------------
x_row_centers = np.array([0.25, 0.75])
y_centers = np.linspace(0.01, 0.99, lamps_per_row)

counts = np.zeros((grid_n, grid_n), dtype=np.int64)

for xc in x_row_centers:
    for yc in y_centers:
        n = rays_per_lamp
        x0 = xc + np.random.uniform(-lamp_length / 2, lamp_length / 2, size=n)
        y0 = np.full(n, yc)

        mu = np.random.uniform(1e-6, 1.0, size=n)
        phi = np.random.uniform(0.0, 2.0 * np.pi, size=n)
        sin_theta = np.sqrt(1.0 - mu ** 2)

        dx = sin_theta * np.cos(phi)
        dy = sin_theta * np.sin(phi)

        x_hit = x0 + dx * (H / mu)
        y_hit = y0 + dy * (H / mu)

        mask = (x_hit >= 0.0) & (x_hit < 1.0) & (y_hit >= 0.0) & (y_hit < 1.0)
        x_valid = x_hit[mask]
        y_valid = y_hit[mask]

        xi = np.floor(x_valid * grid_n).astype(int)
        yi = np.floor(y_valid * grid_n).astype(int)
        np.add.at(counts, (yi, xi), 1)

cell_area = (1.0 / grid_n) * (1.0 / grid_n)
p_cell = counts / counts.sum()
weight_map = p_cell / cell_area

q_eff_mean = float(np.mean(q_effective_samples))
q_map = q_eff_mean * weight_map

# -----------------------------
# 保存图像
# -----------------------------
plt.rcParams["font.sans-serif"] = ["DejaVu Sans"]
plt.rcParams["axes.unicode_minus"] = False

out_dir = "/mnt/data"
cloud_path = os.path.join(out_dir, "石英灯热流密度云图_代码输出.png")
hist_path = os.path.join(out_dir, "石英灯有效热流密度分布_代码输出.png")

fig = plt.figure(figsize=(8.5, 7.2))
ax = fig.add_subplot(111)
im = ax.imshow(q_map, extent=[0, 1, 0, 1], origin="lower", aspect="equal")
fig.colorbar(im, ax=ax, label="Heat flux density (MW/m²)")
ax.set_title("Heat flux density map")
ax.set_xlabel("X / m")
ax.set_ylabel("Y / m")
fig.tight_layout()
fig.savefig(cloud_path, dpi=300, bbox_inches="tight")
plt.close(fig)

fig = plt.figure(figsize=(8.2, 5.8))
ax = fig.add_subplot(111)
ax.hist(q_effective_samples, bins=60)
ax.set_title("Monte Carlo distribution of effective heat flux density")
ax.set_xlabel("Effective heat flux density (MW/m²)")
ax.set_ylabel("Count")
fig.tight_layout()
fig.savefig(hist_path, dpi=300, bbox_inches="tight")
plt.close(fig)

print("理论最大热流密度:", round(q_theoretical, 3), "MW/m²")
print("平均有效热流密度:", round(q_eff_mean, 3), "MW/m²")
print("代码已运行并输出图像：")
print(cloud_path)
print(hist_path)