雨滴谱gamma分布函数

数据：

 

 

 

目标主要是用求出N0、μ、Dm

 

 

把之前那个雨水含量的代码改一改

#!usr/bin/env python
# -*- coding:utf-8 -*-
"""
@author: Suyue
@file: raincontent.py
@time: 2024/05/23
@desc:
"""
import numpy as np
import pandas as pd
import xlwt
import math

df1 = pd.read_excel('G:/lianxi/20240510五原数浓度.xls')

N = df1['数浓度']
d = df1['粒径']
# 0阶矩
d_0 = pow(d,0)
# 3阶矩
d_3 = pow(d,3)
# 4阶矩
d_4 = pow(d,4)

delta_d = df1['变化直径']
# 计算0阶矩
D0 = N*d_0*delta_d
df2 = pd.DataFrame(D0)
D_sum_0 = df2.sum()

# 计算3阶矩
D3 = N*d_3*delta_d
df3 = pd.DataFrame(D3)
D_sum_3 = df3.sum()

# 计算4阶矩
D4 = N*d_4*delta_d
df4 = pd.DataFrame(D4)
D_sum_4 = df4.sum()
print(D_sum_0,D_sum_3,D_sum_4)

结果：

尝试用已知的M0、M3、M4求μ：

然后会发现中间差着阶呢，最后会出现μ的三次方这样，不好求μ的结果，所以得换挨着相近的阶距。

采用M2、M3、M4阶距计算μ：

求λ和Nw的方法

 

 计算形状因子u，斜率参数L,截距N0，gamma分布

 

#!usr/bin/env python
# -*- coding:utf-8 -*-
"""
@author: Suyue
@file: raincontent.py
@time: 2024/05/23
@desc:
"""
import numpy as np
import pandas as pd
import xlwt
import math
from scipy.special import gamma
import matplotlib
import matplotlib.pyplot as plt
import matplotlib.ticker as ticker
matplotlib.rc("font",family='YouYuan')

df1 = pd.read_excel('20240510五原数浓度.xls')

N = df1['数浓度']
d = df1['粒径']
# 粒径平方，2可以理解为2阶矩
d_2 = pow(d,2)
# 粒径3次方，3可以理解为3阶矩
d_3 = pow(d,3)
# 粒径4次方，4可以理解为4阶矩
d_4 = pow(d,4)

delta_d = df1['变化直径']
# 计算2阶矩
D2 = N*d_2*delta_d
df2 = pd.DataFrame(D2)
D_sum_2 = df2.sum()

# 计算3阶矩
D3 = N*d_3*delta_d
df3 = pd.DataFrame(D3)
D_sum_3 = df3.sum()

# 计算4阶矩
D4 = N*d_4*delta_d
df4 = pd.DataFrame(D4)
D_sum_4 = df4.sum()

# 求μ
u = (3*D_sum_4*D_sum_2-4*D_sum_3*D_sum_3)/(D_sum_3*D_sum_3-D_sum_4*D_sum_2)

# 计算雨水含量W
pi = math.pi
p = 1
W = (pi*p*D_sum_3)/(6000)

# 计算雨滴的质量加权平均直径Dm
Dm = D_sum_4/D_sum_3

# 计算斜率参数Lambda
Lam = (D_sum_3*(u+4))/D_sum_4

# 计算截距N0
N0 = (D_sum_2*pow(Lam,u+3))/gamma(u+3)

# 输出形状因子u，截距N0和斜率参数L
# print(u,N0,Lam)

# 建立雨滴谱gamma分布
# 建立x的值，生成从0.3到8的值，步长为0.01
x = np.arange(0.3,8,0.01)
# 得把u,N0,Lam变成array才能进行加减乘除计算
x1 = np.array(x)
La = np.array(Lam)
N_0 = np.array(N0)
u1 = np.array(u)
eLx = np.exp(-La*x1)
x1_u = np.power(x1,u1)
r = np.multiply(x1_u,eLx)
y = np.multiply(N_0,r)

# 绘制折线图
plt.plot(x1, y, ls="-", alpha=0.5, linewidth=1, label='abc')
# 设置X轴刻度间隔和标签旋转角度
plt.xlim(0,9)

plt.legend
plt.xlabel('直径')
plt.ylabel('数浓度')

plt.show()