雨滴谱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()
