Python函数汇总

Posted on | In

简介

汇总自行创建的Python函数，主要关于MINEOS和SPECNM相关输入、输出文件的创建与读取。

汇总函数包位于最后一章：函数包中。

模型创建相关

使用MINEOS模型创建相关的nd模型与平均化模型。

使用Obspy.taup创建npz模型以供Taup绘图使用

1
2
3
4
5
6
7
8
9
10
11
12
def create_taupnd_model(model_dir, model_name, output_dir):
    '''
    :param model_dir: nd模型所在的绝对路径，例如：'/home/shizhiyuanubuntu/PycharmProjects/Seismology/models/'。注意：nd模型必须存在 `mantle` `outer-core` `inner-core` 标识。
    :param model_name: 模型名称，不带后缀
    :param output_dir: ndk文件输出绝对路径
    :return: 返回taup model对象，model具有get_travel_times等属性
    '''
    model_dir = model_dir + model_name + '.nd'
    taup_create.build_taup_model(model_dir,
                                 output_folder=output_dir)
    model = TauPyModel(model_name + '.npz')
    return model

通过已有模型进行平均化

需要个性化的内容：

1、原始模型文件的路径。该模型必须为MINEOS所需格式。路径名包含绝对路径及文件名称。例如：
/home/shizhiyuanubuntu/PycharmProjects/Seismology/models/old_models/prem_noocean.txt

2、需要平均化的边界R值列表。
该列表中的所有值为平均化之后各个层的边界半径值，单位为m。注意：需要从小到大排列，且没有重复，并包含地表与球心半径值。
例如：

1
boundaries = [0., 1221500, 3480000, 6291000, 6371000]

3、输出模型名称。
不包含后缀的名称。

1
model_name = 'prem_mean_LVZ_350'

4、ICB与CMB的半径R值。
需要通过该值查找输出响应索引，因此需要指明ICB与CMB的R值。
例如：

1
2
ICB = 1221500.
CMB = 3480000.

5、是否需要添加低速带。
此选项在平均化模型之后进行，手动设置低速带顶底边界和vp、vs速度。
例如：

1
2
3
4
5
ifLVZ = True
boundary1=3480000
boundary2=3480000 + 350000
LVZvp = 11785.04 - 0.075*11785.04
LVZvs = 6414.27 - 0.15*6414.27

代码：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
from numpy import pi
'''
创建平均模型，并且可以自由改变添加低速带。
'''
def find_indices_of_values(lst, value):
    return [index for index, item in enumerate(lst) if item == value]

def rho_integer(bound1, bound2, radius, rho):
    '''
    计算一个区域内的密度平均值
    :param bound1: 底界面半径，单位m，即R。
    :param bound2: 顶界面半径，同上
    :param radius: 整个模型的半径参数
    :param rho: 整个模型的密度参数
    :return: 标量 密度平均值
    '''
    for i in range(0, len(radius)):
        if radius[i] == bound1 and radius[i+1] != bound1: # 寻找bound1边界, 如果为第一个界面，则进行下述步骤：
            index_bound1 = i  # 记录bound1边界的索引。
            rho_sig = 0.
            for j in range(index_bound1, len(radius)):
                if radius[j] != bound2:
                    rho_sig += 4 * pi * (pow(radius[j + 1], 3) - pow(radius[j], 3)) * rho[j] / 3
                elif radius[j] == bound2 and radius[j-1] != bound2:
                    index_bound2 = j
                    rho_mean = (rho_sig / (4 * pi * (pow(radius[index_bound2], 3) - pow(radius[index_bound1], 3)) / 3))
                    return rho_mean
def v_integer(bound1, bound2, radius, v):
    '''
    计算一个区域内的v、q、eta平均值
    :param bound1: 底界面半径，单位m，即R。
    :param bound2: 顶界面半径，同上
    :param radius: 整个模型的半径参数
    :param rho: 整个模型的某个v q eta参数
    :return: 标量 相应平均值
    '''
    for i in range(0, len(radius)):
        if radius[i] == bound1 and radius[i+1] != bound1: # 寻找bound1边界, 如果为第一个界面，则进行下述步骤：
            index_bound1 = i  # 记录bound1边界的索引。
            v_sig = 0.
            for j in range(index_bound1, len(radius)):
                if radius[j] != bound2:
                    v_sig += (radius[j + 1] - radius[j]) * v[j]
                elif radius[j] == bound2 and radius[j-1] != bound2:
                    index_bound2 = j
                    v_mean = (v_sig / (radius[index_bound2] - radius[index_bound1]))
                    return v_mean

def read_mineos_model(modelfile):
    with open(modelfile, 'r') as file:
        lines = file.readlines()
        Radius = []
        rho = []
        vp = []
        vs = []
        qp = []
        qs = []
        eta = []

        for i in range(0, len(lines)):
            Heads = []
            if i == 0 or i==1 or i==2:
                Heads.append(lines[i])
                continue
            elif lines[i] == '\n':
                continue
            Radius.append(float(lines[i][0:8]))
            rho.append(float(lines[i][8:17]))
            vp.append(float(lines[i][17:26]))
            vs.append(float(lines[i][26:35]))
            qp.append(float(lines[i][35:44]))
            qs.append(float(lines[i][44:53]))
            eta.append(float(lines[i][73:80]))
        file.close()
    return Heads, Radius, rho, vp, vs, qp, qs, eta

def enlarge_models(Radius, rho, vp, vs, qp, qs, eta, enlarge_N=300):
    step = int(max(Radius)/enlarge_N)
    Radius_enlarged = []
    rho_enlarged = []
    vp_enlarged = []
    vs_enlarged = []
    qp_enlarged = []
    qs_enlarged = []
    eta_enlarged = []

    for i in range(0, len(Radius)-1):
        if Radius[i] != Radius[i+1]:
            radius = Radius[i]
            for j in range(0, int((Radius[i+1] - Radius[i]) / step)):
                if radius >= Radius[i+1]:
                    break
                Radius_enlarged.append(radius)
                rho_enlarged.append(rho[i])
                vp_enlarged.append(vp[i])
                vs_enlarged.append(vs[i])
                qp_enlarged.append(qp[i])
                qs_enlarged.append(qs[i])
                eta_enlarged.append(eta[i])
                radius += step
            Radius_enlarged.append(Radius[i+1])
            rho_enlarged.append(rho[i])
            vp_enlarged.append(vp[i])
            vs_enlarged.append(vs[i])
            qp_enlarged.append(qp[i])
            qs_enlarged.append(qs[i])
            eta_enlarged.append(eta[i])

    return Radius_enlarged, rho_enlarged, vp_enlarged, vs_enlarged, qp_enlarged, qs_enlarged, eta_enlarged

def Change_v(boundary1, boundary2, changedvp, changedvs, Radius, rho, vp, vs, qp, qs, eta):
    if Radius[0] == boundary1:
        index_bound1 = 0
    else:
        for i in range(1, len(Radius)-1):
            if Radius[i] == boundary1 and Radius[i+1] == boundary1 : # 寻找bound1边界, 如果为第一个界面，则进行下述步骤：
                index_bound1 = i+1
                continue  # 已经存在这个边界，不需要额外添加。
            elif Radius[i] < boundary1 and Radius[i+1] > boundary1 :
                # 不存在边界，添加一个分界面。
                Radius = Radius[:i+1] + [boundary1, boundary1] + Radius[i+1:]
                rho = rho[:i+1] + [rho[i], rho[i+1]] + rho[i+1:]
                vp = vp[:i+1] + [vp[i], vp[i+1]] + vp[i+1:]
                vs = vs[:i+1] + [vs[i], vs[i+1]] + vs[i+1:]
                qp = qp[:i+1] + [qp[i], qp[i+1]] + qp[i+1:]
                qs = qs[:i+1] + [qs[i], qs[i+1]] + qs[i+1:]
                eta = eta[:i+1] + [eta[i], eta[i+1]] + eta[i+1:]
                index_bound1 = i + 2
            elif Radius[i-1] < boundary1 and Radius[i + 1] > boundary1 and Radius[i] == boundary1:
                # 不存在边界，添加一个分界面。
                Radius = Radius[:i + 1] + [boundary1] + Radius[i + 1:]
                rho = rho[:i + 1] + [rho[i+1]] + rho[i + 1:]
                vp = vp[:i + 1] + [vp[i+1]] + vp[i + 1:]
                vs = vs[:i + 1] + [vs[i+1]] + vs[i + 1:]
                qp = qp[:i + 1] + [qp[i+1]] + qp[i + 1:]
                qs = qs[:i + 1] + [qs[i+1]] + qs[i + 1:]
                eta = eta[:i + 1] + [eta[i+1]] + eta[i + 1:]
                index_bound1 = i + 1

    if Radius[len(Radius)-1] == boundary2:
        index_bound2 = len(Radius)-1
    else:
        for i in range(1, len(Radius) - 1):
            if Radius[i] == boundary2 and Radius[i + 1] == boundary2:  # 寻找bound1边界, 如果为第一个界面，则进行下述步骤：
                index_bound2 = i
                continue  # 已经存在这个边界，不需要额外添加。
            elif Radius[i] < boundary2 and Radius[i + 1] > boundary2:
                # 不存在边界，添加一个分界面。
                Radius = Radius[:i + 1] + [boundary2, boundary2] + Radius[i + 1:]
                rho = rho[:i + 1] + [rho[i], rho[i + 1]] + rho[i + 1:]
                vp = vp[:i + 1] + [vp[i], vp[i + 1]] + vp[i + 1:]
                vs = vs[:i + 1] + [vs[i], vs[i + 1]] + vs[i + 1:]
                qp = qp[:i + 1] + [qp[i], qp[i + 1]] + qp[i + 1:]
                qs = qs[:i + 1] + [qs[i], qs[i + 1]] + qs[i + 1:]
                eta = eta[:i + 1] + [eta[i], eta[i + 1]] + eta[i + 1:]
                index_bound2 = i + 1
            elif Radius[i - 1] < boundary2 and Radius[i + 1] > boundary2 and Radius[i] == boundary2:
                # 不存在边界，添加一个分界面。
                Radius = Radius[:i + 1] + [boundary2] + Radius[i + 1:]
                rho = rho[:i + 1] + [rho[i + 1]] + rho[i + 1:]
                vp = vp[:i + 1] + [vp[i + 1]] + vp[i + 1:]
                vs = vs[:i + 1] + [vs[i + 1]] + vs[i + 1:]
                qp = qp[:i + 1] + [qp[i + 1]] + qp[i + 1:]
                qs = qs[:i + 1] + [qs[i + 1]] + qs[i + 1:]
                eta = eta[:i + 1] + [eta[i + 1]] + eta[i + 1:]
                index_bound2 = i

    for j in range(index_bound1, index_bound2 + 1):
        vp[j] = changedvp
        vs[j] = changedvs
    return Radius, rho, vp, vs, qp, qs, eta


# 原始模型文件
modelfile = '/home/shizhiyuanubuntu/PycharmProjects/Seismology/models/old_models/prem_noocean.txt'
# 读取原始模型
Heads, Radius, rho, vp, vs, qp, qs, eta = read_mineos_model(modelfile)
# 需要平均化的边界R值，单位m, 需要从0开始，Rmax结束
boundaries = [0., 1221500, 3480000, 6291000, 6371000]
# 输出模型名称
model_name = 'prem_mean_LVZ_350'
# ICB与CMB边界半径R值
ICB = 1221500.
CMB = 3480000.

ifLVZ = True
boundary1=3480000
boundary2=3480000 + 350000
# LVZvp = 12510.70 - 0.075*12510.70
LVZvp = 11785.04 - 0.075*11785.04
# LVZvs = 6806.58 - 0.15*6806.58
LVZvs = 6414.27 - 0.15*6414.27

# 存储初始模型平均值的列表
Radius_mean1 = []
rho_mean1 = []
vp_mean1 = []
vs_mean1 = []
qp_mean1 = []
qs_mean1 = []
eta_mean1 = []

# 计算初始模型平均值
for p in range(0, len(boundaries)):
    if p != 0 and p != len(boundaries)-1:
        rho_mean = rho_integer(boundaries[p], boundaries[p+1], Radius, rho)
        vp_mean = v_integer(boundaries[p], boundaries[p+1], Radius, vp)
        vs_mean = v_integer(boundaries[p], boundaries[p + 1], Radius, vs)
        qp_mean = v_integer(boundaries[p], boundaries[p + 1], Radius, qp)
        qs_mean = v_integer(boundaries[p], boundaries[p + 1], Radius, qs)
        eta_mean = v_integer(boundaries[p], boundaries[p + 1], Radius, eta)

        Radius_mean1.append(boundaries[p])
        Radius_mean1.append(boundaries[p])
        rho_mean1.append(rho_mean1[len(rho_mean1)-1])
        rho_mean1.append(rho_mean)
        vp_mean1.append(vp_mean1[len(vp_mean1)-1])
        vp_mean1.append(vp_mean)
        vs_mean1.append(vs_mean1[len(vs_mean1)-1])
        vs_mean1.append(vs_mean)
        qp_mean1.append(qp_mean1[len(qp_mean1)-1])
        qp_mean1.append(qp_mean)
        qs_mean1.append(qs_mean1[len(qs_mean1)-1])
        qs_mean1.append(qs_mean)
        eta_mean1.append(eta_mean1[len(eta_mean1)-1])
        eta_mean1.append(eta_mean)

    elif p == 0:
        rho_mean = rho_integer(boundaries[p], boundaries[p+1], Radius, rho)
        vp_mean = v_integer(boundaries[p], boundaries[p+1], Radius, vp)
        vs_mean = v_integer(boundaries[p], boundaries[p + 1], Radius, vs)
        qp_mean = v_integer(boundaries[p], boundaries[p + 1], Radius, qp)
        qs_mean = v_integer(boundaries[p], boundaries[p + 1], Radius, qs)
        eta_mean = v_integer(boundaries[p], boundaries[p + 1], Radius, eta)

        Radius_mean1.append(boundaries[p])
        rho_mean1.append(rho_mean)
        vp_mean1.append(vp_mean)
        vs_mean1.append(vs_mean)
        qp_mean1.append(qp_mean)
        qs_mean1.append(qs_mean)
        eta_mean1.append(eta_mean)

    elif p==len(boundaries)-1:
        Radius_mean1.append(boundaries[p])
        rho_mean1.append(rho_mean1[len(rho_mean1)-1])
        vp_mean1.append(vp_mean1[len(vp_mean1)-1])
        vs_mean1.append(vs_mean1[len(vs_mean1)-1])
        qp_mean1.append(qp_mean1[len(qp_mean1)-1])
        qs_mean1.append(qs_mean1[len(qs_mean1)-1])
        eta_mean1.append(eta_mean1[len(eta_mean1)-1])


# 扩充模型
Radius_enlarged, rho_enlarged, vp_enlarged, vs_enlarged, qp_enlarged, qs_enlarged, eta_enlarged = (
    enlarge_models(Radius=Radius_mean1, rho=rho_mean1, vp=vp_mean1, vs=vs_mean1,
                   qp=qp_mean1, qs=qs_mean1, eta=eta_mean1, enlarge_N=300))
