SPECNM 使用方法

Posted on | In

本文简介Python包 SPECNM使用方法

参考：https://gitlab.com/JohKem1/specnm

安装

使用anaconda安装。

配置conda源

官方：Please also setup the conda-forge channel to allow for installation of all required packages:

1
2
conda config --add channels conda-forge 
conda config --set channel_priority strict

选择一个空目录，复制specnm包

1
git clone https://gitlab.com/JohKem1/specnm.git

选择specnm所在目录

1
cd specnm

创建一个anaconda环境，用于安装specnm。可以使用已有的环境，但是不确定其他包是否兼容。

1
2
conda env create -f specnm.yml
conda activate specnm

此时可能会出现依赖包 percache 不存在的情况，可能是conda源不存在percache，因此可以手动使用pip安装：

1
pip install percache

如果有其他包不存在也可以用此方式手动安装。

安装specnm；进行conda develop。

1
conda develop .

运行测试

1
pytest

若无法通过测试，可以尝试：

git下载的specnm文件夹中还有一个名为specnm的文件夹，cd进入这个文件夹后能够看到love.py等函数文件，先进入这个文件夹再进行conda develop

1
2
3
cd specnm
conda develop .
pytest

计算环型振型相关信息

注意：无法限定radial order（径向阶数），只能使用angular order（角度阶数）和频率范围限定计算的振型边界范围。在计算PREM模型时，由于缺少0T1振型，因此n(radial order)的赋值存在偏差，真正的n值=赋值数+1。下文示例中存在该问题。

2024.5.28 修改了振型n值创建方式，去除了理论上不存在的0T1振型。

1
2
3
from specnm import love
ray = love(model=model_file, fmax=fmax)
ray_out = ray.love_problem(attenuation_mode='elastic', lmax=lmax)

ray_out为存储相关振型信息的字典，其中：

  • angular orders
  • wave numbers
  • frequencies
  • angular frequencies ：2pifrequencies
  • eigenfunctions ：为一个二维列表，每一个列表为对应频率振型的本征函数，即ray_out[‘eigenfunctions’][0]为ray_out[‘angular orders’]列表中第一个振型的本征函数W分量。
  • eigenfunctions_dr ：本征函数的导数
  • radius eigenfunctions对应的半径列表
  • phase velocities
  • group velocities
  • gamma
  • Q
  • epsilons
  • polynomial order 多项式

输出本征函数和振型信息的函数示例：
Love_output.py

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
from specnm import love
import numpy as np
from matplotlib import pyplot as plt

