使用Pyrocko Manual插件（python）计算P波到时（没做出来版本，中途改用taup了）

Pyrocko Manual是地震学开源工具箱
本人只想计算理论P波到时，所以其他功能也没来得及学习。没啥逻辑，当个笔记为自己做个记录。也不一定能成功。
%%%%%%%%%%%%%%
做了3天，鼻青脸肿的回来了，果然，啥资料没找到，尝试了无数次，还是失败了。但是也不完全算失败吧，就是找到了替代的方法，反正最后用别的软件实现了。这个软件呢，我只用明白了一半，中间一些代码我反复试验，感觉摸了一点点门道吧。也记录一下。
最后用这个生成的数据，改使用taup了。都写这里，对看客有帮助那最好了，如果没有，有做过的朋友看看我的笔记可以来指导我这个大龄傻子一下吗。

安装工具包

最基本的，在终端（Terminal）直接敲

pip install pyrocko

完成以后，使用该命令查看相关版本信息

>pip show pyrocko

计算P波到时，S波到时

这个工具箱对比obspy工具箱，我暂时觉得的好处就是，可视化更方便，而且可以改变速度模型（个人理解，因为这两个插件我都不太会用）。初步尝试感觉功能最强大的还是JAVA环境下的taup。但是这借个我都好陌生。
主要参考教程cake manual
pyrocko.cake

所有速度模型

先看看所有速度模型，在终端输入下面的

cake list-models

会常出现所有的速度模型，这里每个模型有三个文件，都是表示同一地球模型的不同文件格式，可以根据需求选择相应的文件来加载地球模型。

• .m：通常表示模型文件，其中包含地球内部结构的层和不连续性的定义。这些文件通常包含模型的详细参数，如速度、密度等。
• .f：通常表示速度模型文件，其中仅包含速度值，通常用于绘制速度剖面或进行速度分析。
• .l：可能表示模型的层文件，提供关于每个层的详细信息，如层厚度、材料参数等。
  用法如下

cake <subcommand> [options]

比如（不要怀疑，直接扔到终端）：是不是很耐斯，这么高端的图一下子就出来了，BIU!〜(￣▽￣〜) (〜￣▽￣)〜

cake plot-model --model=prem-no-ocean.m

自定义速度模型

你可以通过以下几种方式初始化此类的实例：

• 使用模块函数 load_model() 从文件中读取模型。
• 使用默认构造函数创建一个空模型，并使用 append() 方法（从顶部到底部）添加层和不连续性。
• 使用构造函数 LayeredModel.from_scanlines()，从给定的速度剖面自动创建层和不连续性对象。

主要参数说明

参数:

• vp: P 波速度 [m/s]
• vs: S 波速度 [m/s]
• rho: 密度 [kg/m^3]
• qp: P 波衰减 Qp
• qs: S 波衰减 Qs
• poisson: 泊松比
• lame: Lame 参数 lambda 和剪切模量的元组 [Pa]
  qk: 体波衰减 Qk
  qmu: 剪切波衰减 Qmu
  burgers（元组）：Burgers 流变学参数，包括瞬时粘度 [Pa]（小于等于 0 表示无穷大值）、稳态粘度 [Pa] 和 alpha，即有效和未释放剪切模量之比，mu1/(mu1 + mu2)。
  如果未提供速度和 Lame 参数，则使用标准地壳值：vp = 5800 m/s 和 vs = 3200 m/s。如果未提供 Q 值，则使用标准地壳值：qp = 1456 和 qs = 600。如果未提供 Burgers 材料参数，则瞬时和稳态粘度为 0，alpha=1。
  所有单位均为国际单位制（SI）（m/s、Pa、kg/m^3），除非明确说明。

主要材料属性被视为独立的，并且可以作为属性访问（可以对这些属性进行赋值）：

vp、vs、rho、qp、qs
其他材料属性被视为依赖属性，并可以通过实例方法查询。

实际应用