# 增加低速带
if ifLVZ == True:
    Radius_enlarged, rho_enlarged, vp_enlarged, vs_enlarged, qp_enlarged, qs_enlarged, eta_enlarged = Change_v(boundary1,  boundary2,
             changedvp=LVZvp, changedvs=LVZvs,
             Radius=Radius_enlarged, rho=rho_enlarged, vp=vp_enlarged, vs=vs_enlarged,
             qp=qp_enlarged, qs=qs_enlarged, eta=eta_enlarged)

# 平均化的模型名称
num = len(Radius_enlarged)
index_ICB = find_indices_of_values(Radius_enlarged, ICB)
index_CMB = find_indices_of_values(Radius_enlarged, CMB)
with open(model_name + '.txt', 'w') as file:
    file.write(f'{model_name}   THE MEAN MODE\n')
    file.write('   1    1.0000   1\n')
    file.write('{:>6d}{:>4d}{:>4d}\n'.format(num, index_ICB[0]+1, index_CMB[0]+1))
    for i in range(0, len(Radius_enlarged)):
        file.write(
            '{:>7.0f}.{:>9.2f}{:>9.2f}{:>9.2f}{:>9.1f}{:>9.1f}{:>9.2f}{:>9.2f}{:>9.5f}'.format(Radius_enlarged[i],
                                                                                                       rho_enlarged[i],
                                                                                                       vp_enlarged[i],
                                                                                                       vs_enlarged[i],
                                                                                                       qp_enlarged[i],
                                                                                                       qs_enlarged[i],
                                                                                                       vp_enlarged[i],
                                                                                                       vs_enlarged[i],
                                                                                                     eta_enlarged[i]))
        file.write('\n')
    file.write('\n')
    file.close()

读取MINEOS模型的相关内容

简介：读取MINEOS模型

输入：

  • modelfile MINEOS模型文件绝对路径

输出：

  • 一个列表，包含：
    [Radius, rho, vp, vs, qp, qs, eta]
    其中每一个元素为一个列表，如Radius为一个列表，其他元素对应Radius相同索引所在的值。
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
import os

def get_filename(path):
    if os.path.exists(path):
        # 文件存在，可以进行操作，例如打开读取
        ifexist = True
        return ifexist
    else:
        # 文件不存在，可以选择跳过或者处理错误
        ifexist = False
        return ifexist

def MINEOS_model_read(modelfile):
    '''
    :param modelfile: Absolute file path of model
    :return: A 2D list have:
    [Radius, rho, vp, vs, qp, qs, eta]
    '''
    if get_filename(modelfile) == False:
        print(f'NO such file of {modelfile}')
        return []
    with open(modelfile, 'r') as file:
        lines = file.readlines()
        Radius = []
        rho = []
        vp = []
        vs = []
        qp = []
        qs = []
        eta = []

        for i in range(0, len(lines)):
            if i == 0 or i==1 or i==2 or lines[i]=='\n':
                continue
            Radius.append(float(lines[i][0:8]))
            rho.append(float(lines[i][8:17]))
            vp.append(float(lines[i][17:26]))
            vs.append(float(lines[i][26:35]))
            qp.append(float(lines[i][35:44]))
            qs.append(float(lines[i][44:53]))
            eta.append(float(lines[i][73:80]))
        file.close()
    out = []
    out.append(Radius)
    out.append(rho)
    out.append(vp)
    out.append(vs)
    out.append(qp)
    out.append(qs)
    out.append(eta)
    return out

将MINEOS模型转化为nd模型

需要引入读取MINEOS模型的相关内容中的两个函数

个性化修改项：

  • model_name 输出的nd模型的名称，不带后缀。程序会自动添加nd后缀
  • model_dir 输入的MINEOS模型的绝对路径。
    这里的示例中使用相同的模型名称，因此直接调用更加方便。
    第一项为模型所在文件夹路径，第二项为模型名称，第三项为MINEOS模型文件格式后缀。
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
model_name = 'prem_mean_LVZ_300'
model_dir = '/home/shizhiyuanubuntu/PycharmProjects/Seismology/models/' + model_name + '.txt'
out = MINEOS_model_read(model_dir)

with open(model_name + '.nd', 'w') as file:
    Depth = np.flip((6371000. - np.array(out[0])) / 1000.)
    rho = np.flip(np.array(out[1]) / 1000.)
    vp = np.flip(np.array(out[2]) / 1000.)
    vs = np.flip(np.array(out[3]) / 1000.)
    qp = np.flip(np.array(out[4]))
    qs = np.flip(np.array(out[5]))
    for i in range(0, len(Depth)):
        file.write(
            '{:>7.1f}{:>8.1f}{:>8.1f}{:>9.2f}{:>9.1f}{:>9.1f}'.format(Depth[i],
                                                                        vp[i],
                                                                       vs[i],
                                                                       rho[i],
                                                                       qp[i],
                                                                       qs[i]))
        file.write('\n')
    file.close()

绘制MINEOS模型

模型绘制可以自行进行，这里提供一个绘图示例，每一种参数单独绘制在子图中。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
import matplotlib.pyplot as plt
import numpy as np
model_name1 = '/home/shizhiyuanubuntu/PycharmProjects/Seismology/models/prem_mean_noLVZ.txt'
model_name2 = '/home/shizhiyuanubuntu/PycharmProjects/Seismology/models/prem_mean_LVZ_100.txt'
model_name3 = '/home/shizhiyuanubuntu/PycharmProjects/Seismology/models/prem_mean_LVZ_150.txt'
model_name4 = '/home/shizhiyuanubuntu/PycharmProjects/Seismology/models/prem_mean_LVZ_200.txt'
model_name5 = '/home/shizhiyuanubuntu/PycharmProjects/Seismology/models/prem_mean_LVZ_300.txt'
models = [model_name1, model_name2, model_name3, model_name4, model_name5]
plt.rcParams.update({'font.size':20})
layout = [
    ["A", "B", "C", "."],
]
fig, axs = plt.subplot_mosaic(layout, figsize=[18,9])
fig.subplots_adjust(wspace=0.4, hspace=0.4)
fig.suptitle('Earth models for different thickness LVZ above CMB')
axs['A'].set_xlabel('v$_{p}$(km/s)')
axs['A'].set_ylabel('Radius(km)')
axs['B'].set_xlabel('v$_{s}$(km/s)')
axs['B'].set_ylabel('Radius(km)')
axs['C'].set_xlabel('Denisity(g/cm$^3$)')
axs['C'].set_ylabel('Radius(km)')
for n in range(0, len(models)):
    out = MINEOS_model_read(models[n])
    axs["A"].plot(np.array(out[2]) / 1000, np.array(out[0]) / 1000)  # vp
for n in range(0, len(models)):
    out = MINEOS_model_read(models[n])
    axs["B"].plot(np.array(out[3]) / 1000, np.array(out[0]) / 1000 )  # vs
for n in range(0, len(models)):
    out = MINEOS_model_read(models[n])
    axs["C"].plot(np.array(out[1]) / 1000, np.array(out[0]) / 1000 )  # rho

axs['A'].set_yticks([0, 1000, 2000, 3000, 4000, 5000, 6000, 6371])
axs['B'].set_yticks([0, 1000, 2000, 3000, 4000, 5000, 6000, 6371])
axs['C'].set_yticks([0, 1000, 2000, 3000, 4000, 5000, 6000, 6371])
axs['A'].set_xlim([0, 14])
axs['A'].set_ylim([0, 6371])
axs['B'].set_ylim([0, 6371])
axs['C'].set_ylim([0, 6371])

plt.legend(['noLVZ', '100km LVZ', '150km LVZ', '200km LVZ', '300km LVZ'], loc='lower left',  bbox_to_anchor=(1, 0.5))
plt.show()

MINEOS输出内容读取

振型文件读取

振型文件指由MINEOS min_bran输出的文件，一般命名为模型名称+_振型类型，振型类型为S、T或R。

函数

put_freq_in_same：将MINEOS输出的振型文件的一个参数整合为一个列表，按照n、l排列

read_mineos_modefile：读取MINEOS输出的振型文件中的相关信息并输出

put_0_velocity：取同一个径向阶数值的相关参数。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
import matplotlib.pyplot as plt
import obspy
import numpy as np
from obspy import taup

def put_freq_in_same(n, l, freq, nmax=10000, lmax=101):
    '''
    将MINEOS输出的振型文件的一个参数整合为一个列表，按照n、l排列
    :param n: list 列表 读取输出的径向阶数列表n
    :param l: list 列表  读取输出的角度阶数列表l
    :param freq: list 列表  读取输出的频率列表，输入的freq列表为其他项也可以，例如群速度或者相速度
    :param nmax: int 整数  读取的输出径向阶数列表的最大长度，可以尽量设置大一些以保证所有内容存在
    :param lmax: int 整数  读取的输出角度阶数列表的最大长度，可以尽量设置大一些以保证所有内容存在
    :return: list 列表  freqs_all 索引值为[n, l]的频率列表，如freqs_all[1,2]为n=1,l=2振型的本征频率
    '''
    freqs_all = np.zeros([nmax, lmax])
    for i in range(0, len(n)):
        freqs_all[int(n[i]), int(l[i])] = freq[i]
    freqs_all[freqs_all == 0] = np.nan
    return freqs_all

def read_mineos_modefile(modefile_dir, model_name, mode_type):
    '''
    读取MINEOS输出的振型文件中的相关信息并输出。
    :param modefile_dir: string 字符串 读取的振型文件所在文件夹绝对路径，
        例如：'/home/shizhiyuanubuntu/Documents/mineos-1.0.2/DEMO/prem_sim_nolvz',
        此路径末尾不包含'/'，振型文件就在该文件夹下
    :param model_name: string 字符串 MINEOS计算使用的模型名称，例如：prem_mean_noLVZ
    :param mode_type: string 字符串 振型类型，如S、T、R。表征球型、环型、径向振型的振型输出文件
    :return: n, l, frequency, phase_velocity, group_velocity
            输出为五个列表(list)，分别为 径向阶数、角度阶数、频率、相速度、群速度
    '''
    modefile = model_name + '_' + mode_type
    modefile_alldir = modefile_dir + '/' + modefile
    print(modefile_alldir)
    with open(modefile_alldir, 'r') as file:
        lines = file.readlines()
        n = []
        l = []
        frequency = []
        group_velocity = []
        phase_velocity = []
        for i in range(0, len(lines)):
            if lines[i][6:10] == 'mode':
                index = i+2
        for i in range(index, len(lines)):
            n.append(float(lines[i][0:5]))
            l.append(float(lines[i][8:12]))
            phase_velocity.append(float(lines[i][13:24]))
            frequency.append(float(lines[i][25:40]))
            group_velocity.append(float(lines[i][57:72]))
    return n, l, frequency, phase_velocity, group_velocity

def put_0_velocity(n,frequency, ph_vel, grop_vel, ni = 0):
    '''
    取同一个径向阶数值的相关参数。
    :param n: list 振型文件输出的径向阶数
    :param frequency: list 振型文件输出的频率
    :param ph_vel: list 振型文件输出的相速度
    :param grop_vel: list 振型文件输出的群速度
    :param ni: int 输出的径向阶数的值。
    :return: freq, ph_vel0, grop_vel0
        三个列表(list)，分别对应ni值对应的本征频率、相速度、群速度
    '''
    ph_vel0 = []
    grop_vel0 = []
    freq = []
    for i in range(0, len(n)):
        if n[i] == ni:
            freq.append(frequency[i])
            ph_vel0.append(ph_vel[i])
            grop_vel0.append(grop_vel[i])
    return freq, ph_vel0, grop_vel0

振型文件读取内容绘图示例

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
# 读取振型文件
n_nolvz, l_nolvz, frequency_nolvz, phase_velocity_nolvz, group_velocity_nolvz = (
    read_mineos_modefile(modefile_dir='/home/shizhiyuanubuntu/Documents/mineos-1.0.2/DEMO/prem_sim_nolvz',
                     model_name='prem_mean_noLVZ', mode_type='S'))
lvz_no_freq = put_freq_in_same(n_nolvz, l_nolvz, frequency_nolvz, nmax=40, lmax=500)

n_lvz_100, l_lvz_100, frequency_lvz_100, phase_velocity_lvz_100, group_velocity_lvz_100 = (
    read_mineos_modefile(modefile_dir='/home/shizhiyuanubuntu/Documents/mineos-1.0.2/DEMO/prem_sim_lvz100',
                     model_name='prem_mean_LVZ_100', mode_type='S'))
lvz_100_freq = put_freq_in_same(n_lvz_100, l_lvz_100, frequency_lvz_100, nmax=40, lmax=500)

n_lvz_150, l_lvz_150, frequency_lvz_150, phase_velocity_lvz_150, group_velocity_lvz_150 = (
    read_mineos_modefile(modefile_dir='/home/shizhiyuanubuntu/Documents/mineos-1.0.2/DEMO/prem_sim_lvz150',
                     model_name='prem_mean_LVZ_150', mode_type='S'))
lvz_150_freq = put_freq_in_same(n_lvz_150, l_lvz_150, frequency_lvz_150, nmax=40, lmax=500)

n_lvz_200, l_lvz_200, frequency_lvz_200, phase_velocity_lvz_200, group_velocity_lvz_200 = (
    read_mineos_modefile(modefile_dir='/home/shizhiyuanubuntu/Documents/mineos-1.0.2/DEMO/prem_sim_lvz200',
                     model_name='prem_mean_LVZ_200', mode_type='S'))
lvz_200_freq = put_freq_in_same(n_lvz_200, l_lvz_200, frequency_lvz_200, nmax=40, lmax=500)

n_lvz_300, l_lvz_300, frequency_lvz_300, phase_velocity_lvz_300, group_velocity_lvz_300 = (
    read_mineos_modefile(modefile_dir='/home/shizhiyuanubuntu/Documents/mineos-1.0.2/DEMO/prem_sim_lvz300',
                     model_name='prem_mean_LVZ_300', mode_type='S'))
lvz_300_freq = put_freq_in_same(n_lvz_300, l_lvz_300, frequency_lvz_300, nmax=40, lmax=500)

# 设置颜色
clolor_zise = (0.267, 0.16, 0.353)
clolor_zanglan = (0.255, 0.243, 0.522)
clolor_shenlan = (0.188, 0.408, 0.553)
clolor_lanlv = (0.122, 0.573, 0.545)
clolor_lvse = (0.208, 0.718, 0.467)
clolor_lvhuang = (0.569, 0.835, 0.259)
clolor_huangse = (0.973, 0.902, 0.125)