def Love_output(output_dir, model_dir, model, lmax, fmax, boundaries=[], l0=1, outNmax=100, attenuation_mode='elastic', needplot=True):    # 定义不同的模型文件类型
    '''
    !!!!!! 注意：需要提前在output_dir下创建文件夹：eigenfunctions、eigenfunctions_dr
    :param output_dir:  输出文件路径
    :param model_dir: 输入模型路径
    :param model: 输入模型名称
    :param l0: angular order最小值，默认为1
    :param lmax: angular order最大值
    :param fmax: 频率最大值
    :param boundaries: 是否需要绘制边界，如CMB等。例如需要，请输入：boundaries= [['ICB', 'CMB', 'Crust'], [1221.5, 3480., 6291.]]
    :param outNmax: 输出本征函数文件的径向阶数n的最大值范围
    :param attenuation_mode: love.love_problem中的衰减模式，'elastic' 'first order'
    :param needplot: 是否需要绘制n=0，l=2-10的频散曲线
    :return: 无。输出为本征函数文件、对应的Radius文件、振型相关信息文件。
    '''
    if model[-4:] == '.txt':
        output_file = model[0:-4] + '_T_output.txt'
        output_file2 = model[0:-4] + '_T_output_onlyfreq.txt'
    else:
        output_file = model + '_T_output.txt'
        output_file2 = model[0:-4] + '_T_output_onlyfreq.txt'


    model_file = model_dir + model
    output_dir1 = output_dir + 'eigenfunctions/'  # 输出本征函数
    output_dir2 = output_dir + 'eigenfunctions_dr/'  # 输出本征函数的导数

    ray = love(model=model_file, fmax=fmax)
    ray_out = ray.love_problem(attenuation_mode=attenuation_mode, lmax=lmax, l0=l0)

    angular_orders = ray_out['angular orders']
    # wave_numbers = ray_out['wave numbers']
    frequencies = ray_out['frequencies']*1000  # 单位转换成mhz
    periods = 1./ray_out['frequencies']  # 周期
    # angular_frequencies = ray_out['angular frequencies']
    eigenfunctions = ray_out['eigenfunctions']
    eigenfunctions_dr = ray_out['eigenfunctions_dr']
    phase_velocities = ray_out['phase velocities']
    group_velocities = ray_out['group velocities']
    gamma = ray_out['gamma']
    Q = ray_out['Q']
    epsilons = ray_out['epsilons']
    radius = ray_out['radius']
    polynomial_order = ray_out['polynomial order']

    radius_file_name = output_dir1 + 'Radius.txt'
    with open(radius_file_name, 'w') as file:
        radius = np.array(radius)
        for i in range(0, len(radius)):
            file.write(str(radius[i]))
            file.write('\n')

    for i in range(1, lmax+1):
        N_i = 0  # 当前angular order的索引
        N = 0  # 总数索引
        eigen_funi = []  # 当前angular order的eigenfunction
        eigen_dr = []  # 当前angular order的eigenfunction_dr

        # 找出当前angular order 的本征函数
        for q in angular_orders:
            if q == i:
                N_i = N_i + 1
                eigen_funi.append(eigenfunctions[N])
                eigen_dr.append(eigenfunctions_dr[N])
            N = N + 1

    # j 表示radial order
        for j in range(0, len(eigen_funi)):
            # 输出文件eigenfunction
            output_file_name = output_dir1 + 'W_' + f'radord{j}_angord{i}.txt'
            with open(output_file_name, 'w') as file:
                # 输出eigenfunction
                eigen_funi = np.array(eigen_funi)
                for m in range(0, len(eigen_funi[j])):
                    # file.write(str(eigen_funi[j][m]) + ' ' + str(eigen_dr[j][m]))
                    file.write(str(eigen_funi[j][m]))
                    file.write('\n')

        # 输出eigenfunction_dr
            output_file_name2 = output_dir2 + 'dr_W_' + f'radord{j}_angord{i}.txt'
            if j <= outNmax :
                if i == 1:  # 没有0T1振型，因此将n右移
                    j = j + 1
                with open(output_file_name2, 'w') as file:
                    # 输出radial order
                    # 输出eigenfunction
                    eigen_dr = np.array(eigen_dr)
                    for m in range(0, len(eigen_dr[j])):
                        # file.write(str(eigen_funi[j][m]) + ' ' + str(eigen_dr[j][m]))
                        file.write(str(eigen_dr[j][m]))
                        file.write('\n')

    # plt.show()
    if needplot==True:
        # 绘制0Wl 位移曲线
        from plot_eigenfunctions import open_file_and_plot_eigenfunctions_signal
        fig = plt.figure(figsize=(8, 6), dpi=600)
        ax_dr = fig.add_subplot(111)

        for l in range(2, 11):
            open_file_and_plot_eigenfunctions_signal(radial_order=0, angular_order=l, eigenfunction_type='W',
                                                     data_dir=output_dir1, ifnormalized=True)
        # 绘制网格，便于识别
        plt.plot([0, max(radius)/1000], [0.3, 0.3], 'k-', linewidth=0.5)
        plt.plot([0, max(radius) / 1000], [0.2, 0.2], 'b--', alpha=0.5, linewidth=0.5)  # alpha透明度 linewidth
        plt.plot([0, max(radius) / 1000], [0.4, 0.4], 'b--', alpha=0.5, linewidth=0.5)  # alpha透明度 linewidth
        plt.plot([0, max(radius) / 1000], [0.6, 0.6], 'b--', alpha=0.5, linewidth=0.5)  # alpha透明度 linewidth
        plt.plot([0, max(radius) / 1000], [0.8, 0.8], 'b--', alpha=0.5, linewidth=0.5)  # alpha透明度 linewidth
        plt.plot([0, max(radius) / 1000], [1.0, 1.0], 'b--', alpha=0.5, linewidth=0.5)  # alpha透明度 linewidth

        # 获取当前ylim
        current_ylim = ax_dr.get_ylim()

        # 绘制模型边界
        if boundaries == []:
            mmmmm=1
        else:
            for boundare_name, boundare_radius in zip(boundaries[0], boundaries[1]):
                ax_dr.axvline(x=boundare_radius, color='r', linestyle='--', linewidth=0.3, alpha=0.7)
                # plt.plot([boundare_radius, boundare_radius], [0, 1], 'r--', linewidth=0.3, alpha=0.7)
                plt.text(boundare_radius, current_ylim[0], boundare_name)
        ax_dr.set_ylim(current_ylim)
        #########################
        plt.xlabel("Radius(km)")
        plt.ylabel("0Wl")
        plt.xlim([0, max(radius)/1000])
        plt.title(f"Eigenfunctions of {model} model")

        figname = output_dir + 'W_eigenfunctions.png'
        plt.savefig(figname)
        plt.close(fig)

    # plt.show()

    # 重整mode相关参数，进行输出
    mode_out_put = []
    for l in range(1, lmax+1):
        l_dictionary = {'n': [], 'frequency': [], 'phase_velocity': [], 'group_velocity': [], 'gamma': [],
                        'Q': [], 'epsilons': [], 'periods': []}
        if l == 1:  # 没有0T1振型，因此去除该情况
            n = 1
        else:
            n = 0
        for i in range(0, len(angular_orders)):
            if angular_orders[i] == l:
                l_dictionary['n'].append(n)
                n = n + 1
                l_dictionary['frequency'].append(frequencies[i])
                l_dictionary['phase_velocity'].append(phase_velocities[i])
                l_dictionary['group_velocity'].append(group_velocities[i])
                l_dictionary['gamma'].append(gamma[i])
                l_dictionary['Q'].append(Q[i])
                l_dictionary['epsilons'].append(epsilons[i])
                l_dictionary['periods'].append(periods[i])
        mode_out_put.append(l_dictionary)


