python 雨滴谱gamma分布飞机作业前中后变化情况

数据：

目标做雨滴谱gamma分布，作业中、作业后1h、作业后2h

代码如下：

#!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('数据.xlsx')

N_ing = df1['作业中']
N_aft1 = df1['作业后1h']
N_aft2 = df1['作业后2h']
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_ing = N_ing*d_2*delta_d
D2_aft1 = N_aft1*d_2*delta_d
D2_aft2 = N_aft2*d_2*delta_d

df2_ing = pd.DataFrame(D2_ing)
df2_aft1 = pd.DataFrame(D2_aft1)
df2_aft2 = pd.DataFrame(D2_aft2)

D_sum_2_ing = df2_ing.sum()
D_sum_2_aft1 = df2_aft1.sum()
D_sum_2_aft2 = df2_aft2.sum()

# 计算3阶矩
D3_ing = N_ing*d_3*delta_d
D3_aft1 = N_aft1*d_3*delta_d
D3_aft2 = N_aft2*d_3*delta_d

df3_ing = pd.DataFrame(D3_ing)
df3_aft1 = pd.DataFrame(D3_aft1)
df3_aft2 = pd.DataFrame(D3_aft2)

D_sum_3_ing = df3_ing.sum()
D_sum_3_aft1 = df3_aft1.sum()
D_sum_3_aft2 = df3_aft2.sum()

# 计算4阶矩
D4_ing = N_ing*d_4*delta_d
D4_aft1 = N_aft1*d_4*delta_d
D4_aft2 = N_aft2*d_4*delta_d

df4_ing = pd.DataFrame(D4_ing)
df4_aft1 = pd.DataFrame(D4_aft1)
df4_aft2 = pd.DataFrame(D4_aft2)

D_sum_4_ing = df4_ing.sum()
D_sum_4_aft1 = df4_aft1.sum()
D_sum_4_aft2 = df4_aft2.sum()

# 求μ
u_ing = (3*D_sum_4_ing*D_sum_2_ing-4*D_sum_3_ing*D_sum_3_ing)/(D_sum_3_ing*D_sum_3_ing-D_sum_4_ing*D_sum_2_ing)
u_aft1 = (3*D_sum_4_aft1*D_sum_2_aft1-4*D_sum_3_aft1*D_sum_3_aft1)/(D_sum_3_aft1*D_sum_3_aft1-D_sum_4_aft1*D_sum_2_aft1)
u_aft2 = (3*D_sum_4_aft2*D_sum_2_aft2-4*D_sum_3_aft2*D_sum_3_aft2)/(D_sum_3_aft2*D_sum_3_aft2-D_sum_4_aft2*D_sum_2_aft2)

# 计算雨水含量W
pi = math.pi
p = 1
W_ing = (pi*p*D_sum_3_ing)/(6000)
W_aft1 = (pi*p*D_sum_3_aft1)/(6000)
W_aft2 = (pi*p*D_sum_3_aft2)/(6000)


# 计算雨滴的质量加权平均直径Dm
Dm_ing = D_sum_4_ing/D_sum_3_ing
Dm_aft1 = D_sum_4_aft1/D_sum_3_aft1
Dm_aft2 = D_sum_4_aft2/D_sum_3_aft2


# 计算斜率参数Lambda
Lam_ing = (D_sum_3_ing*(u_ing+4))/D_sum_4_ing
Lam_aft1 = (D_sum_3_aft1*(u_aft1+4))/D_sum_4_aft1
Lam_aft2 = (D_sum_3_aft2*(u_aft2+4))/D_sum_4_aft2


# 计算截距N0
N0_ing = (D_sum_2_ing*pow(Lam_ing,u_ing+3))/gamma(u_ing+3)
N0_aft1 = (D_sum_2_aft1*pow(Lam_aft1,u_aft1+3))/gamma(u_aft1+3)
N0_aft2 = (D_sum_2_aft2*pow(Lam_aft2,u_aft2+3))/gamma(u_aft2+3)


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

# 建立雨滴谱gamma分布
# 建立x的值，生成从0.3到8的值，步长为0.01
x = np.arange(0.3,8,0.01)
# 得把u,N0,Lam变成array才能进行加减乘除计算
x1_ing = np.array(x)
x1_aft1 = np.array(x)
x1_aft2 = np.array(x)

La_ing = np.array(Lam_ing)
La_aft1 = np.array(Lam_aft1)
La_aft2 = np.array(Lam_aft2)

N_0_ing = np.array(N0_ing)
N_0_aft1 = np.array(N0_aft1)
N_0_aft2 = np.array(N0_aft2)

u1_ing = np.array(u_ing)
u1_aft1 = np.array(u_aft1)
u1_aft2 = np.array(u_aft2)

eLx_ing = np.exp(-La_ing*x1_ing)
eLx_aft1 = np.exp(-La_aft1*x1_aft1)
eLx_aft2 = np.exp(-La_aft2*x1_aft2)

x1_u_ing = np.power(x1_ing,u1_ing)
x1_u_aft1 = np.power(x1_aft1,u1_aft1)
x1_u_aft2 = np.power(x1_aft2,u1_aft2)

r_ing = np.multiply(x1_u_ing,eLx_ing)
r_aft1 = np.multiply(x1_u_aft1,eLx_aft1)
r_aft2 = np.multiply(x1_u_aft2,eLx_aft2)

y_ing = np.multiply(N_0_ing,r_ing)
y_aft1 = np.multiply(N_0_aft1,r_aft1)
y_aft2 = np.multiply(N_0_aft2,r_aft2)
# print(y_aft2)

# 绘制折线图
plt.plot(x1_ing, y_ing, ls="-", alpha=0.5, linewidth=1, label='abc')
plt.plot(x1_aft1, y_aft1, ls="-", alpha=0.5, linewidth=1, label='abc')
plt.plot(x1_aft2, y_aft2, ls="-", alpha=0.5, linewidth=1, label='abc')

# 设置X轴刻度间隔和标签旋转角度
plt.xlim(0,9)

plt.legend
plt.xlabel('直径')
plt.ylabel('数浓度')
plt.legend(['作业中','作业后1h','作业后2h'])

plt.show()