# 创建fig
layout = [
    ["A", "A", "A", "B", "B"],
    ["A", "A","A", "C", "C"],
    ["A", "A","A", ".", "."],
]
fig, axs = plt.subplot_mosaic(layout, figsize=[18,9])
fig.subplots_adjust(wspace=0.5, hspace=0.6)
fig.suptitle('Rayleigh wave dispersion curve for different LVZ thickness')
# 设置坐标
axs['A'].set_xlabel('Angular order')
axs['A'].set_ylabel('frequency(mHz)')
axs['B'].set_xlabel('Angular order')
axs['B'].set_ylabel('frequency(mHz)')
axs['C'].set_xlabel('Angular order')
axs['C'].set_ylabel('frequency(mHz)')
# 绘制n=0-10、l=0-101范围内的频散曲线
for n in range(0, 10):
    axs['A'].plot([i for i in range(0, 101)], lvz_no_freq[n, 0:101], color = 'k')
    axs['A'].plot([i for i in range(0, 101)], lvz_100_freq[n, 0:101], color = 'grey')
    axs['A'].plot([i for i in range(0, 101)], lvz_150_freq[n, 0:101], color = 'lightcoral')
    axs['A'].plot([i for i in range(0, 101)], lvz_200_freq[n, 0:101], color = 'brown')
    axs['A'].plot([i for i in range(0, 101)], lvz_300_freq[n, 0:101], color='maroon')
    axs['A'].text(100, np.nanmax(lvz_no_freq[n, 0:101]), '{}'.format(n))
# 绘制n=0、l=0-101范围内的频散曲线
for n in range(0, 1):
    axs['B'].plot([i for i in range(0, 101)], lvz_no_freq[n, 0:101], color = 'k')
    axs['B'].plot([i for i in range(0, 101)], lvz_100_freq[n, 0:101], color = 'grey')
    axs['B'].plot([i for i in range(0, 101)], lvz_150_freq[n, 0:101], color = 'lightcoral')
    axs['B'].plot([i for i in range(0, 101)], lvz_200_freq[n, 0:101], color = 'brown')
    axs['B'].plot([i for i in range(0, 101)], lvz_300_freq[n, 0:101], color='maroon')
    # axs['B'].text(100, np.nanmax(lvz_no_freq[n, 0:101]), '{}'.format(n))
    axs['B'].set_title('Radial order=0', fontsize=14)
# 绘制n=1、l=0-101范围内的频散曲线
for n in range(1, 2):
    axs['C'].plot([i for i in range(0, 101)], lvz_no_freq[n, 0:101], color = 'k')
    axs['C'].plot([i for i in range(0, 101)], lvz_100_freq[n, 0:101], color = 'grey')
    axs['C'].plot([i for i in range(0, 101)], lvz_150_freq[n, 0:101], color = 'lightcoral')
    axs['C'].plot([i for i in range(0, 101)], lvz_200_freq[n, 0:101], color = 'brown')
    axs['C'].plot([i for i in range(0, 101)], lvz_300_freq[n, 0:101], color='maroon')
    # axs['C'].text(100, np.nanmax(lvz_no_freq[n, 0:101]), '{}'.format(n))
    axs['C'].set_title('Radial order=1', fontsize=14)
axs['C'].legend(['noLVZ', '100km LVZ', '150km LVZ', '200km LVZ', '300km LVZ'], loc='lower left',  bbox_to_anchor=(-0., -1.86))

本征函数文件读取

函数

  • mineos_file_read
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
import numpy as np
import matplotlib.pyplot as plt
from plot_single_eigenfunction import plot_W_eigenfunction, plot_U_eigenfunction
from plot_eigenfunctions import open_file_and_plot_eigenfunctions_signal
import os

def get_filename(path):
    if os.path.exists(path):
        # 文件存在，可以进行操作，例如打开读取
        ifexist = True
        return ifexist
    else:
        # 文件不存在，可以选择跳过或者处理错误
        ifexist = False
        return ifexist

def mineos_file_read(MINEOS_eigenfunction_dir, Modetype, radial_order, angular_order, Mineos_file=[]):
    '''
     读取MINEOS单个本征函数文件；
    :param MINEOS_eigenfunction_dir: string MINEOS 文件绝对路径，
        例如：'/home/shizhiyuanubuntu/PycharmProjects/Seismology/MINEOS_eigenfunctions/'
        程序内部根据Modetype自动添加'test_T_ASC/'或'test_S_ASC/'等文件夹名。
    :param Modetype: string 振型类型，S、T、R
    :param radial_order: int 径向阶数
    :param angular_order: int 角度阶数
    :param Mineos_file: Mineos本征函数文件的名称.ASC，可以不指定该名称。
    如果指定则按照改名称选择，忽略radial_oder和angular_order参数；
    如果不指定则需要指定radial_order、angular_order
    :return: 返回值均为list
    若振型类型为S，返回：Radius, U_eigenfunction, V_eigenfunction, P_eigenfunction
    分别表示半径列表、U分量、V分量、P分量（考虑重力相关参数）
    若振型类型为T，返回：Radius, W_eigenfunction
    若振型类型为R，返回：Radius, R_eigenfunction
    '''
    radial_orderstr, angular_orderstr = MINEOS_eigenfunction_file_name(radial_order, angular_order)
    if Mineos_file == []:
        if Modetype == 'T' or Modetype == 't' or Modetype == 'W':
            modename = 'T.' + radial_orderstr + '.' + angular_orderstr + '.ASC'
            print(modename)
            T_MINEOS_eigenfunction_file = MINEOS_eigenfunction_dir + 'test_T_ASC/' + modename
            if get_filename(T_MINEOS_eigenfunction_file) == False:
                print(f'NO such file of {T_MINEOS_eigenfunction_file}')
                return [], []
            with open(T_MINEOS_eigenfunction_file, 'r') as file:
                lines = file.readlines()
                W_eigenfunction = []
                Radius = []
                for i in range(0, len(lines)):
                    if i == 0:
                        continue
                    Radius.append(float(lines[i][0:8]))
                    W_eigenfunction.append(float(lines[i][8:23]))
                file.close()
            return Radius, W_eigenfunction

        elif Modetype == 'S' or Modetype == 's' or Modetype == 'U'  or Modetype == 'V':
            modename = 'S.' + radial_orderstr + '.' + angular_orderstr + '.ASC'
            print(modename)
            S_MINEOS_eigenfunction_file = MINEOS_eigenfunction_dir + 'test_S_ASC/' + modename
            if get_filename(S_MINEOS_eigenfunction_file) == False:
                print(f'NO such file of {S_MINEOS_eigenfunction_file}')
                return [], [], [], []
            with open(S_MINEOS_eigenfunction_file, 'r') as file:
                lines = file.readlines()
                U_eigenfunction = []
                V_eigenfunction = []
                P_eigenfunction = []
                Radius = []
                for i in range(0, len(lines)):
                    if i == 0:
                        continue
                    else:
                        Radius.append(float(lines[i][0:8]))
                        U_eigenfunction.append(float(lines[i][8:23]))
                        V_eigenfunction.append(float(lines[i][38:53]))
                        P_eigenfunction.append(float(lines[i][68:83]))
                file.close()
            return Radius, U_eigenfunction, V_eigenfunction, P_eigenfunction
        elif Modetype == 'R' or Modetype == 'r':
            modename = 'S.' + radial_orderstr + '.' + angular_orderstr + '.ASC'
            print(modename)
            R_MINEOS_eigenfunction_file = MINEOS_eigenfunction_dir + 'test_R_ASC/' + modename
            if get_filename(R_MINEOS_eigenfunction_file) == False:
                print(f'NO such file of {R_MINEOS_eigenfunction_file}')
                return [], []
            with open(R_MINEOS_eigenfunction_file, 'r') as file:
                lines = file.readlines()
                R_eigenfunction = []
                Radius = []
                for i in range(0, len(lines)):
                    if i == 0:
                        continue
                    Radius.append(float(lines[i][0:8]))
                    R_eigenfunction.append(float(lines[i][8:23]))
                file.close()
            return Radius, R_eigenfunction
        else:
            if Modetype == 'T' or Modetype == 't'  or Modetype == 'W':
                T_MINEOS_eigenfunction_file = MINEOS_eigenfunction_dir + 'test_T_ASC/' + Mineos_file
                if get_filename(T_MINEOS_eigenfunction_file) == False:
                    print(f'NO such file of {T_MINEOS_eigenfunction_file}')
                    return [], []
                with open(T_MINEOS_eigenfunction_file, 'r') as file:
                    lines = file.readlines()
                    W_eigenfunction = []
                    Radius = []
                    for i in range(0, len(lines)):
                        if i == 0:
                            continue
                        Radius.append(float(lines[i][0:8]))
                        W_eigenfunction.append(float(lines[i][8:23]))
                    file.close()
                return Radius, W_eigenfunction

            elif Modetype == 'S' or Modetype == 's' or Modetype == 'U' or Modetype == 'V':
                S_MINEOS_eigenfunction_file = MINEOS_eigenfunction_dir + 'test_S_ASC/' + Mineos_file
                if get_filename(S_MINEOS_eigenfunction_file) == False:
                    print(f'NO such file of {S_MINEOS_eigenfunction_file}')
                    return [], [], [], []
                with open(S_MINEOS_eigenfunction_file, 'r') as file:
                    lines = file.readlines()
                    U_eigenfunction = []
                    V_eigenfunction = []
                    P_eigenfunction = []
                    Radius = []
                    for i in range(0, len(lines)):
                        if i == 0:
                            continue
                        else:
                            Radius.append(float(lines[i][0:8]))
                            U_eigenfunction.append(float(lines[i][8:23]))
                            V_eigenfunction.append(float(lines[i][38:53]))
                            P_eigenfunction.append(float(lines[i][68:83]))
                    file.close()
                return Radius, U_eigenfunction, V_eigenfunction, P_eigenfunction
            elif Modetype == 'R' or Modetype == 'r':
                R_MINEOS_eigenfunction_file = MINEOS_eigenfunction_dir + 'test_R_ASC/' + Mineos_file
                if get_filename(R_MINEOS_eigenfunction_file) == False:
                    print(f'NO such file of {R_MINEOS_eigenfunction_file}')
                    return [], []
                with open(R_MINEOS_eigenfunction_file, 'r') as file:
                    lines = file.readlines()
                    R_eigenfunction = []
                    Radius = []
                    for i in range(0, len(lines)):
                        if i == 0:
                            continue
                        Radius.append(float(lines[i][0:8]))
                        R_eigenfunction.append(float(lines[i][8:23]))
                    file.close()
                return Radius, R_eigenfunction

def numofzero(numb):
    '''
    添加0,因为MINEOS输出的本征函数文件命名规则需要添加很多0,因此使用该函数创建本征函数名以供读取使用。
    :param numb: int 输入的数字
    :return: string 增加0之后的数字
    '''
    if numb < 10:
        itsstr = '000000' + str(numb)
    elif numb < 100:
        itsstr = '00000' + str(numb)
    elif numb < 1000:
        itsstr = '0000' + str(numb)
    elif numb < 10000:
        itsstr = '000' + str(numb)
    return itsstr


def MINEOS_eigenfunction_file_name(radial_order, angular_order):
    '''
    创建MINEOS本征函数名
    :param radial_order: int 径向阶数 
    :param angular_order: int 角度阶数
    :return: string
        radial_orderstr: 径向阶数补0后的字符串类型
        angular_orderstr: 径向阶数补0后的字符串类型
    '''
    radial_orderstr = numofzero(radial_order)
    angular_orderstr = numofzero(angular_order)
    return radial_orderstr, angular_orderstr

SPECNM输出内容读取

振型相关文件的读取

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
def read_SPECNM_mode_frequency_output(mode_file_freq):
    '''
    读取SPECNM振型输出文件；仅输出频率信息
    :param mode_file_freq: string 振型相关文件绝对路径，即该文件所在绝对路径及文件名——完整的路径地址
    :return: list n, l, freq, period
        返回 径向阶数、角度阶数、频率、周期 的列表
    '''
    with open(mode_file_freq, 'r') as file:
        lines = file.readlines()
        n = []
        l = []
        freq = []
        period = []
        for i in range(0, len(lines)):
            if i == 0 or i == 1 or lines[i] == '\n':
                continue
            n.append(float(lines[i][0:6]))
            l.append(float(lines[i][6:12]))
            freq.append(float(lines[i][12:27]))
            period.append(float(lines[i][27:39]))
        file.close()
    return n, l, freq, period

def read_SPECNM_mode_output(mode_file, modetype):
    '''
    读取SPECNM振型输出文件；相关函数为 Love_output 等SPECNM振型输出函数
    :param mode_file: string 振型相关文件绝对路径，即该文件所在绝对路径及文件名——完整的路径地址
    例如：'/home/shizhiyuanubuntu/PycharmProjects/Seismology/prem_noocean/Rayleigh_output/prem_noocean_S_output.txt'
    :param modetype: string 振型类型，为S、T、R
    :return: list
    若振型类型为S或T，输出列表 n, l, freq, period, ph_v, grp_v, gama, Q, epsilon
    为径向阶数、角度阶数、频率、周期、相速度、群速度、gama、品质因子Q、epsilon
    '''
    if modetype == 'S' or modetype == 'T':
        with open(mode_file, 'r') as file:
            lines = file.readlines()
            n = []
            l = []
            freq = []
            period = []
            ph_v = []
            grp_v = []
            gama = []
            Q = []
            epsilon = []
            '{:>6d}{:>6d}{:>15.5f}{:>12.3f}{:>12.3f}{:>12.3f}{:>12.3f}{:>12.3f}{:>12.3f}'
            for i in range(0, len(lines)):
                if i == 0 or i == 1 or lines[i] == '\n':
                    continue
                n.append(float(lines[i][0:6]))
                l.append(float(lines[i][6:12]))
                freq.append(float(lines[i][12:27]))
                period.append(float(lines[i][27:39]))
                ph_v.append(float(lines[i][39:51]))
                grp_v.append(float(lines[i][51:63]))
                gama.append(float(lines[i][63:75]))
                Q.append(float(lines[i][75:87]))
                epsilon.append(float(lines[i][87:99]))
            file.close()
        return n, l, freq, period, ph_v, grp_v, gama, Q, epsilon
    elif modetype=='R':
        with open(mode_file, 'r') as file:
            lines = file.readlines()
            n = []
            l = []
            freq = []
            period = []
            gama = []
            Q = []
            epsilon = []
            for i in range(0, len(lines)):
                if i == 0 or i == 1 or lines[i] == '\n':
                    continue
                n.append(float(lines[i][0:6]))
                l.append(float(lines[i][6:12]))
                freq.append(float(lines[i][12:27]))
                period.append(float(lines[i][27:39]))
                gama.append(float(lines[i][39:51]))
                Q.append(float(lines[i][51:63]))
                epsilon.append(float(lines[i][63:75]))
            file.close()
        return n, l, freq, period,[], [], gama, Q, epsilon

def find_freq(mode_file, modetype, angular_order, radial_order):
    '''
    寻找SPECNM特定一系列振型的频率，已知振型参数n、l。
    :param mode_file: string 振型相关文件绝对路径，即该文件所在绝对路径及文件名——完整的路径地址
    :param modetype: 振型类型，为S、T、R
    :param angular_order:  1D array, 角度阶数列表，
    :param radial_order:  1D array, 径向阶数列表
    :return: list 频率列表
    '''
    n, l, freq, period = read_SPECNM_mode_frequency_output(mode_file_freq=mode_file)
    out_freq = []
    for i in range(0, len(angular_order)):
        for j in range(0, len(n)):
            if radial_order[i] == n[j] and angular_order[i] == l[j]:
                out_freq.append(freq[j])
    return out_freq