# 输出mode参数到文件。
    output_file = output_dir + output_file
    with open(output_file, 'w') as file:
        file.write(f'MODEL is {model}, polynomial order is {polynomial_order} \n')
        # file.write('n l freq(mHz) period(s) ph_vel grp_vel gamma Q epsilon \n')
        file.write(
            '{:>6s}{:>6s}{:>15s}{:>12s}{:>12s}{:>12s}{:>12s}{:>12s}{:>12s}'.format(
                'n',
                'l',
                'freq(mHz)',
                'period(s)',
                'ph_v(m/s)',
                'grp_v(m/s)',
                'gamma',
                'Q',
                'epsilon'
            ))
        file.write('\n')
        for l in range(1, lmax+1):
            # n = mode_out_put[l-1]['n']
            for pp in range(0, len(mode_out_put[l-1]['n'])):
                file.write(
                    '{:>6d}{:>6d}{:>15.5f}{:>12.3f}{:>12.3f}{:>12.3f}{:>12.3f}{:>12.3f}{:>12.3f}'.format(
                        mode_out_put[l - 1]['n'][pp],
                        l,
                        mode_out_put[l - 1]['frequency'][pp],
                        mode_out_put[l - 1]['periods'][pp],
                        mode_out_put[l - 1]['phase_velocity'][pp],
                        mode_out_put[l - 1]['group_velocity'][pp],
                        mode_out_put[l - 1]['gamma'][pp],
                        mode_out_put[l - 1]['Q'][pp],
                        mode_out_put[l - 1]['epsilons'][pp]
                    ))
                file.write('\n')


# 只输出mode参数到文件。
    output_file2 = output_dir + output_file2
    with open(output_file2, 'w') as file:
        file.write(f'MODEL is {model}, polynomial order is {polynomial_order} \n')
        # file.write('n l freq(mHz) period(s) \n')
        file.write(
            '{:>6s}{:>6s}{:>15s}{:>12s}'.format(
               'n',
                'l',
                'freq(mHz)',
                'period(s)'
            ))
        file.write('\n')
        for l in range(1, lmax+1):
            # n = mode_out_put[l-1]['n']
            for pp in range(0, len(mode_out_put[l-1]['n'])):
                file.write(
                    '{:>6d}{:>6d}{:>15.5f}{:>12.3f}'.format(
                        mode_out_put[l - 1]['n'][pp],
                        l,
                        mode_out_put[l - 1]['frequency'][pp],
                        mode_out_put[l - 1]['periods'][pp]
                    ))
                file.write('\n')

计算球型振型(l=0)相关信息

由于l=0时称为radial mode，因此单独计算。相比于Love，其没有lmax选项；没有相速度和群速度输出。

1
2
3
from specnm import rayleigh
ray = radial(model=model_file, fmax=fmax)
rayout = ray.radial_problem(self, fmin=0.0001388888888888889, fmax=None, attenuation_mode=None, **multishift_kwargs)

radial_problem函数只有两个参数：fmin和fmax

  • fmin 浮点型 最小输出频率
  • fmax 浮点型 最大输出频率
  • attentuation_mode 字符串格式 默认为None， 可以输入’elastic’。
  • multishift_kwargs 字典格式 目前未知，官方介绍：additional arguments for the multishift

输出本征函数和振型信息的函数示例：
Rayleigh_output.py

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
from specnm import radial
import numpy as np
from matplotlib import pyplot as plt

def Radial_output(output_dir, model_dir, model, fmax, boundaries=[], outNmax=100, attenuation_mode='elastic'):    # 定义不同的模型文件类型
    '''
    !!!!!! 注意：需要提前在output_dir下创建文件夹：eigenfunctions、eigenfunctions_dr
    :param output_dir:  输出文件路径
    :param model_dir: 输入模型路径
    :param model: 输入模型名称
    :param lmax: angular order最大值
    :param fmax: 频率最大值
    :param boundaries: 是否需要绘制边界，如CMB等。例如需要，请输入：boundaries= [['ICB', 'CMB', 'Crust'], [1221.5, 3480., 6291.]]
    :param outNmax: 输出本征函数文件的径向阶数n的最大值范围
    :return: 无。输出为本征函数文件、对应的Radius文件、振型相关信息文件。
    '''

    if model[-4:] == '.txt':
        output_file = model[0:-4] + '_R_output.txt'
        output_file2 = model[0:-4] + '_R_output_onlyfreq.txt'
    else:
        output_file = model + '_R_output.txt'
        output_file2 = model + '_R_output_onlyfreq.txt'

    model_file = model_dir + model
    output_dir1 = output_dir + 'eigenfunctions/'  # 输出本征函数
    output_dir2 = output_dir + 'eigenfunctions_dr/'  # 输出本征函数的导数


    ray = radial(model=model_file, fmax=fmax)
    ray_out = ray.radial_problem(attenuation_mode=attenuation_mode)
    print(ray_out)
    angular_orders = ray_out['angular orders']
    # wave_numbers = ray_out['wave numbers']
    frequencies = ray_out['frequencies']*1000  # 单位转换成mhz
    periods = 1./ray_out['frequencies']  # 周期
    # angular_frequencies = ray_out['angular frequencies']
    eigenfunctions = ray_out['eigenfunctions']
    eigenfunctions_dr = ray_out['eigenfunctions_dr']
    phase_velocities = ray_out['phase velocities']
    group_velocities = ray_out['group velocities']
    gamma = ray_out['gamma']
    Q = ray_out['Q']
    epsilons = ray_out['epsilons']
    radius = ray_out['radius']
    polynomial_order = ray_out['polynomial order']

    radius_file_name = output_dir1 + 'Radius.txt'
    with open(radius_file_name, 'w') as file:
        radius = np.array(radius)
        for i in range(0, len(radius)):
            file.write(str(radius[i]))
            file.write('\n')

    # j 表示radial order
    for j in range(0, len(angular_orders)):
        if j <= outNmax:
# 输出文件eigenfunction
            output_file_name = output_dir1 + 'U_' + f'radord{j}_angord0.txt'
            with open(output_file_name, 'w') as file:
                # 输出eigenfunction
                for m in range(0, len(eigenfunctions[j])):
                    # file.write(str(eigen_funi[j][m]) + ' ' + str(eigen_dr[j][m]))
                    file.write(str(eigenfunctions[j][m]))
                    file.write('\n')

            # 输出eigenfunction_dr
            output_file_name2 = output_dir2 + 'dr_U_' + f'radord{j}_angord{i}.txt'
            with open(output_file_name2, 'w') as file:
                # 输出eigenfunction
                for m in range(0, len(eigenfunctions_dr[j])):
                    # file.write(str(eigen_funi[j][m]) + ' ' + str(eigen_dr[j][m]))
                    file.write(str(eigenfunctions_dr[j][m]))
                    file.write('\n')

# 输出mode参数到文件。
    output_file = output_dir + output_file
    with open(output_file, 'w') as file:
        file.write(f'MODEL is {model}, polynomial order is {polynomial_order} \n')
        # file.write('n l freq(mHz) period(s)  gamma Q epsilon \n')
        file.write(
            '{:>6s}{:>6s}{:>15s}{:>12s}{:>12s}{:>12s}{:>12s}'.format(
                'n',
                'l',
                'freq(mHz)',
                'period(s)',
                'gamma',
                'Q',
                'epsilons'
            ))
        file.write('\n')
        for i in range(0, len(angular_orders)):
            n = i
            '''
            file.write(str(n) + ' ' + '0'
                       + ' ' + str(frequencies[n])
                       + ' ' + str(periods[n])
                       + ' '  + str(gamma[n]) + ' ' + str(Q[n]) + ' '
                       + str(epsilons[n]))
            '''
            file.write(
                '{:>6d}{:>6d}{:>15.5f}{:>12.3f}{:>12.3f}{:>12.3f}{:>12.3f}'.format(
                    n,
                    0,
                    frequencies[n],
                    periods[n],
                    gamma[n],
                    Q[n],
                    epsilons[n]
                ))
            file.write('\n')


# 只输出mode参数到文件。
    output_file2 = output_dir + output_file2
    with open(output_file2, 'w') as file:
        file.write(f'MODEL is {model}, polynomial order is {polynomial_order} \n')
        # file.write('n l freq(mHz) period(s)  gamma Q epsilon \n')
        file.write(
            '{:>6s}{:>6s}{:>15s}{:>12s}'.format(
                'n',
                'l',
                'freq(mHz)',
                'period(s)'
            ))
        file.write('\n')
        for i in range(0, len(angular_orders)):
            n = i
            '''
            file.write(str(n) + ' ' + '0'
                       + ' ' + str(frequencies[n])
                       + ' ' + str(periods[n])
                       + ' '  + str(gamma[n]) + ' ' + str(Q[n]) + ' '
                       + str(epsilons[n]))
            '''
            file.write(
                '{:>6d}{:>6d}{:>15.5f}{:>12.3f}'.format(
                    n,
                    0,
                    frequencies[n],
                    periods[n]
                ))
            file.write('\n')

计算球型振型相关信息

有两个函数：

  • rayleigh 不考虑重力影响，本征函数具有U、V分量。
  • rayleigh_fg 计算具有full gravity的球型振型，本征函数具有U、V、P分量。

rayleigh

1
2
3
from specnm import rayleigh
ray = rayleigh(model=model_file, fmax=fmax)
ray_out = ray.rayleigh_problem(attenuation_mode='elastic',lmax=lmax)

ray_out为存储相关振型信息的字典，其中：

  • angular orders
  • wave numbers
  • frequencies
  • angular frequencies ：2pifrequencies
  • eigenfunctions ：为一个二维列表，每一个列表为对应频率振型的本征函数，即ray_out[‘eigenfunctions’][0]为ray_out[‘angular orders’]列表中第一个振型的本征函数的两个分量。其中，两个分量按照顺序排列索引分别为：0、2、4……为U分量，1、3、5……为V分量。
    • U_eigen_funi = eigenfunctions[i][0::2][0:] 获取U分量
    • V_eigen_funi = eigenfunctions[i][1::2][0:] 获取V分量
      注意：官方test中去除了第一个数值（radius=0）的本征函数数值，但是radius是从0开始，因此不知其去除的含义。如下方式获取某一振型两个分量的本征函数。i为angular orders中的振型索引，[0::2]表示从第0个开始至末尾获取数据，步长为2；[0:]表示对获取到的本征函数进行切片，下文中未进行切片，但是官方test中去除了第一个数据，因此使用[1:]，未知原因。
  • eigenfunctions_dr ：本征函数的导数
  • radius eigenfunctions对应的半径列表
  • phase velocities
  • group velocities
  • gamma
  • Q
  • epsilons
  • polynomial order 多项式

