import matplotlib.pyplot as plt
import numpy as np
import math
fig = plt.figure() # 建立图形窗口
ax = fig.gca(projection='3d') # 使用三维坐标
w25 = 5.880 * 10 ** (-8)
w26 = 1.781 * 10 ** (-8)
w27 = 3.635 * 10 ** (-8)
w28 = 2.468 * 10 ** (-7)
w29 = -1.227 * 10 ** (-8)
w30 = -2.883 * 10 ** (-9)
w31 = -1.361 * 10 ** (-6)
w32 = 5.844 * 10 ** (-7)
w33 = -1.199 * 10 ** (-7)
w34 = 3.169 * 10 ** (-9)
w9 = -1.801 * 10 ** (-5)
w10 = -4.832 * 10 ** (-6)
w11 = -3.321 * 10 ** (-6)
w12 = -5.468 * 10 ** (-5)
w13 = 2.610 * 10 ** (-5)
w14 = -1.919 * 10 ** (-6)
w19 = -3.055 * 10 ** (-4)
w20 = -8.582 * 10 ** (-5)
w21 = 3.410 * 10 ** (-5)
w22 = 1.372 * 10 ** (-4)
w23 = -2.527 * 10 ** (-5)
w24 = -1.084 * 10 ** (-6)
w1 = 3.658 * 10 ** (-1)
w2 = 1.223 * 10 ** (-3)
w3 = -1.452 * 10 ** (-3)
w4 = 4.107 * 10 ** (-4)
w6 = 4.411 * 10 ** (-2)
w7 = 2.410 * 10 ** (-3)
w8 = -9.462 * 10 ** (-4)
w16 = -1.999 * 10 ** (-2)
w17 = 6.975 * 10 ** (-3)
w18 = 3.062 * 10 ** (-3)
w5 = -5.268 * 10 ** (-1)
w15 = -4.708 * 10 ** (-1)
w0 = -3.573 * 10 ** (-2)
l = 300
Ld = 60
Qd = 1.5
n = np.arange(0.1, 0.9, 0.01)
Lt = np.arange(10, 30, 0.25)
n, Lt = np.meshgrid(n, Lt)
# n = np.arange(0.2, 0.3, 0.01)
# Lt = np.arange(23, 25, 0.2)
H = (w25 * (Ld ** 3) + w26 * (Ld ** 2) * Lt + w27 * (
Ld ** 2) * l + w28 * Ld * Lt ** 2 + w29 * Ld * Lt * l + w30 * Ld * l ** 2 + w31 * Lt ** 3 + w32 * Lt ** 2 * l + w33 * Lt * l ** 2 + w34 * l ** 3) * Lt
M = (w9 * (Ld ** 2) + w10 * Ld * Lt + w11 * Ld * l + w12 * Lt ** 2 + w13 * Lt * l + w14 * l ** 2) * Lt + (
w19 * (
Ld ** 2) + w20 * Ld * Lt + w21 * Ld * l + w22 * Lt ** 2 + w23 * Lt * l + w24 * l ** 2) * Lt * Qd / Ld
P = (w2 * Ld + w3 * Lt + w4 * l) * Lt + (w6 * Ld + w7 * Lt + w8 * l) * Lt * Qd / Ld + (
w16 * Ld + w17 * Lt + w18 * l) * Lt * Qd ** 2 / Ld ** 2
S = w1 * Lt * Qd / Ld + (w5 * Lt * Qd ** 2 / Ld ** 2) + w15 * Lt * Qd ** 3 / Ld ** 3 + w0 * Lt
z = H * n ** (3 / 2) + M * n + P * n ** (1 / 2) + S
#X, Y = np.meshgrid(Lt, n)
Z = np.array(z)
# Z_result = []
# for i in range(0, 10):
# for j in range(1):
# if j == 0:
# Z_result.append([])
# Z_result[i].append(Z[i])
#
# Z_result1 = np.array(Z_result)
# print(Z_result)
# print(X)
# print(Y)
# print(Z)
#print(max(Z))
ax.plot_surface(n, Lt, Z) # 绘制
ax.set_xlabel('${ n_p}$ [$10^{16} cm^{-3}]$')
ax.set_ylabel('$L_t$')
ax.set_zlabel('$Q_t$ [nC]')
# ax.view_init(20, 45)
plt.show()