读取并绘制本征函数

绘制单个本征函数图像

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
from specnm import rayleigh_fg
from matplotlib import pyplot as plt
from specnm import seismogram
import numpy as np

def open_file_and_plot_eigenfunctions_signal(radial_order, angular_order, data_dir, eigenfunction_type='U', linestyle='k-',ifdr = False, ifnormalized=0):
    '''
    输出单个本征函数文件，绘制图像。输出为plot的一条曲线，需要自行设置figure以及plot.show()
    :param radial_order: 径向阶数n
    :param angular_order: 角度阶数l
    :param data_dir: 本征函数文件所在绝对路径，例如：
    该本征函数文件为radord9angord1.txt，其路径位于：
    '/home/shizhiyuanubuntu/PycharmProjects/Seismology/prem_noocean/Rayleigh_output/eigenfunctions/'
    :param eigenfunction_type: 本征函数的类型，有U、V、W三种类型，分别对应S振型的两种分量和T振型一种分量。
    :param linestyle: 颜色和线型，例如'k-'
    :param ifdr: 该文件是否是本征函数的导数（关于radius的导数），默认为False
    :param ifnormalized: 本征函数是否进行归一化 0:不归一化 1：以地表位移为参考 2：以最大值为参考
    :return: 无，仅绘制本征函数曲线，为一个plot函数。需要在函数外部进行figure创建和plt.show
    '''
    if ifdr == True: # 如果绘制的是eigenfunction_dr.
        eigenfunction_file = 'dr_' + eigenfunction_type + '_' + f'radord{radial_order}_angord{angular_order}.txt'
    else:
        eigenfunction_file = eigenfunction_type + '_' + f'radord{radial_order}_angord{angular_order}.txt'
    file = data_dir + eigenfunction_file
    radius_file = data_dir + 'Radius.txt'
    with open(file, 'r') as file:
        lines = file.readlines()
        eigenfunction = []
        eigenfunction_dr = []
        for line in lines:
            eigenfunction.append(float(line))
        file.close()

    with open(radius_file, 'r') as file:
        lines = file.readlines()
        radius = []
        for line in lines:
            radius.append(float(line))
        file.close()

    # print(eigenfunction)
    eigenfunction = np.array(eigenfunction)
    radius = np.array(radius) / 1000
    if ifnormalized == 1:
        plt.plot(radius, eigenfunction / eigenfunction[len(eigenfunction) - 1], linestyle)  # 归一化后
    elif ifnormalized == 2:
        plt.plot(radius, eigenfunction / abs(max(eigenfunction)), linestyle)  # 按照最大值归一化
    else:
        plt.plot(radius, eigenfunction, linestyle)  # 非归一化

绘制本征函数图像

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
import specnm
from matplotlib import pyplot as plt
import numpy as np

def plot_W_eigenfunction(wdir, radial_order, angular_order, ifnormalized=0, 
                         ifnegative=False, mode_type='W', linecolorsytle='k-', 
                         linewidth=1, if_radius_normalized = False, ifaxis=False, 
                         ifreverse=False, iflines=False):
    '''
    绘制特定振型的本征函数
    :param wdir:  本征函数文件所在路径，例如W_radord1angord2.txt位于该路径下：
    '/home/shizhiyuanubuntu/PycharmProjects/Seismology/prem_noocean/Love_output/eigenfunctions/'
    :param radial_order: int 径向阶数
    :param angular_order: int 角度阶数
    :param ifnormalized: 是否进行归一化。
        0：不进行归一化；1：使用地表位移作为基准进行归一化； 2: 使用绝对值最大值作为基准进行归一化 
        默认为0
    :param ifnegative: 是否对本征函数值进行反转，即加一个负号；
    :param mode_type: string 本征函数振型类型，尽管在路径中已经指出不同的振型类型存放于不同的文件夹下，但是依然需要指定
        振型类型以方便创建需要读取的文件的名称。分为U、V、W三种类型；Radial_output输出仍为U、V分量。
    :param linecolorsytle: string 绘制的本征函数的曲线颜色和线型，默认'k-'
    :param linewidth: float 绘制的本征函数的曲线宽度，默认为1
    :param if_radius_normalized: 输出的图像坐标轴中半径是否归一化；归一化后半径范围为相对值0-1。
    :param ifaxis: 是否在地表处显示一条横线表示在边界面；以及是否显示地球660间断面和CMB。
    :param ifreverse: False: 半径为y轴；True：半径为x轴。
    :return: 只进行plot，需要在函数外围设置figure和plt.show()
    '''

    weigenfunction_file = mode_type + '_' + f'radord{radial_order}_angord{angular_order}.txt'
    radial_file = 'Radius.txt'
    weigenfunction_file = wdir + weigenfunction_file
    wradial_file = wdir + radial_file
    with open(weigenfunction_file, 'r') as file:
        lines = file.readlines()
        weigenfunction = []
        for line in lines:
            weigenfunction.append(float(line))
        file.close()
    with open(wradial_file, 'r') as file:
        lines = file.readlines()
        wradius = []
        for line in lines:
            wradius.append(float(line))
        file.close()

    if ifnegative==True:
        if weigenfunction[np.argmax(np.abs(np.array(weigenfunction)))] < 0:
            weigenfunction = -np.array(weigenfunction)

    if ifnormalized == 1:
        if ifreverse == False:
            if if_radius_normalized == False:
                plt.plot(np.array(weigenfunction)/weigenfunction[len(weigenfunction)-1], np.array(wradius)/1000, linecolorsytle, linewidth=linewidth)
            else:
                plt.plot(np.array(weigenfunction) / weigenfunction[len(weigenfunction) - 1], np.array(wradius) / max(np.array(wradius)),
                         linecolorsytle, linewidth=linewidth)
        else:
            if if_radius_normalized == False:
                plt.plot(np.array(wradius) / 1000, np.array(weigenfunction) / weigenfunction[len(weigenfunction) - 1],
                     linecolorsytle, linewidth=linewidth)
            else:
                plt.plot(np.array(wradius) / max(np.array(wradius)), np.array(weigenfunction) / weigenfunction[len(weigenfunction) - 1],
                         linecolorsytle, linewidth=linewidth)

        if ifaxis == 1:
            if weigenfunction[0] < 0:
                plt.axhline(y=wradius[0]/1000, xmin=weigenfunction[0] / weigenfunction[len(weigenfunction)-1], xmax=0, color='k',
                            linestyle='-', linewidth=1, alpha=1)
            else:
                plt.axhline(y=wradius[0] / 1000, xmin=0, xmax=weigenfunction[0] / weigenfunction[len(weigenfunction)-1], color='k',
                            linestyle='-', linewidth=1, alpha=1)
            plt.text(-1, wradius[0]/1000, 'CMB')
            plt.text(-1, max(wradius)/1000-670, '670')
            plt.plot([-1, -0.7], [max(wradius) / 1000 - 670, max(wradius) / 1000 - 670], 'k-')

            if weigenfunction[len(weigenfunction) - 1] < 0:
                plt.axhline(y=wradius[len(wradius) - 1] / 1000, xmin=weigenfunction[len(weigenfunction) - 1] / max(weigenfunction),
                            xmax=0,
                            color='k', linestyle='-', linewidth=1, alpha=1)
            else:
                plt.axhline(y=wradius[len(wradius) - 1] / 1000, xmin=0,
                            xmax=weigenfunction[len(weigenfunction) - 1] / weigenfunction[len(weigenfunction)-1],
                            color='k', linestyle='-', linewidth=1, alpha=1)
        if ifaxis == 1:
            plt.text(-1, wradius[len(wradius)-1]/1000, '0')
    elif ifnormalized == 2:  # 当使用绝对值最大值作为归一化条件时
        if ifreverse == False:
            if if_radius_normalized == False:
                plt.plot(np.array(weigenfunction)/max(abs(np.array(weigenfunction))), np.array(wradius)/1000, linecolorsytle, linewidth=linewidth)
            else:
                plt.plot(np.array(weigenfunction) / max(abs(np.array(weigenfunction))), np.array(wradius) / max(np.array(wradius)),
                         linecolorsytle, linewidth=linewidth)
        else:
            if if_radius_normalized == False:
                plt.plot(np.array(wradius) / 1000, np.array(weigenfunction) / max(abs(np.array(weigenfunction))),
                     linecolorsytle, linewidth=linewidth)
            else:
                plt.plot(np.array(wradius) / max(np.array(wradius)), np.array(weigenfunction) / max(abs(np.array(weigenfunction))),
                         linecolorsytle, linewidth=linewidth)

        if ifaxis == 1:
            if weigenfunction[0] < 0:
                # plt.axhline(y=wradius[0]/1000, xmin=weigenfunction[0]/max(abs(np.array(weigenfunction))), xmax=0, color='k',
                #            linestyle='-', linewidth=1, alpha=1)
                plt.plot([weigenfunction[0] / max(abs(np.array(weigenfunction))), 0],
                         [wradius[0] / 1000, wradius[0] / 1000], 'k-')
            else:
                # plt.axhline(y=wradius[0]/1000, xmin=0, xmax=weigenfunction[0] / max(abs(np.array(weigenfunction))), color='k',
                #            linestyle='-', linewidth=1, alpha=1)
                plt.plot([0, weigenfunction[0] / max(abs(np.array(weigenfunction)))],
                         [wradius[0] / 1000, wradius[0] / 1000], 'k-')
            plt.text(-1, wradius[0]/1000, 'CMB')
            plt.text(-1, max(wradius)/1000 - 670, '670')
            plt.plot([-1,-0.7], [max(wradius) / 1000 - 670, max(wradius) / 1000 - 670], 'k-')

            if weigenfunction[len(weigenfunction) - 1] < 0:
                # plt.axhline(y=wradius[len(wradius)-1]/1000, xmin=weigenfunction[len(weigenfunction) - 1] / max(abs(np.array(weigenfunction))),
                #             xmax=0, color='k', linestyle='-', linewidth=1, alpha=1)
                plt.plot([weigenfunction[len(weigenfunction) - 1] / max(abs(np.array(weigenfunction))), 0],
                         [wradius[len(wradius)-1]/1000, wradius[len(wradius)-1]/1000], 'k-')
            else:
                # plt.axhline(y=wradius[len(wradius)-1]/1000, xmin=0,
                 #            xmax=weigenfunction[len(weigenfunction) - 1] / max(abs(np.array(weigenfunction))),
                 #            color='k', linestyle='-', linewidth=1, alpha=1)\
                plt.plot([0, weigenfunction[len(weigenfunction) - 1]]/ max(abs(np.array(weigenfunction))),
                         [wradius[len(wradius)-1]/1000, wradius[len(wradius)-1]/1000], 'k-')
        if ifaxis == 1:
            plt.text(-1, wradius[len(wradius)-1]/1000, '0')

    else: # 不进行归一化
        if ifreverse == False:
            if if_radius_normalized == False:
                plt.plot(np.array(weigenfunction), np.array(wradius)/1000, linecolorsytle, linewidth=linewidth)
            else:
                plt.plot(np.array(weigenfunction), np.array(wradius) / max(np.array(wradius)), linecolorsytle, linewidth=linewidth)
        else:
            if if_radius_normalized == False:
                plt.plot(np.array(wradius) / 1000, np.array(weigenfunction),
                     linecolorsytle, linewidth=linewidth)
            else:
                plt.plot(np.array(wradius) / max(np.array(wradius)), np.array(weigenfunction),
                         linecolorsytle, linewidth=linewidth)

        if ifaxis == 1:
            plt.text(-1, wradius[0]/1000, 'CMB')
            plt.text(-1, max(wradius)/1000 - 670, '670')
            plt.plot([-1,-0.7], [max(wradius) / 1000 - 670, max(wradius) / 1000 - 670], 'k-')

            if weigenfunction[0] < 0:
                # plt.axhline(y=wradius[0]/1000, xmin=weigenfunction[0]/max(abs(np.array(weigenfunction))), xmax=0, color='k',
                #            linestyle='-', linewidth=1, alpha=1)
                plt.plot([weigenfunction[0] / max(abs(np.array(weigenfunction))), 0],
                         [wradius[0] / 1000, wradius[0] / 1000], 'k-')
            else:
                # plt.axhline(y=wradius[0]/1000, xmin=0, xmax=weigenfunction[0] / max(abs(np.array(weigenfunction))), color='k',
                #            linestyle='-', linewidth=1, alpha=1)
                plt.plot([0, weigenfunction[0] / max(abs(np.array(weigenfunction)))],
                         [wradius[0] / 1000, wradius[0] / 1000], 'k-')

            if weigenfunction[len(weigenfunction) - 1] < 0:
                # plt.axhline(y=wradius[len(wradius)-1]/1000, xmin=weigenfunction[len(weigenfunction) - 1] / max(abs(np.array(weigenfunction))),
                #             xmax=0, color='k', linestyle='-', linewidth=1, alpha=1)
                plt.plot([weigenfunction[len(weigenfunction) - 1] / max(abs(np.array(weigenfunction))), 0],
                         [wradius[len(wradius)-1]/1000, wradius[len(wradius)-1]/1000], 'k-')
            else:
                # plt.axhline(y=wradius[len(wradius)-1]/1000, xmin=0,
                 #            xmax=weigenfunction[len(weigenfunction) - 1] / max(abs(np.array(weigenfunction))),
                 #            color='k', linestyle='-', linewidth=1, alpha=1)\
                plt.plot([0, weigenfunction[len(weigenfunction) - 1]]/ max(abs(np.array(weigenfunction))),
                         [wradius[len(wradius)-1]/1000, wradius[len(wradius)-1]/1000], 'k-')
        if ifaxis == 1:
            plt.text(-1, wradius[len(wradius)-1]/1000, '0')

函数包

SPECNM_read.py

简介：关于SPECNM软件相关输出内容的读取、绘图、和一些数学操作。

函数目录：

  • get_filename(path) 对文件是否存在进行判断与处理
  • read_SPECNM_mode_frequency_output(mode_file_freq) 读取SPECNM振型输出文件；仅输出频率信息
  • find_specnm_mode_freq(mode_file, radial_order, angular_order)
  • SPECNM_mode_name_make(all_dir, model_name, mode_type)
  • read_SPECNM_mode_output(mode_file, modetype) 读取SPECNM振型输出文件；相关函数为 Love_output 等SPECNM振型输出函数
  • find_freq(mode_file, modetype, angular_order, radial_order) 寻找SPECNM特定一系列振型的频率，已知振型参数n(radial order)、l(angular order)
  • open_file_and_plot_eigenfunctions_signal(radial_order, angular_order, data_dir, eigenfunction_type=‘U’, linestyle=‘k-’,ifdr = False, ifnormalized=0) 输出单个本征函数文件，绘制图像。输出为plot的一条曲线，需要自行设置figure以及plot.show()
  • specnm_file_read(wdir, mode_type, radial_order, angular_order) 读取SPECNM本征函数文件
  • plot_specnm_eigenfunction_with_eig(wradius, weigenfunction , ifnormalized=0,
    ifnegative=False, mode_type=‘W’, linecolorsytle=‘k-’,
    linewidth=1, if_radius_normalized = False, ifaxis=False, ifboundary=False,
    ifreverse=False, iflines=False) 绘制特定振型的本征函数，输入为本征函数值
  • plot_W_eigenfunction(wdir, radial_order, angular_order, ifnormalized=0,
    ifnegative=False, mode_type=‘W’, linecolorsytle=‘k-’,
    linewidth=1, if_radius_normalized = False, ifaxis=False, ifboundary=False,
    ifreverse=False, iflines=False) 绘制特定振型的本征函数， 输入为本征函数文件路径