输出本征函数和振型信息的函数示例：

修改日志：
2024.6.18 修改了地球模型需要去除不存在的${1}$S${1}$振型的问题
Rayleigh_output.py

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
from specnm import love
from specnm import seismogram
from specnm import spheroidal_base
from specnm import specnm_base

import sys

import warnings
import inspect
import os
import pathlib
import numpy as np
from specnm import rayleigh_fg
from specnm import rayleigh
import matplotlib.pyplot as plt

# 定义不同的模型文件类型
def Ray_output(output_dir, model_dir, model, lmax, fmax, boundaries, l0=1, outNmax=100, attenuation_mode='elastic',
               epsilon_spurious=1e-3, needplot=True, planet_type='Earth'):
    '''
     !!!!!! 注意：需要提前在output_dir下创建文件夹：eigenfunctions、eigenfunctions_dr
     :param output_dir:  输出文件路径
     :param model_dir: 输入模型路径
     :param model: 输入模型名称
     :param l0: angular order最小值，默认为1
     :param lmax: angular order最大值
     :param fmax: 频率最大值
     :param boundaries: 是否需要绘制边界，如CMB等。例如需要，请输入：boundaries= [['ICB', 'CMB', 'Crust'], [1221.5, 3480., 6291.]]
     :param outNmax: 输出本征函数文件的径向阶数n的最大值范围
     :param attenuation_mode: love.love_problem中的衰减模式，'elastic' 'first order'
     :param needplot: 是否需要绘制n=0，l=2-10的频散曲线
     :param planet_type: 默认为'Earth'，因为地球模型需要去除0S1、1S1振型，而其他模型不确定是否需要去除1S1,目前保留这一问题
     :return: 无。输出为本征函数文件、对应的Radius文件、振型相关信息文件。
     '''

    if model[-4:] == '.txt':
        output_file = model[0:-4] + '_S_output.txt'
        output_file2 = model[0:-4] + '_S_output_onlyfreq.txt'
    else:
        output_file = model + '_S_output.txt'
        output_file2 = model + '_S_output_onlyfreq.txt'

    model_file = model_dir + model
    output_dir1 = output_dir + 'eigenfunctions/'  # 输出本征函数
    output_dir2 = output_dir + 'eigenfunctions_dr/'  # 输出本征函数的导数
    # 计算rayleigh模型相关信息
    ray = rayleigh(model=model_file, fmax=fmax)
    ray_out = ray.rayleigh_problem(attenuation_mode=attenuation_mode,lmax=lmax,l0=l0, epsilon_spurious=epsilon_spurious)

    angular_orders = ray_out['angular orders']
    wave_numbers = ray_out['wave numbers']
    frequencies = ray_out['frequencies'] * 1000   # 单位转换成mhz
    periods = 1./ray_out['frequencies']   # 周期
    angular_frequencies = ray_out['angular frequencies']
    eigenfunctions = ray_out['eigenfunctions']
    eigenfunctions_dr = ray_out['eigenfunctions_dr']
    phase_velocities = ray_out['phase velocities']
    group_velocities = ray_out['group velocities']
    gamma = ray_out['gamma']
    Q = ray_out['Q']
    epsilons = ray_out['epsilons']
    radius = ray_out['radius']
    polynomial_order = ray_out['polynomial order']

    # 输出本征函数对应的半径信息
    radius_file_name = output_dir1 + 'Radius.txt'
    with open(radius_file_name, 'w') as file:
        radius = np.array(radius)
        for i in range(0, len(radius)):
            file.write(str(radius[i]))
            file.write('\n')

    # 输出本征函数相关内容
    j_r = 0
    for i in range(1, lmax+1):
        N_i = 0  # 当前angular order的索引
        N = 0  # 总数索引
        eigen_funi = []  # 当前angular order的eigenfunction
        eigen_dr = []  # 当前angular order的eigenfunction_dr


        # 找出当前angular order 的本征函数
        for q in angular_orders:
            if q == i:
                N_i = N_i + 1
                eigen_funi.append(eigenfunctions[N])
                eigen_dr.append(eigenfunctions_dr[N])
            N = N + 1

    # j 表示radial order
        for j in range(0, len(eigen_funi)):
        # 输出文件eigenfunction
        # 绘制图像
            if j == 0:
                # print(f'j={j}')
                j_r = j_r+1
                eigen_funi = np.array(eigen_funi) # eigen_funi是一个包含不同radial order，相同angular order的本征函数组。
                # rayleigh中包含u v两个分量，因此需要分离两个分量。
    #   官方test中去除了第一个数值（radius=0）的本征函数数值，但是radius是从0开始，因此不知其去除的含义。
            U_eigen_funi = eigen_funi[j, 0::2][0:]
            V_eigen_funi = eigen_funi[j, 1::2][0:]

            if j <= outNmax:
                if i == 1:  # 没有0S1振型，因此将n右移;地球模型不存在1S1
                    if planet_type == 'Earth':
                        j=j+2
                    else:
                        j = j + 1
                output_file_nameU = output_dir1 + 'U_' + f'radord{j}_angord{i}.txt'
                with open(output_file_nameU, 'w') as file:
                    for m in range(0, len(U_eigen_funi)):
                    # file.write(str(eigen_funi[j][m]) + ' ' + str(eigen_dr[j][m]))
                        file.write(str(U_eigen_funi[m]))
                        file.write('\n')
                output_file_nameV = output_dir1 + 'V_' + f'radord{j}_angord{i}.txt'
                with open(output_file_nameV, 'w') as file:
                    for m in range(0, len(V_eigen_funi)):
                    # file.write(str(eigen_funi[j][m]) + ' ' + str(eigen_dr[j][m]))
                        file.write(str(V_eigen_funi[m]))
                        file.write('\n')
        # 输出radial order
        # 输出eigenfunction

    print(f'jr={j_r}')

    if needplot==True:
        # 绘制本征函数曲线
        from plot_eigenfunctions import open_file_and_plot_eigenfunctions_signal
        fig = plt.figure(figsize=(8, 6), dpi=600)
        # 绘制0Ul 位移曲线
        ax_dr = fig.add_subplot(111)
        for l in range(2, 11):
            open_file_and_plot_eigenfunctions_signal(radial_order=0, angular_order=l, eigenfunction_type='U', data_dir=output_dir1, ifnormalized=True)
        # 绘制网格，便于识别
        plt.plot([0, max(radius) / 1000], [0.3, 0.3], 'k-', linewidth=0.5)
        plt.plot([0, max(radius) / 1000], [0.2, 0.2], 'b--', alpha=0.5, linewidth=0.5)  # alpha透明度 linewidth
        plt.plot([0, max(radius) / 1000], [0.4, 0.4], 'b--', alpha=0.5, linewidth=0.5)  # alpha透明度 linewidth
        plt.plot([0, max(radius) / 1000], [0.6, 0.6], 'b--', alpha=0.5, linewidth=0.5)  # alpha透明度 linewidth
        plt.plot([0, max(radius) / 1000], [0.8, 0.8], 'b--', alpha=0.5, linewidth=0.5)  # alpha透明度 linewidth
        plt.plot([0, max(radius) / 1000], [1.0, 1.0], 'b--', alpha=0.5, linewidth=0.5)  # alpha透明度 linewidth

        # 获取当前ylim
        current_ylim = ax_dr.get_ylim()
        # 绘制模型边界
        if boundaries == []:
            mmmmm=1
        else:
            for boundare_name, boundare_radius in zip(boundaries[0], boundaries[1]):
                ax_dr.axvline(x=boundare_radius, color='r', linestyle='--', linewidth=0.3, alpha=0.7)
                # plt.plot([boundare_radius, boundare_radius], [0, 1], 'r--', linewidth=0.3, alpha=0.7)
                plt.text(boundare_radius, current_ylim[0], boundare_name)
        ax_dr.set_ylim(current_ylim)
        #########################

        plt.plot([0, max(radius)/1000], [0.3, 0.3], 'k-')
        plt.xlabel("Radius(km)")
        plt.ylabel("0Ul")
        plt.xlim([0, max(radius)/1000])
        # plt.ylim([-0.2, 1.2])
        plt.title(f"Eigenfunctions of {model} model")
        figname = output_dir + 'U_eigenfunctions.png'
        plt.savefig(figname)
        plt.close()

        fig = plt.figure(figsize=(8, 6), dpi=600)
        # 绘制0Vl 位移曲线
        ax_dr = fig.add_subplot(111)
        for l in range(2, 11):
            open_file_and_plot_eigenfunctions_signal(radial_order=0, angular_order=l, eigenfunction_type='V', data_dir=output_dir1, ifnormalized=True)
        # 绘制网格，便于识别
        plt.plot([0, max(radius) / 1000], [0.3, 0.3], 'k-', linewidth=0.5)
        plt.plot([0, max(radius) / 1000], [0.2, 0.2], 'b--', alpha=0.5, linewidth=0.5)  # alpha透明度 linewidth
        plt.plot([0, max(radius) / 1000], [0.4, 0.4], 'b--', alpha=0.5, linewidth=0.5)  # alpha透明度 linewidth
        plt.plot([0, max(radius) / 1000], [0.6, 0.6], 'b--', alpha=0.5, linewidth=0.5)  # alpha透明度 linewidth
        plt.plot([0, max(radius) / 1000], [0.8, 0.8], 'b--', alpha=0.5, linewidth=0.5)  # alpha透明度 linewidth
        plt.plot([0, max(radius) / 1000], [1.0, 1.0], 'b--', alpha=0.5, linewidth=0.5)  # alpha透明度 linewidth

        # 获取当前ylim
        current_ylim = ax_dr.get_ylim()
        # 绘制模型边界
        if boundaries == []:
            mmmmm=1
        else:
            for boundare_name, boundare_radius in zip(boundaries[0], boundaries[1]):
                ax_dr.axvline(x=boundare_radius, color='r', linestyle='--', linewidth=0.3, alpha=0.7)
                # plt.plot([boundare_radius, boundare_radius], [0, 1], 'r--', linewidth=0.3, alpha=0.7)
                plt.text(boundare_radius, current_ylim[0], boundare_name)
        ax_dr.set_ylim(current_ylim)
        #########################
        plt.plot([0, max(radius)/1000], [0.3, 0.3], 'k-')
        plt.xlabel("Radius(km)")
        plt.ylabel("0Vl")
        plt.xlim([0, max(radius)/1000])
        # plt.ylim([0, 1.2])
        plt.title(f"Eigenfunctions of {model} model")
        figname = output_dir + 'V_eigenfunctions.png'
        plt.savefig(figname)
        plt.close(fig)

    # 重整mode相关参数，进行输出
    mode_out_put = []
    for l in range(1, lmax+1):
        l_dictionary = {'n': [], 'frequency': [], 'phase_velocity': [], 'group_velocity': [], 'gamma': [], 'Q': [], 'epsilons': [], 'periods': []}
        if l == 1:  # 没有0S1振型，因此去除该情况  地球模型没有1S1振型
            if planet_type == 'Earth':
                n = 2
            else:
                n = 1
        else:
            n = 0
        for i in range(0, len(angular_orders)):
            if angular_orders[i] == l:
                l_dictionary['n'].append(n)
                n = n + 1
                l_dictionary['frequency'].append(frequencies[i])
                l_dictionary['phase_velocity'].append(phase_velocities[i])
                l_dictionary['group_velocity'].append(group_velocities[i])
                l_dictionary['gamma'].append(gamma[i])
                l_dictionary['Q'].append(Q[i])
                l_dictionary['epsilons'].append(epsilons[i])
                l_dictionary['periods'].append(periods[i])
        mode_out_put.append(l_dictionary)


    # 输出mode参数到文件。
    output_file = output_dir + output_file
    with open(output_file, 'w') as file:
        file.write(f'MODEL is {model}, polynomial order is {polynomial_order} \n')
        file.write(
            '{:>6s}{:>6s}{:>15s}{:>12s}{:>12s}{:>12s}{:>12s}{:>12s}{:>12s}'.format(
                'n',
                'l',
                'freq(mHz)',
                'period(s)',
                'ph_v(m/s)',
                'grp_v(m/s)',
                'gamma',
                'Q',
                'epsilon'
            ))
        file.write('\n')
        for l in range(1, lmax+1):
            #　n = mode_out_put[l-1]['n']
            for pp in range(0, len(mode_out_put[l-1]['n'])):
                file.write(
                    '{:>6d}{:>6d}{:>15.5f}{:>12.3f}{:>12.3f}{:>12.3f}{:>12.7f}{:>12.3f}{:>12.9f}'.format(
                        mode_out_put[l - 1]['n'][pp],
                        l,
                        mode_out_put[l - 1]['frequency'][pp],
                        mode_out_put[l - 1]['periods'][pp],
                        mode_out_put[l - 1]['phase_velocity'][pp],
                        mode_out_put[l - 1]['group_velocity'][pp],
                        mode_out_put[l - 1]['gamma'][pp],
                        mode_out_put[l - 1]['Q'][pp],
                        mode_out_put[l - 1]['epsilons'][pp]
                    ))

                file.write('\n')

