# Python气象辐射光谱能量平衡模型

{% embed url="https://viadean.notion.site/Python-10a1ae7b9a3280d487d9f64519b462dc" %}

{% content-ref url="../zhi-dian-mi-jin-brief" %}
[zhi-dian-mi-jin-brief](../zhi-dian-mi-jin-brief)
{% endcontent-ref %}

### 🎯要点

1. 根据温室模型，计算不同情景下辐射通量和评估能量平衡，构建复杂温室模型计算
2. 计算和绘图大气、海洋、陆地表面和海冰复合模型数据
3. 建立简单能量平衡情景模型，并根据模型计算释放温度和时滞，计算并绘制地面辐射和吸收光谱
4. 计算和模拟非散射辐射传输，求解史瓦西微分方程

### 🍪语言内容分比

{% @mermaid/diagram content="pie title 语言分比
 "Python":90
 "C/C++":30" %}

{% @mermaid/diagram content="pie title 内容分比
 "算法模型":90
 "气候学":80
 "数学,微分方程":50
 "物理学":40
 "环境学,温室效应":40
 "热传递和辐射":30" %}

### 🍇Python微分方程

如果因变量具有恒定的变化率：

$$
\frac{d y}{d t}=C
$$

其中 C 是某个常数，您可以在 f 函数中提供微分方程，然后使用以下代码使用此模型计算答案。该代码假设 0 和 10 之间有 100 个均匀间隔的时间， y 的初始值为 6 ，变化率为 1.2 ：

```python
 import numpy as np
 import matplotlib.pyplot as plt
 from scipy.integrate import solve_ivp
 from ode_helpers import state_plotter
 ​
 def f(t, y, c):
     dydt = [c[0]]
     return dydt
 ​
 tspan = np.linspace(0, 10, 100)
 yinit = [6]
 c = [1.2]
 ​
 sol = solve_ivp(lambda t, y: f(t, y, c), 
                 [tspan[0], tspan[-1]], yinit, t_eval=tspan, rtol = 1e-5)
 ​
 state_plotter(sol.t, sol.y, 1)
```

如果因变量的变化率是时间的函数，则可以轻松编码。例如，如果微分方程是某个二次函数，如下所示：

$$
\frac{d y}{d t}=\alpha t^2+\beta t+\gamma
$$

您可以使用此模型通过以下代码计算答案，假设 y 的初始值 0 到 4 之间有 20 个均匀间隔的时间 𝑦 是 6，多项式由向量 \[2, -6, 3] 定义：

```python
 import numpy as np
 import matplotlib.pyplot as plt
 from scipy.integrate import solve_ivp
 from ode_helpers import state_plotter
 ​
 def f(t, y, c):
     dydt = np.polyval(c, t)
     return dydt
     
 tspan = np.linspace(0, 4, 20)
 yinit = [6]
 c = [2, -6, 3]
 ​
 sol = solve_ivp(lambda t, y: f(t, y, c), 
                 [tspan[0], tspan[-1]], yinit, t_eval=tspan, rtol = 1e-5)
 ​
 state_plotter(sol.t, sol.y, 1)
```

对于增长型，变化率取决于数量以及一些比例常数：

$$
\frac{d y}{d t}=C \cdot y
$$

其中 C 又是某个常数。以下代码将计算上述模型的 3 秒内 25 个点的增长量，初始为 10，比例常数为 1.02：

```python
 import numpy as np
 import matplotlib.pyplot as plt
 from scipy.integrate import solve_ivp
 from ode_helpers import state_plotter
 ​
 def f(t, y, c):
     dydt = [c[0] * y[0]]
     return dydt
 ​
 tspan = np.linspace(0, 3, 25)
 yinit = [10]
 c = [1.02]
 ​
 sol = solve_ivp(lambda t, y: f(t, y, c), 
                 [tspan[0], tspan[-1]], yinit, t_eval=tspan, rtol = 1e-5)
 ​
 state_plotter(sol.t, sol.y, 1)
```

通过确保微分函数返回每个变量的值，可以解决多变量系统。例如，在以下系统中，第一个变量的变化率仅取决于时间，而第二个变量则取决于时间和第一个变量：

$$
\frac{d y_0}{d t}=\alpha \cos (\beta t) \quad \frac{d y_1}{d t}=\gamma y_0+\delta t
$$

该系统的微分函数 f 将有一个 2 元素列表作为输出。另外，如果您的系统具有多个因变量，请务必将初始条件放入列表中。例如，系统定义为：

$$
\frac{d y_0}{d t}=4 \cos (3 t) \quad \frac{d y_1}{d t}=-2 y_0+0.5 t
$$

您可以使用以下脚本来求解 $$y_0$$ 和 $$y_1$$；代码假设$$y_0$$ 从 0 开始，$$y_1$$ 从 -3 开始

```python
 import numpy as np
 import matplotlib.pyplot as plt
 from scipy.integrate import solve_ivp
 from ode_helpers import state_plotter
 ​
 def f(t, y, c):
     dydt = [c[0]*np.cos(c[1]*t), c[2]*y[0]+c[3]*t]
     return dydt
 ​
 tspan = np.linspace(0, 5, 100)
 yinit = [0, -3]
 c = [4, 3, -2, 0.5]
 ​
 sol = solve_ivp(lambda t, y: f(t, y, c), 
                 [tspan[0], tspan[-1]], yinit, t_eval=tspan, rtol = 1e-5)
 ​
 state_plotter(sol.t, sol.y, 1)
```

系统必须仅以一阶微分方程的形式来写出。要求解具有高阶导数的系统，您首先要编写一个简单的一阶方程级联系统，然后在微分函数中使用它们。例如，假设您有一个以恒定抖动为特征的系统：

$$
j=\frac{d^3 y}{d t^3}=C
$$

首先要做的是写出三个一阶微分方程来表示三阶方程：

$$
\begin{aligned} & y_0=y \\ & \frac{d y_0}{d t}=\frac{d y}{d t}=y_1 \quad \longrightarrow \quad \frac{d y_0}{d t}=y_1 \\ & \frac{d y_1}{d t}=\frac{d^2 y_0}{d t^2}=\frac{d^2 y}{d t^2}=y_2 \quad \longrightarrow \quad \frac{d y_1}{d t}=y_2 \\ & \frac{d y_2}{d t}=\frac{d^2 y_1}{d t^2}=\frac{d^3 y_0}{d t^3}=\frac{d^3 y}{d t^3}=j=C \quad \longrightarrow \quad \frac{d y_2}{d t}=C \end{aligned}
$$

请注意导数如何级联，以便现在可以将恒定加加速度方程写为一组三个一阶方程。请注意，在该系统中，y\[0] 表示位置，y\[1] 表示速度，y\[2] 表示加速度。这种类型的级联系统在运动方程建模时经常出现。

```python
 import numpy as np
 import matplotlib.pyplot as plt
 from scipy.integrate import solve_ivp
 from ode_helpers import state_plotter
 ​
 def f(t, y, c):
     dydt = [y[1], y[2], c[0]]
     return dydt
 ​
 tspan = np.linspace(0, 8, 50)
 yinit = [6, 2, -4]
 c = [1.3]
 ​
 sol = solve_ivp(lambda t, y: f(t, y, c), 
                 [tspan[0], tspan[-1]], yinit, t_eval=tspan, rtol = 1e-5)
 ​
 state_plotter(sol.t, sol.y, 1)
```