2024.6.21 新增关于插值与归一化的内容

  • specnm_eig_interpolate(Radius, eigenfunction, order=5, ext=3) 对某一个本征函数序列进行样条插值，使用scipy.interpolate.InterpolatedUnivariateSpline函数。
  • specnm_eig_interpolate_R(splfunctions, split_radius, R) 返回某一半径值位置的本征函数插值得到的值
  • specnm_eig_interpolate_Rs(splfunctions, split_radius, Rs) 使用新的Radius序列计算插值后的本征函数。
  • normalize_toroidal(model_radius, rho, angular_order, eigen_radius, W, mode_freq, W_order=5, W_ext=1, rho_order=1, rho_ext=1) 对环型振型本征函数分量W进行归一化，使用MINEOS的归一化方法
  • normalize_spheroidal(model_radius, rho, angular_order, eigen_radius, U, V, mode_freq, U_order=5, U_ext=1, V_order=5, V_ext=1, rho_order=1, rho_ext=1) 对球型振型本征函数分量W进行归一化，使用MINEOS的归一化方法
  • specnm_eigenfunction_normalize(specnm_output_dir, model_name, model_dir, mode_type, radial_order, angular_order) 通过输入SPECNM本征函数文件进行读取并归一化，给出归一化的本征函数和半径
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
import os
from MINEOS_read import MINEOS_model_read, plot_MINEOS_W_eigenfunction

def get_filename(path):
    if os.path.exists(path):
        # 文件存在，可以进行操作，例如打开读取
        ifexist = True
        return ifexist
    else:
        # 文件不存在，可以选择跳过或者处理错误
        ifexist = False
        return ifexist

def read_SPECNM_mode_frequency_output(mode_file_freq):
    '''
    读取SPECNM振型输出文件；仅输出频率信息
    :param mode_file_freq: string 振型相关文件绝对路径，即该文件所在绝对路径及文件名——完整的路径地址
    :return: list n, l, freq, period
        返回 径向阶数、角度阶数、频率、周期 的列表
    '''
    with open(mode_file_freq, 'r') as file:
        lines = file.readlines()
        n = []
        l = []
        freq = []
        period = []
        for i in range(0, len(lines)):
            if i == 0 or i == 1 or lines[i] == '\n':
                continue
            n.append(float(lines[i][0:6]))
            l.append(float(lines[i][6:12]))
            freq.append(float(lines[i][12:27]))
            period.append(float(lines[i][27:39]))
        file.close()
    return n, l, freq, period

def find_specnm_mode_freq(mode_file, radial_order, angular_order):
    n, l, freq, period = read_SPECNM_mode_frequency_output(mode_file)
    # print(n)
    # print(l)
    freqsss = []
    for i in range(0, len(n)):
        if int(n[i]) == radial_order and int(l[i]) == angular_order:
            freqsss.append(freq[i])
    if freqsss == []:
        print(f'Cannot find this mode n={radial_order}l={angular_order} freq.')
        return 0
    else:
        return freqsss[0]

def SPECNM_mode_name_make(all_dir, model_name, mode_type):
    if mode_type == 'S':
        ppp = 'Rayleigh_output/'
    elif mode_type == 'T':
        ppp = 'Love_output/'
    else:
        ppp = 'Radial_output/'
    mode_dir = all_dir + '/' + model_name + '/' + ppp + model_name + '_' + mode_type + '_output.txt'
    return mode_dir

def read_SPECNM_mode_output(mode_file, modetype):
    '''
    读取SPECNM振型输出文件；相关函数为 Love_output 等SPECNM振型输出函数
    :param mode_file: string 振型相关文件绝对路径，即该文件所在绝对路径及文件名——完整的路径地址
    例如：'/home/shizhiyuanubuntu/PycharmProjects/Seismology/prem_noocean/Rayleigh_output/prem_noocean_S_output.txt'
    :param modetype: string 振型类型，为S、T、R
    :return: list
    若振型类型为S或T，输出列表 n, l, freq, period, ph_v, grp_v, gama, Q, epsilon
    为径向阶数、角度阶数、频率、周期、相速度、群速度、gama、品质因子Q、epsilon
    '''
    if modetype == 'S' or modetype == 'T':
        with open(mode_file, 'r') as file:
            lines = file.readlines()
            n = []
            l = []
            freq = []
            period = []
            ph_v = []
            grp_v = []
            gama = []
            Q = []
            epsilon = []
            '{:>6d}{:>6d}{:>15.5f}{:>12.3f}{:>12.3f}{:>12.3f}{:>12.3f}{:>12.3f}{:>12.3f}'
            for i in range(0, len(lines)):
                if i == 0 or i == 1 or lines[i] == '\n':
                    continue
                n.append(float(lines[i][0:6]))
                l.append(float(lines[i][6:12]))
                freq.append(float(lines[i][12:27]))
                period.append(float(lines[i][27:39]))
                ph_v.append(float(lines[i][39:51]))
                grp_v.append(float(lines[i][51:63]))
                gama.append(float(lines[i][63:75]))
                Q.append(float(lines[i][75:87]))
                epsilon.append(float(lines[i][87:99]))
            file.close()
        return n, l, freq, period, ph_v, grp_v, gama, Q, epsilon
    elif modetype=='R':
        with open(mode_file, 'r') as file:
            lines = file.readlines()
            n = []
            l = []
            freq = []
            period = []
            gama = []
            Q = []
            epsilon = []
            for i in range(0, len(lines)):
                if i == 0 or i == 1 or lines[i] == '\n':
                    continue
                n.append(float(lines[i][0:6]))
                l.append(float(lines[i][6:12]))
                freq.append(float(lines[i][12:27]))
                period.append(float(lines[i][27:39]))
                gama.append(float(lines[i][39:51]))
                Q.append(float(lines[i][51:63]))
                epsilon.append(float(lines[i][63:75]))
            file.close()
        return n, l, freq, period,[], [], gama, Q, epsilon

def find_freq(mode_file, modetype, angular_order, radial_order):
    '''
    寻找SPECNM特定一系列振型的频率，已知振型参数n、l。
    :param mode_file: string 振型相关文件绝对路径，即该文件所在绝对路径及文件名——完整的路径地址
    :param modetype: 振型类型，为S、T、R
    :param angular_order:  1D array, 角度阶数列表，
    :param radial_order:  1D array, 径向阶数列表
    :return: list 频率列表
    '''
    n, l, freq, period = read_SPECNM_mode_frequency_output(mode_file_freq=mode_file)
    out_freq = []
    for i in range(0, len(angular_order)):
        for j in range(0, len(n)):
            if radial_order[i] == n[j] and angular_order[i] == l[j]:
                out_freq.append(freq[j])
    return out_freq

from specnm import rayleigh_fg
from matplotlib import pyplot as plt
from specnm import seismogram
import numpy as np

def open_file_and_plot_eigenfunctions_signal(radial_order, angular_order, data_dir, eigenfunction_type='U', linestyle='k-',ifdr = False, ifnormalized=0):
    '''
    输出单个本征函数文件，绘制图像。输出为plot的一条曲线，需要自行设置figure以及plot.show()
    :param radial_order: 径向阶数n
    :param angular_order: 角度阶数l
    :param data_dir: 本征函数文件所在绝对路径，例如：
    该本征函数文件为radord9angord1.txt，其路径位于：
    '/home/shizhiyuanubuntu/PycharmProjects/Seismology/prem_noocean/Rayleigh_output/eigenfunctions/'
    :param eigenfunction_type: 本征函数的类型，有U、V、P、W四种类型，分别对应S振型的三种分量和T振型一种分量。
    :param linestyle: 颜色和线型，例如'k-'
    :param ifdr: 该文件是否是本征函数的导数（关于radius的导数），默认为False
    :param ifnormalized: 本征函数是否进行归一化 0:不归一化 1：以地表位移为参考 2：以最大值为参考
    :return: 无，仅绘制本征函数曲线，为一个plot函数。需要在函数外部进行figure创建和plt.show
    '''
    if ifdr == True: # 如果绘制的是eigenfunction_dr.
        eigenfunction_file = 'dr_' + eigenfunction_type + '_' + f'radord{radial_order}_angord{angular_order}.txt'
    else:
        eigenfunction_file = eigenfunction_type + '_' + f'radord{radial_order}_angord{angular_order}.txt'
    file = data_dir + eigenfunction_file
    radius_file = data_dir + 'Radius.txt'
    with open(file, 'r') as file:
        lines = file.readlines()
        eigenfunction = []
        eigenfunction_dr = []
        for line in lines:
            eigenfunction.append(float(line))
        file.close()

    with open(radius_file, 'r') as file:
        lines = file.readlines()
        radius = []
        for line in lines:
            radius.append(float(line))
        file.close()

    # print(eigenfunction)
    eigenfunction = np.array(eigenfunction)
    radius = np.array(radius) / 1000
    if ifnormalized == 1:
        plt.plot(radius, eigenfunction / eigenfunction[len(eigenfunction) - 1], linestyle)  # 归一化后
    elif ifnormalized == 2:
        plt.plot(radius, eigenfunction / abs(max(eigenfunction)), linestyle)  # 按照最大值归一化
    else:
        plt.plot(radius, eigenfunction, linestyle)  # 非归一化

import specnm
from matplotlib import pyplot as plt
import numpy as np
def specnm_file_read(wdir, mode_type, radial_order, angular_order):
    '''
    读取SPECNM本征函数文件
    :param wdir:
    :param mode_type:
    :param radial_order:
    :param angular_order:
    :return:
    '''
    weigenfunction_file = mode_type + '_' + f'radord{radial_order}_angord{angular_order}.txt'
    radial_file = 'Radius.txt'
    weigenfunction_file = wdir + weigenfunction_file
    wradial_file = wdir + radial_file
    if get_filename(weigenfunction_file) == False:
        print(f'NO such file of {weigenfunction_file}')
        return [], []
    with open(weigenfunction_file, 'r') as file:
        lines = file.readlines()
        weigenfunction = []
        for line in lines:
            weigenfunction.append(float(line))
        file.close()
    with open(wradial_file, 'r') as file:
        lines = file.readlines()
        wradius = []
        for line in lines:
            wradius.append(float(line))
        file.close()
    return wradius, weigenfunction

def plot_specnm_eigenfunction_with_eig(wradius, weigenfunction , ifnormalized=0,
                         ifnegative=False, mode_type='W', linecolorsytle='k-',
                         linewidth=1, if_radius_normalized = False, ifaxis=False, ifboundary=False,
                         ifreverse=False, iflines=False):
    '''
        绘制特定振型的本征函数
        :param wradius 本征函数对应半径
        :param weigenfunction 本征函数
        :param radial_order: int 径向阶数
        :param angular_order: int 角度阶数
        :param ifnormalized: 是否进行归一化。
            0：不进行归一化；1：使用地表位移作为基准进行归一化； 2: 使用绝对值最大值作为基准进行归一化
            默认为0
        :param ifnegative: 是否对本征函数值进行反转，即加一个负号；
        :param mode_type: string 本征函数振型类型，尽管在路径中已经指出不同的振型类型存放于不同的文件夹下，但是依然需要指定
            振型类型以方便创建需要读取的文件的名称。分为U、V、W三种类型；Radial_output输出仍为U、V分量。
        :param linecolorsytle: string 绘制的本征函数的曲线颜色和线型，默认'k-'
        :param linewidth: float 绘制的本征函数的曲线宽度，默认为1
        :param if_radius_normalized: 输出的图像坐标轴中半径是否归一化；归一化后半径范围为相对值0-1。
        :param ifboundary: 是否在地表处显示一条横线表示在边界面
        :param ifaxis: 是否显示地球660间断面和CMB。
        :param ifreverse: False: 半径为y轴；True：半径为x轴。
        :return: 只进行plot，需要在函数外围设置figure和plt.show()
        '''
    if weigenfunction == []: # 确保不会出现错误
        return

    if ifnegative == True:
        weigenfunction = -np.array(weigenfunction)
    boundary670 = 670
    CMB = 3480
    SURFACE = 6371
    if if_radius_normalized == True:
        wradius = np.array(wradius) / max(np.array(wradius))
        boundary670 = 670 * 1000 / max(np.array(wradius))
        CMB = 3480 * 1000 / max(np.array(wradius))
        SURFACE = 6371 * 1000 / max(np.array(wradius))
    else:
        wradius = np.array(wradius) / 1000

    if ifnormalized == 1:  # 以地表位移作为条件归一化
        if ifreverse == False:
            plt.plot(np.array(weigenfunction) / weigenfunction[len(weigenfunction) - 1], wradius,
                     linecolorsytle, linewidth=linewidth)
        else:
            plt.plot(wradius, np.array(weigenfunction) / weigenfunction[len(weigenfunction) - 1],
                     linecolorsytle, linewidth=linewidth)

        if ifboundary == True:
            if weigenfunction[0] < 0:
                plt.axhline(y=wradius[0], xmin=weigenfunction[0] / weigenfunction[len(weigenfunction) - 1], xmax=0,
                            color='k',
                            linestyle='-', linewidth=1, alpha=1)
            else:
                plt.axhline(y=wradius[0], xmin=0, xmax=weigenfunction[0] / weigenfunction[len(weigenfunction) - 1],
                            color='k',
                            linestyle='-', linewidth=1, alpha=1)
            if weigenfunction[len(weigenfunction) - 1] < 0:
                plt.axhline(y=wradius[len(wradius) - 1],
                            xmin=weigenfunction[len(weigenfunction) - 1] / max(weigenfunction),
                            xmax=0,
                            color='k', linestyle='-', linewidth=1, alpha=1)
            else:
                plt.axhline(y=wradius[len(wradius) - 1], xmin=0,
                            xmax=weigenfunction[len(weigenfunction) - 1] / weigenfunction[len(weigenfunction) - 1],
                            color='k', linestyle='-', linewidth=1, alpha=1)
        if ifaxis == True:
            plt.text(-1, CMB, 'CMB')
            plt.text(-1, max(wradius)-boundary670, '670')
            plt.plot([-1, -0.7], [max(wradius) - boundary670, max(wradius) - boundary670], 'k-')
            plt.text(-1, SURFACE, '0')

    elif ifnormalized == 2:  # 当使用绝对值最大值作为归一化条件时
        if ifreverse == False:
            plt.plot(np.array(weigenfunction) / max(abs(np.array(weigenfunction))), wradius, linecolorsytle,
                     linewidth=linewidth)
        else:
            plt.plot(wradius, np.array(weigenfunction) / max(abs(np.array(weigenfunction))),
                     linecolorsytle, linewidth=linewidth)

        if ifboundary == True:
            if weigenfunction[0] < 0:
                # plt.axhline(y=wradius[0]/1000, xmin=weigenfunction[0]/max(abs(np.array(weigenfunction))), xmax=0, color='k',
                #            linestyle='-', linewidth=1, alpha=1)
                plt.plot([weigenfunction[0] / max(abs(np.array(weigenfunction))), 0],
                         [wradius[0], wradius[0]], 'k-')
            else:
                # plt.axhline(y=wradius[0]/1000, xmin=0, xmax=weigenfunction[0] / max(abs(np.array(weigenfunction))), color='k',
                #            linestyle='-', linewidth=1, alpha=1)
                plt.plot([0, weigenfunction[0] / max(abs(np.array(weigenfunction)))],
                         [wradius[0], wradius[0]], 'k-')

            if weigenfunction[len(weigenfunction) - 1] < 0:
                # plt.axhline(y=wradius[len(wradius)-1]/1000, xmin=weigenfunction[len(weigenfunction) - 1] / max(abs(np.array(weigenfunction))),
                #             xmax=0, color='k', linestyle='-', linewidth=1, alpha=1)
                plt.plot([weigenfunction[len(weigenfunction) - 1] / max(abs(np.array(weigenfunction))), 0],
                         [wradius[len(wradius) - 1], wradius[len(wradius) - 1]], 'k-')
            else:
                # plt.axhline(y=wradius[len(wradius)-1]/1000, xmin=0,
                #            xmax=weigenfunction[len(weigenfunction) - 1] / max(abs(np.array(weigenfunction))),
                #            color='k', linestyle='-', linewidth=1, alpha=1)\
                plt.plot([0, weigenfunction[len(weigenfunction) - 1]] / max(abs(np.array(weigenfunction))),
                         [wradius[len(wradius) - 1], wradius[len(wradius) - 1]], 'k-')
        if ifaxis == True:
            plt.text(-1, CMB, 'CMB')
            plt.text(-1, max(wradius) - boundary670, '670')
            plt.plot([-1,-0.7], [max(wradius) - boundary670, max(wradius) - boundary670], 'k-')
            plt.text(-1, SURFACE, '0')


    else:  # 不进行归一化
        if ifreverse == False:
            plt.plot(np.array(weigenfunction), wradius, linecolorsytle, linewidth=linewidth)
        else:
            plt.plot(wradius, np.array(weigenfunction),
                     linecolorsytle, linewidth=linewidth)

        if ifboundary == True:

            if weigenfunction[0] < 0:
                plt.axhline(y=wradius[0], xmin=weigenfunction[0],
                            xmax=0, color='k',
                            linestyle='-', linewidth=1, alpha=1)
            else:
                plt.axhline(y=wradius[0], xmin=0,
                            xmax=weigenfunction[0], color='k',
                            linestyle='-', linewidth=1, alpha=1)
            if weigenfunction[len(weigenfunction) - 1] < 0:
                plt.axhline(y=wradius[len(wradius) - 1],
                            xmin=weigenfunction[len(weigenfunction) - 1],
                            xmax=0,
                            color='k', linestyle='-', linewidth=1, alpha=1)
            else:
                plt.axhline(y=wradius[len(wradius) - 1], xmin=0,
                            xmax=weigenfunction[len(weigenfunction) - 1],
                            color='k', linestyle='-', linewidth=1, alpha=1)
        if ifaxis == True:
            plt.text(-1.2 * abs(min(weigenfunction)), CMB, 'CMB')
            plt.text(-1.2 * abs(min(weigenfunction)), max(wradius) - boundary670, '670')
            plt.plot([-1.2 * abs(min(weigenfunction)), -1 * abs(min(weigenfunction))], [max(wradius) - boundary670, max(wradius) - boundary670], 'k-')
            plt.text(-1.2 * abs(min(weigenfunction)), SURFACE, '0')