# 只输出mode参数到文件。
    output_file2 = output_dir + output_file2
    with open(output_file2, 'w') as file:
        file.write(f'MODEL is {model}, polynomial order is {polynomial_order} \n')
        # file.write('n l freq(mHz) period(s) \n')
        file.write(
            '{:>6s}{:>6s}{:>15s}{:>12s}'.format(
                'n',
                'l',
                'freq(mHz)',
                'period(s)'
            ))
        file.write('\n')
        for l in range(1, lmax+1):
            # n = mode_out_put[l-1]['n']
            for pp in range(0, len(mode_out_put[l-1]['n'])):
                file.write(
                    '{:>6d}{:>6d}{:>15.5f}{:>12.3f}'.format(
                        mode_out_put[l - 1]['n'][pp],
                        l,
                        mode_out_put[l - 1]['frequency'][pp],
                        mode_out_put[l - 1]['periods'][pp]
                    ))
                file.write('\n')

rayleigh_fg

同rayleigh，eigenfunctions的索引不同：
0、3、6……为U分量，1、4、7……为V分量，2、5、8为P分量。

使用上述函数输出振型信息

输出内容为：振型相关信息、每个振型的本征函数文件、基阶振型的三分量本征函数图，每个分量保存在对应的output_dir（输出路径）中：例如W分量保存在Love_output文件夹下。