看了两天教程，终于敢下手了，好的，让我们拭目以待，是大成功还是翻车呢！（博士论文写巨慢，还搞这些有的没的，我真是饿了，啊啊啊啊）
好的，我想要这样的模型，一般呢我只用Vp，当然我好像看到了这里有自动的Vp和Vs之间的相互转换。vst=cake.castagna_vp_to_vs(vpt)

我比较想选择第三种方式，就是在原有模型的基础上进行魔改，好的，那么我们看看原有模型


从入门到放弃，对不起大家 就做到这里吧，再做就不礼貌了。吭哧瘪肚做两天，还是没做出来。
下面是我的残余代码

from pyrocko import cake
import numpy as np
mod1 = cake.load_model(fn='ak135-f-continental.m', format='nd', crust2_profile=None)
deep1 = np.array([[0], [4], [5], [9],[12], [15], [16], [17], [20], [25], [34], [36]])
vptp1 = np.array([[5.2], [5.6], [5.75], [5.85],[5.95], [6.05], [6.15], [6.25], [6.35], [6.55], [7.2], [7.8]])
deep1*=1000
vptp1*=1000

modcust=[]
nam=[]
s=0
for s, (depth, vpt) in enumerate(zip(deep1, vptp1)):
    #depth=deep1[0]
    #vpt=vptp1[0]
    print(depth,vpt)
    dept_cal=depth+500
    sss=mod1.layer(dept_cal, direction=4)
    vst=cake.castagna_vp_to_vs(vpt)
    try:
        rhot = sss.m.rho
    except AttributeError:
        rhot = sss.mtop.rho
    mt = cake.Material(vp=vpt,vs=vst, rho=rhot, qp=None)
    modcust.append(mt)
    nam1='layer'+str(s)
    nam.append(nam1)


producer = [(depth, ss, namee) for depth,ss,namee in zip(deep1, modcust,nam)]
mod2 = cake.LayeredModel.from_scanlines(producer)

for s, (depth, namee) in enumerate(zip(deep1, nam)):
    if s != 0:
        mod2.Interface(z=depth, mabove=nam[s-1], mbelow=namee)

#modcust.append(modcust1)
#modcust.append(cake.Discontinuity(depth))
#producer = [(depth, material, 'homogeneous layer') for depth,material in zip(deep1,modcust)]

#producer = [(depth, material, 'homogeneous layer') for depth,ss in zip(deep1,modcust)]
#mod2 = cake.LayeredModel.from_scanlines(producer)

mod1是我们引入的，mod2是我生成的，不一样，而且我不知道怎么给它们拼接上
mod2是这个样子的，我用了很多种方法，生成的东西都很奇怪。我觉得用不了，但是这些生成的数可以用来修改模型，也算是不幸中的万幸吧。【vp vs 密度 不记得 不记得】P波速度是模型中给定的，S波速度是用自带公式计算的。

自定义速度模型-TAUP（java）

坚持了很久，最后不得已还是得手动改，因为实在改不明白。

找到.nd文件

这个直接搜索就行，这个还是pyrocko安装包里的，随便找一个复制到工作目录。

打开以后是这样的

拿着上面pyrocko生成的数据直接修改。现在工作目录里就有一个名字叫ak135_change.m.nd的文件了

为了提高速度，可以生成一个一个.taup的文件，直接在java终端输入，然后你的工作目录就会生成同名的.taup文件

taup_create -nd ak135_change.m.nd -verbose

可以用这个命令计算震源深度为12.5km，震中距离为24.21km 的台站理论P波到时

taup_time -mod ak135_change.m.taup -ph p,P,Pn -h 12.5 -km 24.21

好的，那么接下来怎么批量呢，因为我们是windows下使用，所以用.bat而不是.sh文件。
这是一个简单的批量计算的文件。

@echo off

rem 定义要运行的多个 km 值
set km_values=22.4 30 40 50

rem 删除现有的输出文件
del output.txt

rem 循环遍历每个 km 值
for %%i in (%km_values%) do (
    rem 运行 taup_time 命令，并将输出追加到文件 output.txt
    taup_time -mod ak135_change.m.taup -ph p,P,Pn -h 10 -km %%i >> output.txt
)

更多的课题我准备以后再开发啦。