def plot_W_eigenfunction(wdir, radial_order, angular_order, ifnormalized=0,
                         ifnegative=False, mode_type='W', linecolorsytle='k-',
                         linewidth=1, if_radius_normalized = False, ifaxis=False, ifboundary=False,
                         ifreverse=False, iflines=False):
    '''
    绘制特定振型的本征函数
    :param wdir:  本征函数文件所在路径，例如W_radord1angord2.txt位于该路径下：
    '/home/shizhiyuanubuntu/PycharmProjects/Seismology/prem_noocean/Love_output/eigenfunctions/'
    :param radial_order: int 径向阶数
    :param angular_order: int 角度阶数
    :param ifnormalized: 是否进行归一化。
        0：不进行归一化；1：使用地表位移作为基准进行归一化； 2: 使用绝对值最大值作为基准进行归一化
        默认为0
    :param ifnegative: 是否对本征函数值进行反转，即加一个负号；
    :param mode_type: string 本征函数振型类型，尽管在路径中已经指出不同的振型类型存放于不同的文件夹下，但是依然需要指定
        振型类型以方便创建需要读取的文件的名称。分为U、V、W三种类型；Radial_output输出仍为U、V分量。
    :param linecolorsytle: string 绘制的本征函数的曲线颜色和线型，默认'k-'
    :param linewidth: float 绘制的本征函数的曲线宽度，默认为1
    :param if_radius_normalized: 输出的图像坐标轴中半径是否归一化；归一化后半径范围为相对值0-1。
    :param ifboundary: 是否在地表处显示一条横线表示在边界面
    :param ifaxis: 是否显示地球660间断面和CMB。
    :param ifreverse: False: 半径为y轴；True：半径为x轴。
    :return: 只进行plot，需要在函数外围设置figure和plt.show()
    '''
    wradius, weigenfunction = specnm_file_read(wdir, mode_type, radial_order, angular_order)
    if weigenfunction == []: # 确保不会出现错误
        print(f'NO such file of {radial_order}{mode_type}{angular_order}')
        return
    plot_specnm_eigenfunction_with_eig(wradius=wradius, weigenfunction=weigenfunction, ifnormalized=ifnormalized,
                                       ifnegative=ifnegative, mode_type=mode_type, linecolorsytle=linecolorsytle,
                                       linewidth=linewidth, if_radius_normalized=if_radius_normalized, ifaxis=ifaxis, ifboundary=ifboundary,
                                       ifreverse=ifreverse, iflines=iflines)


from scipy.interpolate import InterpolatedUnivariateSpline
from scipy.interpolate import splprep, splev
def specnm_eig_interpolate(Radius, eigenfunction, order=5, ext=3):
    '''
    对某一个本征函数序列进行样条插值，使用scipy.interpolate.InterpolatedUnivariateSpline函数。

    :param Radius: numpy.array 输入的SPECNM本征函数半径序列
    :param eigenfunction: numpy.array 输入的SPECNM本征函数序列
    :param order: 样条插值阶数，为数字1-5，默认5
    :param ext: ext=3为插值函数在超出输入范围时进行外插值
    :return:
        spls 列表，为每一个split_radius切片对应的插值函数，spls[index](Radius_new)即可调用该函数，
        注意：index为相应Radius_new对应的split_radius切片的索引，因此Radius_new也需要切片
        split_radius, 2D list， 每一个包含根据Raius序列中，根据相同半径切片，将其分为n个间断面。
        split_eigenfunction， 2D list， 每一个包含根据Raius序列中，根据相同半径切片，将其分为n个间断面，这一个间断面对应的原始eigenfunction
    '''


    if eigenfunction == [] or Radius == []:
        print('There is no file.')
        spls=[]
        split_radius=[]
        split_eigenfunction=[]
        return spls, split_radius, split_eigenfunction

    Radius = np.array(Radius)
    eigenfunction = np.array(eigenfunction)

    # 记录Radius相同位置的索引
    radius_index = [0]
    if Radius[2] > Radius[3]:
        print('Radius is not increasing.')
        Radius = np.flipud(Radius)
        eigenfunction = np.flipud(eigenfunction)

    for i in range(0, len(Radius)-1):
        if Radius[i] == Radius[i+1]:
            radius_index.append(i)

    radius_index.append(len(Radius)-1)
    # 切片
    split_radius = []
    split_eigenfunction = []
    for i in range(0, len(radius_index)-1):
        if i == 0:
            split_radius.append(Radius[ radius_index[i] : radius_index[i+1]+1 ])
            split_eigenfunction.append(eigenfunction[ radius_index[i] : radius_index[i+1]+1 ])
        else:
            split_radius.append(Radius[ radius_index[i] + 1 : radius_index[i+1] + 1 ])
            split_eigenfunction.append(eigenfunction[radius_index[i] + 1 : radius_index[i + 1] + 1])
    spls = []
    for i in range(0, len(split_radius)):
        if len(split_radius[i]) != len(split_eigenfunction[i]):
            print(len(split_radius[i]))
            print(len(split_eigenfunction[i]))
            print('danger! split_radius and split_eigenfunction has no same length')
        try:
            spls.append(InterpolatedUnivariateSpline(split_radius[i], split_eigenfunction[i], k=order, ext=ext))
        except ValueError:
            print('In SPECNM_read specnm_eig_interpolate: spls.append(InterpolatedUnivariateSpline(split_radius[i], split_eigenfunction[i], k=order, ext=ext))')
            print('Your spline order: k may be too high, try lower. The range is 1-5.')

    # spl = InterpolatedUnivariateSpline(Radius, eigenfunction, k=5, ext=3)
    # tck, u = splprep([Radius, eigenfunction], s=0)
    # new_points = splev(u, tck)

    return spls, split_radius, split_eigenfunction

def specnm_eig_interpolate_R(splfunctions, split_radius, R):
    '''

    :param splfunctions:
    :param split_radius: 切片半径
    :param R: 半径 单位m
    :return: 返回某一个R位置的本征函数，如果R不在原来的范围内，值为0.
    '''
    if splfunctions == []:
        eigenfunction_in_R = np.nan
        return eigenfunction_in_R
    if R <= split_radius[0][0] or R >= split_radius[len(split_radius)-1][len(split_radius[len(split_radius)-1])-1] :
        eigenfunction_in_R = 0
        return eigenfunction_in_R
    for j in range(0, len(split_radius)):
        if R >= split_radius[j][0] and R <= split_radius[j][len(split_radius[j])-1]:
            eigenfunction_in_R = splfunctions[j](R)
            return eigenfunction_in_R

def specnm_eig_interpolate_Rs(splfunctions, split_radius, Rs):
    '''
    使用新的Radius序列计算插值后的本征函数。
    :param splfunctions: list [spl1, spl2, ...] specnm_eig_interpolate返回的插值函数列表
    :param split_radius: list [radius range1, radius range2, ...] specnm_eig_interpolate返回的插值半径根据不连续增长的半径切片的列表
    （由于有两个相同的半径值，因此无法进行插值，需要进行切片之后才可插值）
    :param Rs: numpy.array 插值所需要的半径列表
    :return: 插值后的序列
    '''
    eigenfunction_spline_inter = []
    for R in Rs:
        if splfunctions == []:
            eigenfunction_spline_inter.append(np.nan)
        else:
            eigenfunction_spline_inter.append(specnm_eig_interpolate_R(splfunctions, split_radius, R))
    return eigenfunction_spline_inter


import numpy as np
from scipy.integrate import quad

def normalize_toroidal(model_radius, rho, angular_order, eigen_radius, W, mode_freq, W_order=5, W_ext=1, rho_order=1, rho_ext=1):
    '''
    对环型振型本征函数分量W进行归一化，使用MINEOS的归一化方法
    :param model_radius: 使用的模型的半径序列，例如MINEOS模型：prem_noocean
    :param rho: 使用的模型的密度序列，例如MINEOS模型：prem_noocean
    :param eigen_radius: 使用的SPECNM计算得到的本征函数对应的半径序列
    :param W: SPECNM计算得到的环型振型本征函数分量W
    :param mode_freq: 这个振型的本征频率
    :param W_order: 对本征函数W分量进行样条插值的平滑度，取值为1-5
    :param W_ext: 对本征函数W分量进行样条插值的外插值方法
        if ext=0 or ‘extrapolate’, return the extrapolated value.
        if ext=1 or ‘zeros’, return 0
        if ext=2 or ‘raise’, raise a ValueError
        if ext=3 of ‘const’, return the boundary value.
    :param rho_order: 对模型的密度进行样条插值的平滑度
    :param rho_ext: 对模型的密度进行样条插值的外插值方法
    :return:
    归一化之后的W分量和对应半径(m)
    '''
    rn = max(model_radius)  # rn为地表半径（地球半径）

    # 最大重力加速度
    G = 6.6723E-11
    tau = 1000.0
    rhobar = 5515.0  # 地球参数
    con = np.pi * G
    gn = con * rhobar * rn
    # bigg, tau, rhobar / 6.6723e-11, 1000.0, 5515.0
    # / pi = 4.0　d0 * datan(1.d0)   con = pi * bigg
    vn = np.sqrt(gn * rn) #　最大速度
    wn = vn / rn  # 最大频率
    accn = 1.e+20 / (rhobar * rn ** 4)  # 最大加速度

    Wspls, Wsplit_radius, Wsplit_eigenfunction = (
        specnm_eig_interpolate(Radius=eigen_radius, eigenfunction=W,
                               order=W_order, ext=W_ext))
    rho_spls, rho_split_radius, rho_split_eigenfunction = (
        specnm_eig_interpolate(Radius=model_radius, eigenfunction=rho,
                               order=rho_order, ext=rho_ext))

    Wintegrand = lambda r: (specnm_eig_interpolate_R(splfunctions=rho_spls,
                    split_radius=rho_split_radius, R=r) / 1000 * angular_order * (angular_order + 1) *
        specnm_eig_interpolate_R(splfunctions=Wspls, split_radius=Wsplit_radius,
                                 R=r)  ** 2 * pow(r/rn, 2))
    Wintegral, _ = quad(Wintegrand, 0, rn)  # 从0到rn进行积分
    Wnormalization_factor = np.sqrt(1 / ( pow(mode_freq/1000, 2) * Wintegral))
    print(Wnormalization_factor)
    W = np.array(W) * Wnormalization_factor
    return W, np.array(eigen_radius)