先将上述函数保存到相应的py文件中，在目录下创建use.py文件，即可。
注意：提前在输出路径下创建好output_dir以及在output_dir下创建文件夹：eigenfunctions、eigenfunctions_dr。

use.py

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
from Love_output import Love_output
from Rayleigh_output import Ray_output

# 模型的边界信息，例如CMB等分界面
prem_noocean_boundary = [['ICB', 'CMB', 'Crust'], [1221.5, 3480., 6291.]]
MoonWeber_boundary = [['ICB', 'CMB', 'Crust'], [240., 330., 1697.1]]
MoonGarcia_boundary = [['CMB', 'Crust'], [380., 1709.1]]
model_dir = '/home/shizhiyuanubuntu/PycharmProjects/Seismology/models/'

# 需要修改的内容
'''
lmax ：angular order计算范围
fmax ：频率计算范围
model ：使用的模型名称（注意是否存在.txt）
boundary ：使用的模型边界信息，在上方给出
# model = 'MoonGarcia.txt'
# model = 'MoonWeber.txt'
output_dir ： 球型振型信息输出路径
output_dir2 ：环形振型信息输出路径
'''
lmax = 100
fmax = 1
model = 'prem_noocean'
boundary = prem_noocean_boundary
# model = 'MoonGarcia.txt'
# model = 'MoonWeber.txt'
output_dir = '/home/shizhiyuanubuntu/PycharmProjects/Seismology/prem_noocean/Rayleigh_output/'
output_dir2 = '/home/shizhiyuanubuntu/PycharmProjects/Seismology/prem_noocean/Love_output/'

