发射度的代码

import numpy as np
from numpy import *
import h5py
import os
import matplotlib.pyplot as plt
import matplotlib.pyplot as plt
import scipy.signal
import time

start = time.time()
ls = []
files = os.listdir(r'D:Raw7')  # 返回路径下的文件和文件名列表

plt.rcParams[r'xtick.direction'] = 'in'  # 设置坐标轴的刻度向内
plt.rcParams['ytick.direction'] = 'in'

for i in files:

    fileName = r'D:Raw7/' + str(i)
    f = h5py.File(fileName, 'r')
    dataname = '/x1'
    dataset = f[dataname]
    x1 = dataset[...]
    # print(mean(x1))
    dataname = '/x2'
    dataset = f[dataname]
    x2 = dataset[...]
    # print(mean(x2))
    dataname = '/x3'
    dataset = f[dataname]
    x3 = dataset[...]

    dataname = '/p1'
    dataset = f[dataname]
    p1 = dataset[...]

    dataname = '/p2'
    dataset = f[dataname]
    p2 = dataset[...]
    b = np.arange(len(x1))

    m = []
    mm = []
    n = []
    nn = []
    a = []
    for j in b:
        if 100 >= x3[j] >= -100:
            m.append(x2[j])
            n.append(p2[j])

    a = m - mean(m)
    z = n - mean(n)
    b = mean(a * a)
    c = mean(z * z)
    e = mean(a * z)
    f = sqrt(b * c - e ** 2)
    ls.append(f / 1.8834)

print(mean(ls))
print(len(ls))
# print(len(n))


# x = np.arange(0, 3.99, 0.0266)
x = np.arange(0, 40.166, 0.266)
# x = np.arange(151)
# plt.ylim(0, 7)
# plt.xlim(0, 60)

plt.xlabel("Propagation Distance(cm)")
plt.ylabel("Proj. Emittance, εN/εN,0")
plt.plot(x, ls, '-', linewidth=1)
end = time.time()
print('运行时间:%s 秒' % (end - start))
plt.show()

三个片粒子：sigama_v       5.35, 5.45        5.589,5.689     5.828,5.928