def normalize_spheroidal(model_radius, rho, angular_order, eigen_radius, U, V, mode_freq, U_order=5, U_ext=1,  V_order=5, V_ext=1, rho_order=1, rho_ext=1):
    '''
    与normalize_toroidal相同
    :param model_radius:
    :param rho:
    :param angular_order:
    :param eigen_radius:
    :param U:
    :param V:
    :param mode_freq:
    :param U_order:
    :param U_ext:
    :param V_order:
    :param V_ext:
    :param rho_order:
    :param rho_ext:
    :return:
    '''
    rn = max(model_radius)  # rn为地表半径（地球半径）
    # 最大重力加速度
    G = 6.6723E-11
    tau = 1000.0
    rhobar = 5515.0  # 地球参数
    con = np.pi * G
    gn = con * rhobar * rn
    # bigg, tau, rhobar / 6.6723e-11, 1000.0, 5515.0
    # / pi = 4.0　d0 * datan(1.d0)   con = pi * bigg
    vn = np.sqrt(gn * rn) #　最大速度
    wn = vn / rn  # 最大频率
    accn = 1.e+20 / (rhobar * rn ** 4)  # 最大加速度

    Uspls, Usplit_radius, Usplit_eigenfunction = (
        specnm_eig_interpolate(Radius=eigen_radius, eigenfunction=U,
                               order=U_order, ext=U_ext))

    V = np.array(V) / np.sqrt(angular_order * (angular_order + 1))

    Vspls, Vsplit_radius, Vsplit_eigenfunction = (
        specnm_eig_interpolate(Radius=eigen_radius, eigenfunction=np.array(V),
                               order=V_order, ext=V_ext))
    rho_spls, rho_split_radius, rho_split_eigenfunction = (
        specnm_eig_interpolate(Radius=model_radius, eigenfunction=rho,
                               order=rho_order, ext=rho_ext))

    UVintegrand = lambda r: (specnm_eig_interpolate_R(splfunctions=rho_spls,
                    split_radius=rho_split_radius, R=r) / 1000 *
                             (pow(specnm_eig_interpolate_R(splfunctions=Uspls, split_radius=Usplit_radius, R=r), 2) +
                               angular_order * (angular_order + 1)* pow(specnm_eig_interpolate_R(splfunctions=Vspls, split_radius=Vsplit_radius, R=r), 2) ) * pow(r/rn, 2))
    UVintegral, _ = quad(UVintegrand, 0, rn)  # 从0到rn进行积分
    print(mode_freq)
    UVnormalization_factor = np.sqrt(1 / ( pow(mode_freq/1000, 2) * UVintegral))
    print(UVnormalization_factor)
    U = np.array(U) * UVnormalization_factor
    V = np.array(V) * UVnormalization_factor
    return U, V, np.array(eigen_radius)

def specnm_eigenfunction_normalize(specnm_output_dir, model_name, model_dir, mode_type, radial_order, angular_order):
    '''
    通过输入SPECNM本征函数文件进行读取并归一化，给出归一化的本征函数和半径
    :param specnm_output_dir:  SPECNM输出的路径，这个是SPECNM输出总路径，如：
        '/home/shizhiyuanubuntu/PycharmProjects/Seismology/Suart_Russell/'
        举例旗下包含目录： 注意！！ 这个目录中的文件夹需要提前自行创建，并且严格按照下列名称命名。 /名称  表示文件夹，/名称.txt 表示文件
            -/prem_noocean  由模型名称命名的文件夹
                -/Rayleigh_output   球型振型输出
                    -/eigenfunctions 存储本征函数的文件夹，U，V分量位于这一个文件夹里
                    -/eigenfunctions_dr 存储本征函数的导数的文件夹
                    -/prem_noocean_S_output.txt 存储振型的相关信息的文件
                    -/prem_noocean_S_output_onlyfreq.txt  只存储振型的相关本征频率信息的文件
                -/Love_output  同上，为环型振型输出
                    -/eigenfunctions
                    -/eigenfunctions_dr
                    -/prem_noocean_T_output.txt
                    -/prem_noocean_T_output_onlyfreq.txt
                -/Radial_output 同上，为l=0时的输出
                    -/eigenfunctions
                    -/eigenfunctions_dr
                    -/prem_noocean_R_output.txt
                    -/prem_noocean_R_output_onlyfreq.txt
    :param model_name: 使用的模型名称，不要添加后缀。例如‘prem_noocean’，且模型必须为.txt文件
    :param model_dir: 模型所位于的路径，例如'/home/shizhiyuanubuntu/PycharmProjects/Seismology/models/old_models/'
    :param mode_type: 振型类型，有 U，V，W三种
    :param radial_order: 径向阶数
    :param angular_order: 角度阶数
    :return:
        mode_type='U'，输出 本征函数对应半径radius, 归一化之后的U, 归一化之后的V
        mode_type='W'，输出 本征函数对应半径radius, 归一化之后的W
    '''
    # 使用的模型路径
    model_dir = model_dir + model_name + '.txt'
    model_pro = MINEOS_model_read(model_dir)
    # 使用的模型的半径序列和密度序列
    model_radius = model_pro[0]
    rho = model_pro[1]
    if mode_type == 'U' or mode_type == 'V':
        wdir_UV = specnm_output_dir + model_name + '/Rayleigh_output/eigenfunctions/'
        mode_file = specnm_output_dir + model_name + '/Rayleigh_output/' + model_name + '_S_output.txt'
        # 读取振型本征频率信息
        mode_freq = find_specnm_mode_freq(mode_file, radial_order=radial_order, angular_order=angular_order)
        # 读取U振型本征函数信息
        Uradius, Ueigenfunction = specnm_file_read(wdir=wdir_UV, mode_type='U', radial_order=radial_order,
                                                angular_order=angular_order)
        # 读取V振型本征函数信息
        Vradius, Veigenfunction = specnm_file_read(wdir=wdir_UV, mode_type='V', radial_order=radial_order,
                                                   angular_order=angular_order)
        # 进行归一化，返回的是归一化之后的U、V、radius分量
        U, V, radius = normalize_spheroidal(model_radius=model_radius, rho=rho, eigen_radius=Uradius, U=Ueigenfunction,
                                            V=Veigenfunction, mode_freq=mode_freq, angular_order=angular_order)
        return radius, U, V
    elif mode_type == 'W':
        wdir_W = specnm_output_dir + model_name + '/Love_output/eigenfunctions/'
        mode_file = specnm_output_dir + model_name + '/Love_output/'+ model_name + '_T_output.txt'
        # 读取振型本征频率信息
        mode_freq = find_specnm_mode_freq(mode_file, radial_order=radial_order, angular_order=angular_order)
        # 读取W振型本征函数信息
        Wradius, Weigenfunction = specnm_file_read(wdir=wdir_W, mode_type='W', radial_order=radial_order,
                                                   angular_order=angular_order)
        W, radius = normalize_toroidal(model_radius=model_radius, rho=rho, eigen_radius=Wradius, W=Weigenfunction,
                                       mode_freq=mode_freq, angular_order=angular_order)
        return radius, W
    else:
        print('Wrong mode type input.')
        return [] , []

MINEOS_read.py

关于MINEOS软件输入与输出、模型相关的函数：

函数列表：

  • get_filename(path)
  • MINEOS_model_read(modelfile) 读取MINEOS模型，即MINEOS DEMO中给出的标准模型
  • put_freq_in_same(n, l, freq, nmax=10000, lmax=101, ifnan=True) 将MINEOS输出的振型文件的一个参数整合为一个列表，按照n、l排列
  • read_mineos_modefile(modefile_dir, model_name, mode_type) 读取MINEOS输出的振型文件中的相关信息并输出
  • put_0_velocity(n,frequency, ph_vel, grop_vel, ni = 0) 取同一个径向阶数值的相关参数。
  • mineos_file_read(MINEOS_eigenfunction_dir, Modetype, radial_order, angular_order, Mineos_file=[]) 读取MINEOS单个本征函数文件
  • numofzero(numb) 添加0,因为MINEOS输出的本征函数文件命名规则需要添加很多0,因此使用该函数创建本征函数名以供读取使用。
  • MINEOS_eigenfunction_file_name(radial_order, angular_order) 创建MINEOS本征函数文件名
  • plot_MINEOS_W_eigenfunction(wdir, radial_order, angular_order, ifnormalized=0,
    ifnegative=False, mode_type=‘W’, linecolorsytle=‘k-’,
    linewidth=1, if_radius_normalized = False, ifaxis=False, ifboundary=False,
    ifreverse=False, iflines=False) 绘制特定振型的本征函数 2024.6.21 已解决ifnormalized = 0时存在BUG
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
import os

def get_filename(path):
    if os.path.exists(path):
        # 文件存在，可以进行操作，例如打开读取
        ifexist = True
        return ifexist
    else:
        # 文件不存在，可以选择跳过或者处理错误
        ifexist = False
        return ifexist

def MINEOS_model_read(modelfile):
    '''
    :param modelfile: Absolute file path of model
    :return: A 2D list have:
    [Radius, rho, vp, vs, qp, qs, eta]
    '''
    if get_filename(modelfile) == False:
        print(f'NO such file of {modelfile}')
        return []
    with open(modelfile, 'r') as file:
        lines = file.readlines()
        Radius = []
        rho = []
        vp = []
        vs = []
        qp = []
        qs = []
        eta = []

        for i in range(0, len(lines)):
            if i == 0 or i==1 or i==2 or lines[i]=='\n':
                continue
            Radius.append(float(lines[i][0:8]))
            rho.append(float(lines[i][8:17]))
            vp.append(float(lines[i][17:26]))
            vs.append(float(lines[i][26:35]))
            qp.append(float(lines[i][35:44]))
            qs.append(float(lines[i][44:53]))
            eta.append(float(lines[i][73:80]))
        file.close()
    out = []
    out.append(Radius)
    out.append(rho)
    out.append(vp)
    out.append(vs)
    out.append(qp)
    out.append(qs)
    out.append(eta)
    return out

import matplotlib.pyplot as plt
import obspy
import numpy as np
from obspy import taup

def put_freq_in_same(n, l, freq, nmax=10000, lmax=101, ifnan=True):
    '''
    将MINEOS输出的振型文件的一个参数整合为一个列表，按照n、l排列
    :param n: list 列表 读取输出的径向阶数列表n
    :param l: list 列表  读取输出的角度阶数列表l
    :param freq: list 列表  读取输出的频率列表，输入的freq列表为其他项也可以，例如群速度或者相速度
    :param nmax: int 整数  读取的输出径向阶数列表的最大长度，可以尽量设置大一些以保证所有内容存在
    :param lmax: int 整数  读取的输出角度阶数列表的最大长度，可以尽量设置大一些以保证所有内容存在
    :return: list 列表  freqs_all 索引值为[n, l]的频率列表，如freqs_all[1,2]为n=1,l=2振型的本征频率
    '''
    freqs_all = np.zeros([nmax, lmax])
    for i in range(0, len(n)):
        freqs_all[int(n[i]), int(l[i])] = freq[i]
    if ifnan==True:
        freqs_all[freqs_all == 0] = np.nan
    return freqs_all

def read_mineos_modefile(modefile_dir, model_name, mode_type):
    '''
    读取MINEOS输出的振型文件中的相关信息并输出。
    :param modefile_dir: string 字符串 读取的振型文件所在文件夹绝对路径，
        例如：'/home/shizhiyuanubuntu/Documents/mineos-1.0.2/DEMO/prem_sim_nolvz',
        此路径末尾不包含'/'，振型文件就在该文件夹下
    :param model_name: string 字符串 MINEOS计算使用的模型名称，例如：prem_mean_noLVZ
    :param mode_type: string 字符串 振型类型，如S、T、R。表征球型、环型、径向振型的振型输出文件
    :return: n, l, frequency, phase_velocity, group_velocity
            输出为五个列表(list)，分别为 径向阶数、角度阶数、频率、相速度、群速度
    '''
    modefile = model_name + '_' + mode_type
    modefile_alldir = modefile_dir + '/' + modefile
    print(modefile_alldir)
    with open(modefile_alldir, 'r') as file:
        lines = file.readlines()
        n = []
        l = []
        frequency = []
        group_velocity = []
        phase_velocity = []
        for i in range(0, len(lines)):
            if lines[i][6:10] == 'mode':
                index = i+2
        for i in range(index, len(lines)):
            n.append(float(lines[i][0:5]))
            l.append(float(lines[i][8:12]))
            phase_velocity.append(float(lines[i][13:24]))
            frequency.append(float(lines[i][25:40]))
            group_velocity.append(float(lines[i][57:72]))
    return n, l, frequency, phase_velocity, group_velocity

def put_0_velocity(n,frequency, ph_vel, grop_vel, ni = 0):
    '''
    取同一个径向阶数值的相关参数。
    :param n: list 振型文件输出的径向阶数
    :param frequency: list 振型文件输出的频率
    :param ph_vel: list 振型文件输出的相速度
    :param grop_vel: list 振型文件输出的群速度
    :param ni: int 输出的径向阶数的值。
    :return: freq, ph_vel0, grop_vel0
        三个列表(list)，分别对应ni值对应的本征频率、相速度、群速度
    '''
    ph_vel0 = []
    grop_vel0 = []
    freq = []
    for i in range(0, len(n)):
        if n[i] == ni:
            freq.append(frequency[i])
            ph_vel0.append(ph_vel[i])
            grop_vel0.append(grop_vel[i])
    return freq, ph_vel0, grop_vel0

import numpy as np
import matplotlib.pyplot as plt
import os

def get_filename(path):
    if os.path.exists(path):
        # 文件存在，可以进行操作，例如打开读取
        ifexist = True
        return ifexist
    else:
        # 文件不存在，可以选择跳过或者处理错误
        ifexist = False
        return ifexist