# 输出球型振型相关信息和本征函数
Ray_output(output_dir=output_dir, model_dir=model_dir, model=model, lmax=lmax, fmax=fmax, boundaries=boundary)
# 输出环型振型相关信息和本征函数
Love_output(output_dir=output_dir2, model_dir=model_dir, model=model, lmax=lmax, fmax=fmax, boundaries=boundary)

绘制地震图

设置震源和台站信息

使用instaseis包中的Source、Receiver函数：

使用Source模块，使用矩张量信息定义震源参数。

示例：

1
2
3
4
5
6
7
8
9
10
from instaseis import Source, Receiver
import numpy as np

M = np.zeros(6)
M[0] = -1*1e26
M[1] = 1*1e26
# 设置震源
src = Source(latitude=0., longitude=0., depth_in_m=10e4, m_rr=M[0], m_tt=M[1], m_pp=M[2], m_rt=M[3], m_rp=M[4], m_tp=M[5])
# 设置接收台站
rec = Receiver(0., 90.)

创建或读取振型数据文件：

创建文件

使用seismogram模块中的mode_database模块的compute函数创建mode数据文件：

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
import os
import pathlib
from specnm import love
from specnm import seismogram
from specnm.seismogram import mode_database
# 模型所在的路径
MODEL_DIR = '/home/shizhiyuanubuntu/PycharmProjects/Seismology/models'

# 设置模型文件名称，即使用的模型。
model = 'noLVZ'
model_file = model + '.txt'

# 设置输出文件的名称，输出的内容为模型计算的文件，之后将使用它进行绘制地震图
output_file = model + '_output.txt'
# 输出文件的路径
tmpdir = '/home/shizhiyuanubuntu/PycharmProjects/Seismology/seismogram_make/data'

# 计算相关振型信息文件
mdb = mode_database.compute(
    model=os.path.join(MODEL_DIR, model_file),
    fmax=0.15,
    mode_types='RST',
    gravity_mode='full',
    init_kwargs={'verbose': True},
    problem_kwargs={'attenuation_mode': 'first order', 'lmax': 300}
)

# 输出相关振型信息文件
filename = os.path.join(tmpdir, f'{model}data.h5')
print('filetest')
print(mdb.modes)
mdb.write_h5(filename)

读取文件

如果已经计算好模型信息文件，可以直接读取，省去计算过程。可以使用mode_database模块的read()函数读取数据:

1
2
3
tmpdir = '/home/shizhiyuanubuntu/PycharmProjects/Seismology/seismogram_make/data'
filename = os.path.join(tmpdir, f'{model}data.h5')
mdb = mode_database.read(filename=filename)

获取地震图并保存在当前目录

使用mode_database模块创建的对象mdb中的属性get_seismogram函数创建地震图：

其中

1
get_seismograms(source, receiver, mode_types='TRS', dt=1., nt=3600, t0=None, kind='displacement', seismometer_response=False, flattening=0.)

函数中的参数：

  • source 上述instaseis设置的震源信息
  • receiver 上述instaseis设置的台站信息
  • mode_types 使用的振型类型，一般包括：
    • T 环型振型
    • S 球型振型
    • R l=0时的球型振型
  • dt 地震图时间间隔，采样率的倒数
  • nt 地震图采样点数
  • t0 地震图起始时间，默认为1970年
  • kind 地震图类型
    • displacement
    • velocity
    • acceleration
  • seismometer_response 是否加入仪器响应
  • flattening 存疑 ？展平变换
1
2
3
4
5
6
7
streams = mdb.get_seismograms(source=src, receiver=rec, mode_types='TRS', dt=1., nt=30000, t0=None, kind='velocity', seismometer_response=False, flattening=0.)
streams.plot()
trz = streams.select(channel='XZ')[0].write(f'{model}_{distance}_Z.SAC', format('SAC'))
tre = streams.select(channel='XE')[0].write(f'{model}_{distance}_E.SAC', format('SAC'))
trn = streams.select(channel='XN')[0].write(f'{model}_{distance}_N.SAC', format('SAC'))
trr = streams.select(channel='XR')[0].write(f'{model}_{distance}_R.SAC', format('SAC'))
trt = streams.select(channel='XT')[0].write(f'{model}_{distance}_T.SAC', format('SAC'))

streams存在属性select(channel=‘XZ’)
可以选择不同的通道。streams有五个通道，分别为：XZ、XN、XE、XR、XT。