def mineos_file_read(MINEOS_eigenfunction_dir, Modetype, radial_order, angular_order, Mineos_file=[]):
    '''
     读取MINEOS单个本征函数文件；
    :param MINEOS_eigenfunction_dir: string MINEOS 文件绝对路径，
        例如：'/home/shizhiyuanubuntu/PycharmProjects/Seismology/MINEOS_eigenfunctions/'
        程序内部根据Modetype自动添加'test_T_ASC/'或'test_S_ASC/'等文件夹名。
    :param Modetype: string 振型类型，S、T、R
    :param radial_order: int 径向阶数
    :param angular_order: int 角度阶数
    :param Mineos_file: Mineos本征函数文件的名称.ASC，可以不指定该名称。
    如果指定则按照改名称选择，忽略radial_oder和angular_order参数；
    如果不指定则需要指定radial_order、angular_order
    :return: 返回值均为list
    若振型类型为S，返回：Radius, U_eigenfunction, V_eigenfunction, P_eigenfunction
    分别表示半径列表、U分量、V分量、P分量（考虑重力相关参数）
    若振型类型为T，返回：Radius, W_eigenfunction
    若振型类型为R，返回：Radius, R_eigenfunction
    '''
    radial_orderstr, angular_orderstr = MINEOS_eigenfunction_file_name(radial_order, angular_order)
    if Mineos_file == []:
        if Modetype == 'T' or Modetype == 't' or Modetype == 'W':
            modename = 'T.' + radial_orderstr + '.' + angular_orderstr + '.ASC'
            print(modename)
            T_MINEOS_eigenfunction_file = MINEOS_eigenfunction_dir + 'test_T_ASC/' + modename
            if get_filename(T_MINEOS_eigenfunction_file) == False:
                print(f'NO such file of {T_MINEOS_eigenfunction_file}')
                return [], []
            with open(T_MINEOS_eigenfunction_file, 'r') as file:
                lines = file.readlines()
                W_eigenfunction = []
                Radius = []
                for i in range(0, len(lines)):
                    if i == 0:
                        continue
                    Radius.append(float(lines[i][0:8]))
                    W_eigenfunction.append(float(lines[i][8:23]))
                file.close()
            return Radius, W_eigenfunction

        elif Modetype == 'S' or Modetype == 's' or Modetype == 'U'  or Modetype == 'V':
            modename = 'S.' + radial_orderstr + '.' + angular_orderstr + '.ASC'
            print(modename)
            S_MINEOS_eigenfunction_file = MINEOS_eigenfunction_dir + 'test_S_ASC/' + modename
            if get_filename(S_MINEOS_eigenfunction_file) == False:
                print(f'NO such file of {S_MINEOS_eigenfunction_file}')
                return [], [], [], []
            with open(S_MINEOS_eigenfunction_file, 'r') as file:
                lines = file.readlines()
                U_eigenfunction = []
                V_eigenfunction = []
                P_eigenfunction = []
                Radius = []
                for i in range(0, len(lines)):
                    if i == 0:
                        continue
                    else:
                        Radius.append(float(lines[i][0:8]))
                        U_eigenfunction.append(float(lines[i][8:23]))
                        V_eigenfunction.append(float(lines[i][38:53]))
                        P_eigenfunction.append(float(lines[i][68:83]))
                file.close()
            return Radius, U_eigenfunction, V_eigenfunction, P_eigenfunction
        elif Modetype == 'R' or Modetype == 'r':
            modename = 'S.' + radial_orderstr + '.' + angular_orderstr + '.ASC'
            print(modename)
            R_MINEOS_eigenfunction_file = MINEOS_eigenfunction_dir + 'test_R_ASC/' + modename
            if get_filename(R_MINEOS_eigenfunction_file) == False:
                print(f'NO such file of {R_MINEOS_eigenfunction_file}')
                return [], []
            with open(R_MINEOS_eigenfunction_file, 'r') as file:
                lines = file.readlines()
                R_eigenfunction = []
                Radius = []
                for i in range(0, len(lines)):
                    if i == 0:
                        continue
                    Radius.append(float(lines[i][0:8]))
                    R_eigenfunction.append(float(lines[i][8:23]))
                file.close()
            return Radius, R_eigenfunction
        else:
            if Modetype == 'T' or Modetype == 't'  or Modetype == 'W':
                T_MINEOS_eigenfunction_file = MINEOS_eigenfunction_dir + 'test_T_ASC/' + Mineos_file
                if get_filename(T_MINEOS_eigenfunction_file) == False:
                    print(f'NO such file of {T_MINEOS_eigenfunction_file}')
                    return [], []
                with open(T_MINEOS_eigenfunction_file, 'r') as file:
                    lines = file.readlines()
                    W_eigenfunction = []
                    Radius = []
                    for i in range(0, len(lines)):
                        if i == 0:
                            continue
                        Radius.append(float(lines[i][0:8]))
                        W_eigenfunction.append(float(lines[i][8:23]))
                    file.close()
                return Radius, W_eigenfunction

            elif Modetype == 'S' or Modetype == 's' or Modetype == 'U' or Modetype == 'V':
                S_MINEOS_eigenfunction_file = MINEOS_eigenfunction_dir + 'test_S_ASC/' + Mineos_file
                if get_filename(S_MINEOS_eigenfunction_file) == False:
                    print(f'NO such file of {S_MINEOS_eigenfunction_file}')
                    return [], [], [], []
                with open(S_MINEOS_eigenfunction_file, 'r') as file:
                    lines = file.readlines()
                    U_eigenfunction = []
                    V_eigenfunction = []
                    P_eigenfunction = []
                    Radius = []
                    for i in range(0, len(lines)):
                        if i == 0:
                            continue
                        else:
                            Radius.append(float(lines[i][0:8]))
                            U_eigenfunction.append(float(lines[i][8:23]))
                            V_eigenfunction.append(float(lines[i][38:53]))
                            P_eigenfunction.append(float(lines[i][68:83]))
                    file.close()
                return Radius, U_eigenfunction, V_eigenfunction, P_eigenfunction
            elif Modetype == 'R' or Modetype == 'r':
                R_MINEOS_eigenfunction_file = MINEOS_eigenfunction_dir + 'test_R_ASC/' + Mineos_file
                if get_filename(R_MINEOS_eigenfunction_file) == False:
                    print(f'NO such file of {R_MINEOS_eigenfunction_file}')
                    return [], []
                with open(R_MINEOS_eigenfunction_file, 'r') as file:
                    lines = file.readlines()
                    R_eigenfunction = []
                    Radius = []
                    for i in range(0, len(lines)):
                        if i == 0:
                            continue
                        Radius.append(float(lines[i][0:8]))
                        R_eigenfunction.append(float(lines[i][8:23]))
                    file.close()
                return Radius, R_eigenfunction

def numofzero(numb):
    '''
    添加0,因为MINEOS输出的本征函数文件命名规则需要添加很多0,因此使用该函数创建本征函数名以供读取使用。
    :param numb: int 输入的数字
    :return: string 增加0之后的数字
    '''
    if numb < 10:
        itsstr = '000000' + str(numb)
    elif numb < 100:
        itsstr = '00000' + str(numb)
    elif numb < 1000:
        itsstr = '0000' + str(numb)
    elif numb < 10000:
        itsstr = '000' + str(numb)
    return itsstr


def MINEOS_eigenfunction_file_name(radial_order, angular_order):
    '''
    创建MINEOS本征函数名
    :param radial_order: int 径向阶数
    :param angular_order: int 角度阶数
    :return: string
        radial_orderstr: 径向阶数补0后的字符串类型
        angular_orderstr: 径向阶数补0后的字符串类型
    '''
    radial_orderstr = numofzero(radial_order)
    angular_orderstr = numofzero(angular_order)
    return radial_orderstr, angular_orderstr

def plot_MINEOS_W_eigenfunction(wdir, radial_order, angular_order, ifnormalized=0,
                         ifnegative=False, mode_type='W', linecolorsytle='k-',
                         linewidth=1, if_radius_normalized = False, ifaxis=False, ifboundary=False,
                         ifreverse=False, iflines=False):
    # 注意 不进行归一化，ifnormalized = 0时存在BUG
    '''
    绘制特定振型的本征函数
    :param wdir:  本征函数文件所在路径，例如test_S_ASC位于该路径下：
    例如：'/home/shizhiyuanubuntu/PycharmProjects/Seismology/MINEOS_eigenfunctions/'
        程序内部根据Modetype自动添加'test_T_ASC/'或'test_S_ASC/'等文件夹名。
    :param radial_order: int 径向阶数
    :param angular_order: int 角度阶数
    :param ifnormalized: 是否进行归一化。
        0：不进行归一化；1：使用地表位移作为基准进行归一化； 2: 使用绝对值最大值作为基准进行归一化
        默认为0
    :param ifnegative: 是否对本征函数值进行反转，即加一个负号；
    :param mode_type: string 本征函数振型类型，尽管在路径中已经指出不同的振型类型存放于不同的文件夹下，但是依然需要指定
        振型类型以方便创建需要读取的文件的名称。分为U、V、W、P三种类型；Radial_output输出仍为U、V分量。
    :param linecolorsytle: string 绘制的本征函数的曲线颜色和线型，默认'k-'
    :param linewidth: float 绘制的本征函数的曲线宽度，默认为1
    :param if_radius_normalized: 输出的图像坐标轴中半径是否归一化；归一化后半径范围为相对值0-1。
    :param ifboundary: 是否在地表处显示一条横线表示在边界面
    :param ifaxis: 是否显示地球660间断面和CMB。
    :param ifreverse: False: 半径为y轴；True：半径为x轴。
    :return: 只进行plot，需要在函数外围设置figure和plt.show()
    '''

    if mode_type == 'S' or mode_type == 'U' or mode_type == 'V' or mode_type == 'P':
        Radius, U_eigenfunction, V_eigenfunction, P_eigenfunction = mineos_file_read(MINEOS_eigenfunction_dir=wdir, Modetype='S',
                     radial_order=radial_order, angular_order=angular_order, Mineos_file=[])
        if mode_type == 'U':
            weigenfunction = U_eigenfunction
        elif mode_type == 'V':
            weigenfunction = V_eigenfunction
        elif mode_type == 'P':
            weigenfunction = P_eigenfunction
    elif mode_type=='T' or mode_type == 'W':
        Radius, W_eigenfunction = mineos_file_read(MINEOS_eigenfunction_dir=wdir, Modetype='T',
                         radial_order=radial_order, angular_order=angular_order, Mineos_file=[])
        weigenfunction = W_eigenfunction
    elif mode_type=='R' or mode_type == 'r':
        Radius, R_eigenfunction = mineos_file_read(MINEOS_eigenfunction_dir=wdir, Modetype='R',
                         radial_order=radial_order, angular_order=angular_order, Mineos_file=[])
        weigenfunction = R_eigenfunction
    else:
        print('Wrong type of mode input.')

    wradius = Radius

    if Radius == []: # 确保不会出现错误
        print(f'NO such file of {radial_order}{mode_type}{angular_order}')
        return

    if ifnegative==True:
        weigenfunction = -np.array(weigenfunction)
    boundary670 = 670
    CMB = 3480
    SURFACE = 6371
    if if_radius_normalized == True:
        wradius = np.array(wradius) / max(np.array(wradius))
        boundary670 = 670 * 1000 / max(np.array(wradius))
        CMB = 3480 * 1000 / max(np.array(wradius))
        SURFACE = 6371 * 1000 / max(np.array(wradius))
    else:
        wradius = np.array(wradius) / 1000

    if ifnormalized == 1:  # 以地表位移作为条件归一化
        if ifreverse == False:
            plt.plot(np.array(weigenfunction) / weigenfunction[len(weigenfunction) - 1], wradius,
                         linecolorsytle, linewidth=linewidth)
        else:
            plt.plot(wradius, np.array(weigenfunction) / weigenfunction[len(weigenfunction) - 1],
                         linecolorsytle, linewidth=linewidth)

        if ifboundary == True:
            if weigenfunction[0] < 0:
                plt.axhline(y=wradius[0], xmin=weigenfunction[0] / weigenfunction[len(weigenfunction)-1], xmax=0, color='k',
                            linestyle='-', linewidth=1, alpha=1)
            else:
                plt.axhline(y=wradius[0], xmin=0, xmax=weigenfunction[0] / weigenfunction[len(weigenfunction)-1], color='k',
                            linestyle='-', linewidth=1, alpha=1)
            if weigenfunction[len(weigenfunction) - 1] < 0:
                plt.axhline(y=wradius[len(wradius) - 1], xmin=weigenfunction[len(weigenfunction) - 1] / max(weigenfunction),
                            xmax=0,
                            color='k', linestyle='-', linewidth=1, alpha=1)
            else:
                plt.axhline(y=wradius[len(wradius) - 1], xmin=0,
                            xmax=weigenfunction[len(weigenfunction) - 1] / weigenfunction[len(weigenfunction)-1],
                            color='k', linestyle='-', linewidth=1, alpha=1)
        if ifaxis == True:
            plt.text(-1, CMB, 'CMB')
            plt.text(-1, max(wradius)-boundary670, '670')
            plt.plot([-1, -0.7], [max(wradius) - boundary670, max(wradius) - boundary670], 'k-')
            plt.text(-1, SURFACE, '0')

    elif ifnormalized == 2:  # 当使用绝对值最大值作为归一化条件时
        if ifreverse == False:
            plt.plot(np.array(weigenfunction)/max(abs(np.array(weigenfunction))), wradius, linecolorsytle, linewidth=linewidth)
        else:
            plt.plot(wradius, np.array(weigenfunction) / max(abs(np.array(weigenfunction))),
                         linecolorsytle, linewidth=linewidth)

        if ifboundary == True:
            if weigenfunction[0] < 0:
                # plt.axhline(y=wradius[0]/1000, xmin=weigenfunction[0]/max(abs(np.array(weigenfunction))), xmax=0, color='k',
                #            linestyle='-', linewidth=1, alpha=1)
                plt.plot([weigenfunction[0] / max(abs(np.array(weigenfunction))), 0],
                         [wradius[0], wradius[0]], 'k-')
            else:
                # plt.axhline(y=wradius[0]/1000, xmin=0, xmax=weigenfunction[0] / max(abs(np.array(weigenfunction))), color='k',
                #            linestyle='-', linewidth=1, alpha=1)
                plt.plot([0, weigenfunction[0] / max(abs(np.array(weigenfunction)))],
                         [wradius[0], wradius[0]], 'k-')

            if weigenfunction[len(weigenfunction) - 1] < 0:
                # plt.axhline(y=wradius[len(wradius)-1]/1000, xmin=weigenfunction[len(weigenfunction) - 1] / max(abs(np.array(weigenfunction))),
                #             xmax=0, color='k', linestyle='-', linewidth=1, alpha=1)
                plt.plot([weigenfunction[len(weigenfunction) - 1] / max(abs(np.array(weigenfunction))), 0],
                         [wradius[len(wradius)-1], wradius[len(wradius)-1]], 'k-')
            else:
                # plt.axhline(y=wradius[len(wradius)-1]/1000, xmin=0,
                 #            xmax=weigenfunction[len(weigenfunction) - 1] / max(abs(np.array(weigenfunction))),
                 #            color='k', linestyle='-', linewidth=1, alpha=1)\
                plt.plot([0, weigenfunction[len(weigenfunction) - 1]]/ max(abs(np.array(weigenfunction))),
                         [wradius[len(wradius)-1], wradius[len(wradius)-1]], 'k-')
        if ifaxis == True:
            plt.text(-1, CMB, 'CMB')
            plt.text(-1, max(wradius) - boundary670, '670')
            plt.plot([-1,-0.7], [max(wradius) - boundary670, max(wradius) - boundary670], 'k-')
            plt.text(-1, SURFACE, '0')

    else: # 不进行归一化
        if ifreverse == False:
                plt.plot(np.array(weigenfunction), wradius, linecolorsytle, linewidth=linewidth)
        else:
                plt.plot(wradius, np.array(weigenfunction),
                     linecolorsytle, linewidth=linewidth)

        if ifboundary == True:
            if weigenfunction[0] < 0:
                plt.axhline(y=wradius[0], xmin=weigenfunction[0],
                            xmax=0, color='k',
                            linestyle='-', linewidth=1, alpha=1)
            else:
                plt.axhline(y=wradius[0], xmin=0,
                            xmax=weigenfunction[0], color='k',
                            linestyle='-', linewidth=1, alpha=1)
            if weigenfunction[len(weigenfunction) - 1] < 0:
                plt.axhline(y=wradius[len(wradius) - 1],
                            xmin=weigenfunction[len(weigenfunction) - 1],
                            xmax=0,
                            color='k', linestyle='-', linewidth=1, alpha=1)
            else:
                plt.axhline(y=wradius[len(wradius) - 1] / 1000, xmin=0,
                            xmax=weigenfunction[len(weigenfunction) - 1],
                            color='k', linestyle='-', linewidth=1, alpha=1)
        if ifaxis == True:
            plt.text(-1.2 * abs(min(weigenfunction)), CMB, 'CMB')
            plt.text(-1.2 * abs(min(weigenfunction)), max(wradius) - boundary670, '670')
            plt.plot([-1.2 * abs(min(weigenfunction)), -1 * abs(min(weigenfunction))], [max(wradius) - boundary670, max(wradius) - boundary670], 'k-')
            plt.text(-1.2 * abs(min(weigenfunction)), SURFACE, '0')