diff --git a/.github/workflows/build.yml b/.github/workflows/build.yml index fd2d5d39..d66eadb8 100644 --- a/.github/workflows/build.yml +++ b/.github/workflows/build.yml @@ -277,6 +277,14 @@ jobs: # don't use $(pwd) to get abspath, not valid on Windows echo "EXAMPLE_COPY_PATH=${{ env.PACK_NAME }}/test_tmp" >> $GITHUB_ENV + - name: Test compare_c_py + working-directory: ${{ env.EXAMPLE_COPY_PATH }}/compare_c_py + run: | + # 检查C程序和Python程序的计算结果是否一致 + chmod +x *.sh + ./run_grt.sh + python run_pygrt.py + - name: Test 1 compare_results working-directory: ${{ env.EXAMPLE_COPY_PATH }}/compare_results run: | diff --git a/.gitignore b/.gitignore index 9d23b79f..f63bfec5 100755 --- a/.gitignore +++ b/.gitignore @@ -10,7 +10,6 @@ /for_paper -/docs /pygrt/C_extension/build/ /pygrt/C_extension/lib/ diff --git a/.readthedocs.yaml b/.readthedocs.yaml new file mode 100644 index 00000000..bbb6f60d --- /dev/null +++ b/.readthedocs.yaml @@ -0,0 +1,43 @@ +# Read the Docs configuration file for Sphinx projects +# See https://docs.readthedocs.io/en/stable/config-file/v2.html for details + +# Required +version: 2 + +# Set the OS, Python version and other tools you might need +build: + os: ubuntu-22.04 + tools: + python: "3.10" + # You can also specify other tool versions: + # nodejs: "20" + # rust: "1.70" + # golang: "1.20" + + apt_packages: ["doxygen"] + + commands: + - cd docs && chmod +x *.sh && ./compile.sh + + +# Build documentation in the "docs/" directory with Sphinx +sphinx: + configuration: docs/source/conf.py + # You can configure Sphinx to use a different builder, for instance use the dirhtml builder for simpler URLs + # builder: "dirhtml" + # Fail on all warnings to avoid broken references + # fail_on_warning: true + +# Optionally build your docs in additional formats such as PDF and ePub +# formats: +# - pdf +# - epub + +# Optional but recommended, declare the Python requirements required +# to build your documentation +# See https://docs.readthedocs.io/en/stable/guides/reproducible-builds.html +python: + install: + - requirements: docs/requirements.txt + + diff --git a/docs/Makefile b/docs/Makefile new file mode 100755 index 00000000..8379f040 --- /dev/null +++ b/docs/Makefile @@ -0,0 +1,24 @@ +# Minimal makefile for Sphinx documentation +# + +# You can set these variables from the command line, and also +# from the environment for the first two. +SPHINXOPTS ?= +SPHINXBUILD ?= sphinx-build +SOURCEDIR = source +BUILDDIR = build + +# Put it first so that "make" without argument is like "make help". +help: + @$(SPHINXBUILD) -M help "$(SOURCEDIR)" "$(BUILDDIR)" $(SPHINXOPTS) $(O) + +clean: + rm $(BUILDDIR)/* -rf + rm doxygen_h/* -rf + +.PHONY: help Makefile clean + +# Catch-all target: route all unknown targets to Sphinx using the new +# "make mode" option. $(O) is meant as a shortcut for $(SPHINXOPTS). +%: Makefile + @$(SPHINXBUILD) -M $@ "$(SOURCEDIR)" "$(BUILDDIR)" $(SPHINXOPTS) $(O) diff --git a/docs/compile.sh b/docs/compile.sh new file mode 100755 index 00000000..886b5c3d --- /dev/null +++ b/docs/compile.sh @@ -0,0 +1,16 @@ +#!/bin/bash + +chmod +x *.sh + +# 清空 +make clean + +# 生成api对应的.rst文件 +./create_api_rst.sh + +# 执行doxygen +doxygen doxyfile_h + +# 构建 +make html +# sphinx-autobuild --port 8000 source build \ No newline at end of file diff --git a/docs/create_api_rst.sh b/docs/create_api_rst.sh new file mode 100755 index 00000000..7f1ec767 --- /dev/null +++ b/docs/create_api_rst.sh @@ -0,0 +1,98 @@ +#!/bin/bash + +# 根据源代码中include/路径下的.h文件，自动在source/API/C_extension文件夹下生成.rst文件 + +APIDIR="source/API" + +PDIR="$APIDIR/C_extension" +rm -rf $PDIR +mkdir -p $PDIR + +SDIR="../pygrt/C_extension/include" + +c_api_rst="$APIDIR/c_api.rst" + +cat > $c_api_rst < $PDIR/$dir1.rst < $PDIR/$dir1/${hnm%%.*}.rst <> $PDIR/$dir1.rst + done + + echo " C_extension/$dir1" >> $c_api_rst +done + + + +# 根据源代码中pygrt/路径下的.py文件，自动在source/API/pygrt文件夹下生成.rst文件 +PDIR="$APIDIR/pygrt" +rm -rf $PDIR +mkdir -p $PDIR + +SDIR="../pygrt" + +py_api_rst="$APIDIR/py_api.rst" + +cat > $py_api_rst < $PDIR/${py%%.*}.rst <> $py_api_rst + +done diff --git a/docs/doxyfile_h b/docs/doxyfile_h new file mode 100755 index 00000000..ce7a5bbc --- /dev/null +++ b/docs/doxyfile_h @@ -0,0 +1,32 @@ +# 头文件配置文档 + +# The INPUT tag is used to specify the files and/or directories that contain +# documented source files. You may enter file names like "myfile.cpp" or +# directories like "/usr/src/myproject". Separate the files or directories +# with spaces. + +INPUT = ../pygrt/C_extension/include + +RECURSIVE = YES + +# If the value of the INPUT tag contains directories, you can use the +# FILE_PATTERNS tag to specify one or more wildcard patterns (like *.cpp +# and *.h) to filter out the source-files in the directories. If left +# blank, Doxygen will consider all files in the directories. + +FILE_PATTERNS = *.h + +# The OUTPUT_DIRECTORY tag is used to specify the (relative or absolute) +# base path where the generated documentation will be put. + +OUTPUT_DIRECTORY = ./doxygen_h + +# If the GENERATE_XML tag is set to YES Doxygen will generate an XML file +# that captures the structure of the code including all documentation. + +GENERATE_XML = YES + +# Doxyfile Configuration +USE_MATHJAX = YES + +MARKDOWN_SUPPORT = YES \ No newline at end of file diff --git a/docs/requirements.txt b/docs/requirements.txt new file mode 100755 index 00000000..e71c4a6f --- /dev/null +++ b/docs/requirements.txt @@ -0,0 +1,15 @@ +breathe +commonmark +recommonmark +sphinx +sphinx-rtd-theme +sphinx-copybutton +sphinx_markdown_tables +sphinxcontrib-mermaid +sphinxcontrib-htmlhelp +sphinxcontrib-jsmath +nbsphinx +graphviz +sphinx-tabs +setuptools-scm # 自动从git tag中构建版本 +# git+https://github.com/Dengda98/PyGRT.git # 包本身 diff --git a/docs/source/API/api.rst b/docs/source/API/api.rst new file mode 100644 index 00000000..4467eea4 --- /dev/null +++ b/docs/source/API/api.rst @@ -0,0 +1,13 @@ +API +================== + +:Author: Zhu Dengda +:Email: zhudengda@mail.iggcas.ac.cn + +----------------------------------------------------------- + +.. toctree:: + :maxdepth: 2 + + c_api + py_api \ No newline at end of file diff --git a/docs/source/API/c_api.rst b/docs/source/API/c_api.rst new file mode 100755 index 00000000..3402dafd --- /dev/null +++ b/docs/source/API/c_api.rst @@ -0,0 +1,15 @@ +C extension API +=============================== + +:Author: Zhu Dengda +:Email: zhudengda@mail.iggcas.ac.cn + +----------------------------------------------------------- + +.. toctree:: + :maxdepth: 2 + + C_extension/common + C_extension/dynamic + C_extension/static + C_extension/travt diff --git a/docs/source/API/py_api.rst b/docs/source/API/py_api.rst new file mode 100755 index 00000000..2d8a3b47 --- /dev/null +++ b/docs/source/API/py_api.rst @@ -0,0 +1,18 @@ +**PyGRT** package API +======================== + +:Author: Zhu Dengda +:Email: zhudengda@mail.iggcas.ac.cn + +----------------------------------------------------------- + + +.. toctree:: + :maxdepth: 2 + + pygrt/c_interfaces + pygrt/c_structures + pygrt/pygrn + pygrt/pymod + pygrt/signals + pygrt/utils diff --git a/docs/source/Tutorial/dynamic/gfunc.rst b/docs/source/Tutorial/dynamic/gfunc.rst new file mode 100644 index 00000000..e684aac1 --- /dev/null +++ b/docs/source/Tutorial/dynamic/gfunc.rst @@ -0,0 +1,181 @@ +.. _gfunc_rst: + +计算动态格林函数 +================= + +:Author: Zhu Dengda +:Email: zhudengda@mail.iggcas.ac.cn + +----------------------------------------------------------- + + +Python中计算动态格林函数的主函数为 :func:`compute_grn() ` ，C程序为 :command:`grt`。 + +核心计算逻辑来自 :ref:`初稿 ` ，具体代码可见与C API中对应同名 :file:`*.c` 文件，其中计算格林函数频谱的主函数为 :file:`grt.c` 里的 :c:func:`integ_grn_spec`。 在公式推导过程中选取的柱坐标系三分量Z（垂向）、R（径向）、T（切向）的正方向分别为垂直向下、震源指向台站向外的方向、相对R方向顺时针旋转90度。 + +.. _warning_C_python_Z_direction: +.. warning:: + + **推导过程中定义的Z分量与C/Python输出的Z分量方向相反**。因为从Z轴向下为正的坐标系对于 **半无限水平分层均匀介质** 来讲更便于公式推导；而在结果输出后，会对Z分量取相反数，因为Z轴向上和目前实采数据使用的Z轴方向一致，这更便于后续数据处理。 + + 除非特别指明（例如描述震源机制），且一般使用 **PyGRT** 时也不会涉及底层程序内部， **在实际使用过程中提及的Z轴均以向上为正** ，仅在讨论代码内部细节和公式推导时对此加以区分。 + + + +示例程序 +----------- + +假设在 :file:`milrow` 模型中，震源深度2km，接收点位于地表，震中距为5km，8km和10km，计算格林函数，要求采样点数500个点，采样间隔0.02s。 + +.. tabs:: + + .. tab:: C + + .. literalinclude:: run/run.sh + :language: bash + :start-after: BEGIN GRN + :end-before: END GRN + + 不同震源深度、接收点深度和震中距的格林函数会在 :rst:dir:`GRN/milrow_{depsrc}_{deprcv}_{dist}/`` 路径下，使用SAC格式保存。 + + 一些基本信息（包括源点和场点的物性参数）保存在SAC头段变量中，其中 :c:var:`t0` 和 :c:var:`t1` 分别代表初至P波和初至S波的到时，由于程序中使用0作为参考时间，故其等价于走时。 + + .. tab:: Python + + .. literalinclude:: run/run.py + :language: python + :start-after: BEGIN GRN + :end-before: END GRN + + 多个震中距的格林函数以列表形式返回，其中每个元素为 |Stream| 类。:class:`Trace.stats.sac` 中保存了SAC头段变量，与C程序输出保持一致。 + + + + +离散波数积分 +-------------- + +格林函数频谱的计算本质转化为求以下积分: + +.. math:: + + P_m(\omega) = \int_0^\infty F(k, \omega)J_m(kr)kdk + +其中 :math:`F(k,\omega)` 称为核函数，它是和介质属性相关的量，与震中距无关。我们可以使用离散波数积分法 :ref:`(Bouchon, 1981) ` 将上式积分转变为求和: + +.. math:: + + P_m(\omega) = \Delta k \sum_{j=0}^{\infty} F(k_j,\omega)J_m(k_j r)k_j + +其中 :math:`\Delta k = 2\pi/L, k_j=j\Delta k`，:math:`L` 为特征长度，要求满足: + +.. math:: + + \left\{ + \begin{aligned} + & r < L/2 \\ + & (L-r)^2 + z_s^2 > (\alpha T)^2 + \end{aligned} + \right. + +其中 :math:`\alpha` 为参考P波速度， :math:`T` 为所需计算的理论地震图的总时间长度。常见的保守经验值为 :math:`L=20r` ，但也应依具体情况而定 。为了避开附加源以及奇点的影响，:ref:`(Bouchon, 1981) ` 在频率上添加微小虚部，具体推导过程详见 :ref:`(Bouchon, 1981) ` 和 :ref:`(张海明, 2021) `。 + + + + +积分形式分类 +-------------- + +通过在面谐矢量坐标系中建立波函数进行公式推导，最终理论地震图的三分量频谱 :math:`W_m(\omega), Q_m(\omega), V_m(\omega)` (分别为垂向，径向，切向)可以表达为: + +.. math:: + + \left\{ + \begin{aligned} + W_m(\omega) &= \int_0^\infty w_m J_m(kr)kdk \\ + Q_m(\omega) &= \int_0^\infty (q_m J_m^{\prime}(kr) - v_m \frac{m}{kr} J_m(kr)) kdk \\ + V_m(\omega) &= \int_0^\infty (q_m \frac{m}{kr} J_m(kr) - v_m J_m^{\prime}(kr)) kdk + \end{aligned} + \right. + +.. note:: + + 初次推导该公式可能会对虚数 :math:`i` 及公式中的正负号感到疑惑，但其实这里的设计是将虚数 :math:`i` 和方向因子 :math:`e^{im\theta}` 合并，所以在后续合成理论地震图时你会发现，:math:`m=0,1,2` 阶的 :math:`W_m, Q_m` 的方向因子对 :math:`\theta` 的偏导就是 :math:`V_m` 的方向因子。 + + +公式来自 :ref:`初稿 ` (5.6.22)式，其中阶数 :math:`m=0,1,2`。核函数 :math:`q_m,w_m,v_m` 根据广义反射透射系数矩阵法求得。为了方便程序实现，根据积分形式，我们对待求积分进行如下分类，其中每一阶都分为4类( :math:`p=0,1,2,3` )，除了0阶只需两类，此时 :math:`v_0=0` : + ++ :math:`m=0` + +.. math:: + + \left\{ + \begin{aligned} + p=0 & \rightarrow - \int q_0(k, \omega) J_0(kr)kdk \\ + p=2 & \rightarrow \int w_0(k, \omega) J_1(kr)kdk + \end{aligned} + \right. + + ++ :math:`m=1,2` + +.. math:: + + \left\{ + \begin{aligned} + p=0 & \rightarrow \int q_m(k, \omega) J_{m-1}(kr)kdk \\ + p=1 & \rightarrow - \int (q_m(k, \omega) + v_m(k, \omega)) \frac{m}{kr} J_m(kr)kdk \\ + p=2 & \rightarrow \int w_m(k, \omega) J_m(kr)kdk \\ + p=3 & \rightarrow - \int v_m(k, \omega) J_{m-1}(kr)kdk + \end{aligned} + \right. + + +以上每个积分都形成 :math:`\int_0^\infty F(k, \omega)J_m(kr)kdk` 的形式，便可逐个使用离散波数积分（或Filon积分、峰谷平均法等）求解每个积分。 + + +.. _grn_types: + +格林函数分类 +-------------- + +程序会输出15个格林函数（也可以选择输出哪些震源），但并不是每个震源类型对应的每一阶每种积分类型都存在。以下为15个格林函数定义的名称，以及对应上述的阶数以及积分类型: + ++----------+-------------------+--------------+----------------------+ +| **名称** | **格林函数类型** | **对应阶数** | **对应积分类型** | ++----------+-------------------+--------------+----------------------+ +| EXZ | 爆炸源Z分量 | :math:`m=0` | :math:`p=2` | ++----------+-------------------+--------------+----------------------+ +| EXR | 爆炸源R分量 | :math:`m=0` | :math:`p=0` | ++----------+-------------------+--------------+----------------------+ +| VFZ | 垂直向下力源Z分量 | :math:`m=0` | :math:`p=2` | ++----------+-------------------+--------------+----------------------+ +| VFR | 垂直向下力源R分量 | :math:`m=0` | :math:`p=0` | ++----------+-------------------+--------------+----------------------+ +| HFZ | 水平力源Z分量 | :math:`m=1` | :math:`p=2` | ++----------+-------------------+--------------+----------------------+ +| HFR | 水平力源R分量 | :math:`m=1` | :math:`(p=0)+(p=1)` | ++----------+-------------------+--------------+----------------------+ +| HFT | 水平力源T分量 | :math:`m=1` | :math:`-(p=1)+(p=3)` | ++----------+-------------------+--------------+----------------------+ +| DDZ | 倾角45度倾滑Z分量 | :math:`m=0` | :math:`p=2` | ++----------+-------------------+--------------+----------------------+ +| DDR | 倾角45度倾滑R分量 | :math:`m=0` | :math:`p=0` | ++----------+-------------------+--------------+----------------------+ +| DSZ | 倾角90度倾滑Z分量 | :math:`m=1` | :math:`p=2` | ++----------+-------------------+--------------+----------------------+ +| DSR | 倾角90度倾滑R分量 | :math:`m=1` | :math:`(p=0)+(p=1)` | ++----------+-------------------+--------------+----------------------+ +| DST | 倾角90度倾滑T分量 | :math:`m=1` | :math:`-(p=1)+(p=3)` | ++----------+-------------------+--------------+----------------------+ +| SSZ | 倾角90度走滑Z分量 | :math:`m=2` | :math:`p=2` | ++----------+-------------------+--------------+----------------------+ +| SSR | 倾角90度走滑R分量 | :math:`m=2` | :math:`(p=0)+(p=1)` | ++----------+-------------------+--------------+----------------------+ +| SST | 倾角90度走滑T分量 | :math:`m=2` | :math:`-(p=1)+(p=3)` | ++----------+-------------------+--------------+----------------------+ + + + + + diff --git a/docs/source/Tutorial/dynamic/integ_converg.rst b/docs/source/Tutorial/dynamic/integ_converg.rst new file mode 100644 index 00000000..14efb24d --- /dev/null +++ b/docs/source/Tutorial/dynamic/integ_converg.rst @@ -0,0 +1,13 @@ +积分收敛性 +=================== + +:Author: Zhu Dengda +:Email: zhudengda@mail.iggcas.ac.cn + +----------------------------------------------------------- + +通过输出积分过程文件，直观感受源点和场点深度接近时，积分收敛性的变化，以及峰谷平均法的作用。 + +动态和静态情况类似，这里以动态解为例。 + +TODO \ No newline at end of file diff --git a/docs/source/Tutorial/dynamic/run/milrow b/docs/source/Tutorial/dynamic/run/milrow new file mode 100644 index 00000000..a99f4912 --- /dev/null +++ b/docs/source/Tutorial/dynamic/run/milrow @@ -0,0 +1,9 @@ +0.2 3.4 1.7 2.3 9e10 9e10 +0.6 3.7 1.9 2.4 9e10 9e10 +0.5 4.2 2.1 2.4 9e10 9e10 +0.5 4.6 2.3 2.5 9e10 9e10 +0.7 4.9 2.8 2.6 9e10 9e10 +0.5 5.1 2.9 2.7 9e10 9e10 +6.0 5.9 3.3 2.7 9e10 9e10 +1. 6.9 4.0 2.8 9e10 9e10 +0.0 8.2 4.7 3.2 9e10 9e10 \ No newline at end of file diff --git a/docs/source/Tutorial/dynamic/run/run.py b/docs/source/Tutorial/dynamic/run/run.py new file mode 100644 index 00000000..80bfc62d --- /dev/null +++ b/docs/source/Tutorial/dynamic/run/run.py @@ -0,0 +1,168 @@ +# START BUILD MODEL +import numpy as np +import pygrt + +# option 1: +# modarr = np.loadtxt("milrow") + +# option 2 +modarr = np.array([ + [0.2, 3.4, 1.7, 2.3, 9e10, 9e10], + [0.6, 3.7, 1.9, 2.4, 9e10, 9e10], + [0.5, 4.2, 2.1, 2.4, 9e10, 9e10], + [0.5, 4.6, 2.3, 2.5, 9e10, 9e10], + [0.7, 4.9, 2.8, 2.6, 9e10, 9e10], + [0.5, 5.1, 2.9, 2.7, 9e10, 9e10], + [6.0, 5.9, 3.3, 2.7, 9e10, 9e10], + [28., 6.9, 4.0, 2.8, 9e10, 9e10], + [0., 8.2, 4.7, 3.2, 9e10, 9e10], +]) +# END BUILD MODEL + +from obspy import * +import matplotlib.pyplot as plt + +def plot_syn(stsyn:Stream, out:str|None=None, sigs:np.ndarray|None=None): + figsize = (10, 4) + nrow = 3 + if sigs is not None: + nrow += 1 + figsize = (10, 4.5) + + fig, axs = plt.subplots(nrow, 1, figsize=figsize, gridspec_kw=dict(hspace=0.0), sharex=True) + nt = stsyn[0].stats.npts + dt = stsyn[0].stats.delta + t = np.arange(nt)*dt + + if sigs is not None: + ax = axs[0] + ax.plot(t[:len(sigs)], sigs, 'k-', lw=0.5) + axs = axs[1:] + + travtP = stsyn[0].stats.sac['t0'] + travtS = stsyn[0].stats.sac['t1'] + + for i, tr in enumerate(stsyn): + ax = axs[i] + ax.plot(t, tr.data, c='k', lw=0.5, label=tr.stats.channel) + ax.legend(loc='upper left') + + ylims = ax.get_ylim() + # 绘制到时 + ax.vlines(travtP, *ylims, colors='b') + ax.text(travtP, ylims[1], "P", ha='left', va='top', color='b') + ax.vlines(travtS, *ylims, colors='r') + ax.text(travtS, ylims[1], "S", ha='left', va='top', color='r') + + ax.set_xlim([t[0], t[-1]]) + ax.set_ylim(np.array(ylims)*1.2) + + if out is not None: + fig.savefig(out, dpi=100) + + +# BEGIN GRN +modarr = np.loadtxt("milrow") + +pymod = pygrt.PyModel1D(modarr, depsrc=2.0, deprcv=0.0) + +# 多个震中距的格林函数以列表形式返回，其中每个元素为 |Stream| 类。 +stgrnLst = pymod.compute_grn( + distarr=[5,8,10], + nt=500, dt=0.02 +) + +print(stgrnLst[0]) +# 15 Trace(s) in Stream: +# .SYN..EXZ | 1970-01-01T00:00:00.000000Z - 1970-01-01T00:00:09.980000Z | 50.0 Hz, 500 samples +# .SYN..VFZ | 1970-01-01T00:00:00.000000Z - 1970-01-01T00:00:09.980000Z | 50.0 Hz, 500 samples +# ... +# END GRN + + +# BEGIN SYN EXP +# 接之前的代码 +idx = 2 +stgrn = stgrnLst[idx] # 选择格林函数 + +stsyn = pygrt.utils.gen_syn_from_gf_EXP(stgrn, M0=1e24, az=30) +print(stsyn) +# 3 Trace(s) in Stream: +# .SYN..EXZ | 1970-01-01T00:00:00.000000Z - 1970-01-01T00:00:09.980000Z | 50.0 Hz, 500 samples +# .SYN..EXR | 1970-01-01T00:00:00.000000Z - 1970-01-01T00:00:09.980000Z | 50.0 Hz, 500 samples +# .SYN..EXT | 1970-01-01T00:00:00.000000Z - 1970-01-01T00:00:09.980000Z | 50.0 Hz, 500 samples +# END SYN EXP +plot_syn(stsyn, "syn_exp.png") + +# BEGIN SYN SF +# 接之前的代码 +idx = 2 +stgrn = stgrnLst[idx] # 选择格林函数 + +stsyn = pygrt.utils.gen_syn_from_gf_SF(stgrn, S=1e16, fN=1, fE=-0.5, fZ=2, az=30) +print(stsyn) +# 3 Trace(s) in Stream: +# .SYN..SFZ | 1970-01-01T00:00:00.000000Z - 1970-01-01T00:00:09.980000Z | 50.0 Hz, 500 samples +# .SYN..SFR | 1970-01-01T00:00:00.000000Z - 1970-01-01T00:00:09.980000Z | 50.0 Hz, 500 samples +# .SYN..SFT | 1970-01-01T00:00:00.000000Z - 1970-01-01T00:00:09.980000Z | 50.0 Hz, 500 samples +# END SYN SF +plot_syn(stsyn, "syn_sf.png") + + +# BEGIN SYN DC +# 接之前的代码 +idx = 2 +stgrn = stgrnLst[idx] # 选择格林函数 + +stsyn = pygrt.utils.gen_syn_from_gf_DC(stgrn, M0=1e24, strike=33, dip=50, rake=120, az=30) +print(stsyn) +# 3 Trace(s) in Stream: +# .SYN..DCZ | 1970-01-01T00:00:00.000000Z - 1970-01-01T00:00:09.980000Z | 50.0 Hz, 500 samples +# .SYN..DCR | 1970-01-01T00:00:00.000000Z - 1970-01-01T00:00:09.980000Z | 50.0 Hz, 500 samples +# .SYN..DCT | 1970-01-01T00:00:00.000000Z - 1970-01-01T00:00:09.980000Z | 50.0 Hz, 500 samples +# END SYN DC +plot_syn(stsyn, "syn_dc.png") + + +# BEGIN SYN MT +idx = 2 +stgrn = stgrnLst[idx] # 选择格林函数 + +stsyn = pygrt.utils.gen_syn_from_gf_MT(stgrn, M0=1e24, MT=[0.1,-0.2,1.0,0.3,-0.5,-2.0], az=30) +print(stsyn) +# 3 Trace(s) in Stream: +# .SYN..MTZ | 1970-01-01T00:00:00.000000Z - 1970-01-01T00:00:09.980000Z | 50.0 Hz, 500 samples +# .SYN..MTR | 1970-01-01T00:00:00.000000Z - 1970-01-01T00:00:09.980000Z | 50.0 Hz, 500 samples +# .SYN..MTT | 1970-01-01T00:00:00.000000Z - 1970-01-01T00:00:09.980000Z | 50.0 Hz, 500 samples +# END SYN MT +plot_syn(stsyn, "syn_mt.png") + + +# BEGIN ZNE +# 接之前的代码 +idx = 2 +stgrn = stgrnLst[idx] # 选择格林函数 + +# 设置ZNE=True可返回ZNE分量 +stsyn = pygrt.utils.gen_syn_from_gf_DC(stgrn, M0=1e24, strike=33, dip=50, rake=120, az=30, ZNE=True) +print(stsyn) +# 3 Trace(s) in Stream: +# .SYN..DCZ | 1970-01-01T00:00:00.000000Z - 1970-01-01T00:00:09.980000Z | 50.0 Hz, 500 samples +# .SYN..DCN | 1970-01-01T00:00:00.000000Z - 1970-01-01T00:00:09.980000Z | 50.0 Hz, 500 samples +# .SYN..DCE | 1970-01-01T00:00:00.000000Z - 1970-01-01T00:00:09.980000Z | 50.0 Hz, 500 samples +plot_syn(stsyn, "syn_dc_zne.png") +# END ZNE + + + +# BEGIN TIME FUNC +idx = 2 +stgrn = stgrnLst[idx] # 选择格林函数 + +stsyn = pygrt.utils.gen_syn_from_gf_SF(stgrn, S=1e16, fN=1, fE=-0.5, fZ=2, az=30) +# 生成时间函数 +trig = pygrt.sigs.gen_triangle_wave(0.6, 0.02) +# 卷积，原地修改 +pygrt.utils.stream_convolve(stsyn, trig) +# END TIME FUNC +plot_syn(stsyn, "syn_sf_trig.png", trig) \ No newline at end of file diff --git a/docs/source/Tutorial/dynamic/run/run.sh b/docs/source/Tutorial/dynamic/run/run.sh new file mode 100755 index 00000000..af4e399f --- /dev/null +++ b/docs/source/Tutorial/dynamic/run/run.sh @@ -0,0 +1,41 @@ +#!/bin/bash + +# BEGIN GRN +# 不同震源深度、接收点深度和震中距的格林函数会在 +# GRN/milrow_{depsrc}_{deprcv}_{dist}/ 路径下，使用SAC格式保存。 +grt -Mmilrow -D2/0 -N500/0.02 -OGRN -R5,8,10 +# END GRN + +# BEGIN SYN EXP +# 合成结果在 syn_exp/ 目录下，以SAC格式保存。 +grt.syn -GGRN/milrow_2_0_10 -S1e24 -A30 -Osyn_exp +# END SYN EXP + +# BEGIN SYN SF +# 合成结果在 syn_sf/ 目录下，以SAC格式保存。 +grt.syn -GGRN/milrow_2_0_10 -S1e16 -A30 -F1/-0.5/2 -Osyn_sf +# END SYN SF + +# BEGIN SYN DC +# 合成结果在 syn_dc/ 目录下，以SAC格式保存。 +grt.syn -GGRN/milrow_2_0_10 -S1e24 -A30 -M33/50/120 -Osyn_dc +# END SYN DC + +# BEGIN SYN MT +# 合成结果在 syn_mt/ 目录下，以SAC格式保存。 +grt.syn -GGRN/milrow_2_0_10 -S1e24 -A30 -T0.1/-0.2/1.0/0.3/-0.5/-2.0 -Osyn_mt +# END SYN MT + + +# BEGIN ZNE +# 传入-N可输出ZNE分量 +grt.syn -GGRN/milrow_2_0_10 -S1e24 -A30 -M33/50/120 -Osyn_dc_zne -N +# END ZNE + + +# BEGIN TIME FUNC +# -Dt定义梯形波的三个时间参数，爬升截止时刻/平稳截止时刻/结束时刻 +# 前两个时刻重合则退化为三角波 +grt.syn -GGRN/milrow_2_0_10 -S1e16 -A30 -F1/-0.5/2 -Osyn_sf_trig -Dt/0.3/0.3/0.6 +# END TIME FUNC + diff --git a/docs/source/Tutorial/dynamic/run/syn_dc.png b/docs/source/Tutorial/dynamic/run/syn_dc.png new file mode 100644 index 00000000..35001f24 Binary files /dev/null and b/docs/source/Tutorial/dynamic/run/syn_dc.png differ diff --git a/docs/source/Tutorial/dynamic/run/syn_dc_zne.png b/docs/source/Tutorial/dynamic/run/syn_dc_zne.png new file mode 100644 index 00000000..ab339990 Binary files /dev/null and b/docs/source/Tutorial/dynamic/run/syn_dc_zne.png differ diff --git a/docs/source/Tutorial/dynamic/run/syn_exp.png b/docs/source/Tutorial/dynamic/run/syn_exp.png new file mode 100644 index 00000000..bf447659 Binary files /dev/null and b/docs/source/Tutorial/dynamic/run/syn_exp.png differ diff --git a/docs/source/Tutorial/dynamic/run/syn_mt.png b/docs/source/Tutorial/dynamic/run/syn_mt.png new file mode 100644 index 00000000..44deef09 Binary files /dev/null and b/docs/source/Tutorial/dynamic/run/syn_mt.png differ diff --git a/docs/source/Tutorial/dynamic/run/syn_sf.png b/docs/source/Tutorial/dynamic/run/syn_sf.png new file mode 100644 index 00000000..ef1bb081 Binary files /dev/null and b/docs/source/Tutorial/dynamic/run/syn_sf.png differ diff --git a/docs/source/Tutorial/dynamic/run/syn_sf_trig.png b/docs/source/Tutorial/dynamic/run/syn_sf_trig.png new file mode 100644 index 00000000..599a3fad Binary files /dev/null and b/docs/source/Tutorial/dynamic/run/syn_sf_trig.png differ diff --git a/docs/source/Tutorial/dynamic/strain_stress.rst b/docs/source/Tutorial/dynamic/strain_stress.rst new file mode 100644 index 00000000..fc6815ec --- /dev/null +++ b/docs/source/Tutorial/dynamic/strain_stress.rst @@ -0,0 +1,9 @@ +计算动态应变和应力 +=================== + +:Author: Zhu Dengda +:Email: zhudengda@mail.iggcas.ac.cn + +----------------------------------------------------------- + +TODO diff --git a/docs/source/Tutorial/dynamic/syn.rst b/docs/source/Tutorial/dynamic/syn.rst new file mode 100644 index 00000000..ec169eeb --- /dev/null +++ b/docs/source/Tutorial/dynamic/syn.rst @@ -0,0 +1,168 @@ +合成动态位移 +================= + +:Author: Zhu Dengda +:Email: zhudengda@mail.iggcas.ac.cn + +----------------------------------------------------------- + + +使用上节计算的格林函数，合成动态位移（理论地震图）。方便起见，这里统一使用milrow模型，震源深度2km，场点位于地表，震中距10km的格林函数，方位角30°。 + +不同震源 +-------------- + +爆炸源 +~~~~~~~~~~~~~~~~~ +标量矩 1e24 dyne·cm。 + +.. tabs:: + + .. tab:: C + + .. literalinclude:: run/run.sh + :language: bash + :start-after: BEGIN SYN EXP + :end-before: END SYN EXP + + .. tab:: Python + + .. literalinclude:: run/run.py + :language: python + :start-after: BEGIN SYN EXP + :end-before: END SYN EXP + + +.. image:: run/syn_exp.png + :align: center + + + +单力源 +~~~~~~~~~~~~~~~~~ +北向力 :math:`f_N=1`，东向力 :math:`f_E=-0.5`，垂直向下的力 :math:`f_Z=2`，单位 1e16 dyne。 + +.. tabs:: + + .. tab:: C + + .. literalinclude:: run/run.sh + :language: bash + :start-after: BEGIN SYN SF + :end-before: END SYN SF + + .. tab:: Python + + .. literalinclude:: run/run.py + :language: python + :start-after: BEGIN SYN SF + :end-before: END SYN SF + + +.. image:: run/syn_sf.png + :align: center + + +剪切源 +~~~~~~~~~~~~~~ +断层走向33°，倾角50°，滑动角120°，标量矩 1e24 dyne·cm。 + +.. tabs:: + + .. tab:: C + + .. literalinclude:: run/run.sh + :language: bash + :start-after: BEGIN SYN DC + :end-before: END SYN DC + + .. tab:: Python + + .. literalinclude:: run/run.py + :language: python + :start-after: BEGIN SYN DC + :end-before: END SYN DC + + +.. image:: run/syn_dc.png + :align: center + + +矩张量源 +~~~~~~~~~~~~~~ +:math:`M_{xx}=0.1, M_{xy}=-0.2, M_{xz}=1.0, M_{yy}=0.3, M_{yz}=-0.5, M_{zz}=-2.0`，单位 1e24 dyne·cm， **其中X为北向，Y为东向，Z为垂直向下**。 + +.. tabs:: + + .. tab:: C + + .. literalinclude:: run/run.sh + :language: bash + :start-after: BEGIN SYN MT + :end-before: END SYN MT + + .. tab:: Python + + .. literalinclude:: run/run.py + :language: python + :start-after: BEGIN SYN MT + :end-before: END SYN MT + + +.. image:: run/syn_mt.png + :align: center + + + +分量旋转 +--------------------- +**PyGRT** 计算默认输出为ZRT分量（柱坐标系），可以设置参数以输出ZNE分量，这里以剪切源为例， + +.. tabs:: + + .. tab:: C + + .. literalinclude:: run/run.sh + :language: bash + :start-after: BEGIN ZNE + :end-before: END ZNE + + .. tab:: Python + + .. literalinclude:: run/run.py + :language: python + :start-after: BEGIN ZNE + :end-before: END ZNE + + +.. image:: run/syn_dc_zne.png + :align: center + + + +时间函数 +--------------------- +**PyGRT** 内置了一些震源时间函数，例如抛物波、梯形波、雷克子波或自定义，这里以单力源为例。 + +.. tabs:: + + .. tab:: C + + .. literalinclude:: run/run.sh + :language: bash + :start-after: BEGIN TIME FUNC + :end-before: END TIME FUNC + + 生成的时间函数会以SAC格式保存在对应路径中，文件名为 :file:`sig.sac`。 其它时间函数以及具体参数用法可使用 :command:`grt.syn -h` 查看说明。 + + .. tab:: Python + + .. literalinclude:: run/run.py + :language: python + :start-after: BEGIN TIME FUNC + :end-before: END TIME FUNC + + 其它时间函数以及具体参数用法可在 :py:mod:`pygrt.signals` 模块中查看函数参数。 + +.. image:: run/syn_sf_trig.png + :align: center \ No newline at end of file diff --git a/docs/source/Tutorial/index.rst b/docs/source/Tutorial/index.rst new file mode 100644 index 00000000..038e09a7 --- /dev/null +++ b/docs/source/Tutorial/index.rst @@ -0,0 +1,89 @@ +简易教程 +============= + +:Author: Zhu Dengda +:Email: zhudengda@mail.iggcas.ac.cn + +----------------------------------------------------------- + +**PyGRT** 程序包由C和Python两个编程语言的代码组成，目的是兼并高效性和便捷性。底层复杂运算由C语言编写，编译链接成动态库 ``libgrt.so`` 供Python调用。Python通过 ``ctypes`` 库导入动态库以使用外部函数，以此兼并了C语言的高效和Python语言的便捷。 + +除了Python脚本式运行， **PyGRT** 保留传统命令行式运行C程序，即在建立动态库过程中也编译了一些可执行文件。 **C程序的运行独立于Python，不需要Python环境，以满足更多计算场景。** 每个C程序可使用 ``-h`` 查看帮助。 + +计算 **动态解** 的主要计算流程如下： + ++ Python + +.. mermaid:: + :zoom: + + flowchart LR + + GG(["compute_grn()"]) + SS(["gen_syn_from_gf_*()"]) + EE(["compute_strain()"]) + TT(["compute_stress()"]) + + G["计算格林函数 + (及其空间导数)"] + S["根据震源机制计算位移 + (及其空间导数)"] + E["计算应变"] + T["计算应力"] + + GG --> G + G --> SS --> S + S --> EE --> E + S --> TT --> T + + classDef cmdcls fill:#f9f2d9,stroke:#e8d174,stroke-width:2px,color:#333; + class GG,SS,EE,TT cmdcls + ++ C + +.. mermaid:: + :zoom: + + flowchart LR + + GG(["grt"]) + SS(["grt.syn"]) + EE(["grt.strain"]) + TT(["grt.stress"]) + + G["计算格林函数 + (及其空间导数)"] + S["根据震源机制计算位移 + (及其空间导数)"] + E["计算应变"] + T["计算应力"] + + GG --> G + G --> SS --> S + S --> EE --> E + S --> TT --> T + + classDef cmdcls fill:#f9f2d9,stroke:#e8d174,stroke-width:2px,color:#333; + class GG,SS,EE,TT cmdcls + +计算 **静态解** 的过程相同，只是函数名称和程序名称有变化。 + + +以下是一些简易教程，可快速上手。 + + +.. toctree:: + :maxdepth: 1 + + prepare + dynamic/gfunc + dynamic/syn + static/static_gfunc + static/static_syn + dynamic/strain_stress + static/static_strain_stress + dynamic/integ_converg + + + + diff --git a/docs/source/Tutorial/prepare.rst b/docs/source/Tutorial/prepare.rst new file mode 100644 index 00000000..d5d6cb2c --- /dev/null +++ b/docs/source/Tutorial/prepare.rst @@ -0,0 +1,50 @@ +准备工作 +============= + +:Author: Zhu Dengda +:Email: zhudengda@mail.iggcas.ac.cn + +----------------------------------------------------------- + + +安装 **PyGRT** +-------------------- + +详见 :ref:`install_section` 。 + + +建立模型文件 +-------------------- + +这里的选项卡 **C** 不代表内部是C语言代码，而是 **使用C语言编译出的可执行程序** （与之相关）。后续的选项卡沿用此设定，不再解释。 + +.. tabs:: + + .. tab:: C + + **PyGRT** 以如下自由格式定义模型中每层的物性参数，每列之间以空格隔开 + + .. code-block:: text + + [厚度(km)] [P波速度(km/s)] [S波速度(km/s)] [密度(g/cm^3)] [P波品质因子] [S波品质因子] + + + 例如 :file:`milrow` 模型（假设文本文件名为 `milrow` ） + + .. literalinclude:: dynamic/run/milrow + :language: text + + + .. tab:: Python + + 模型格式与C一致，在Python中可以使用 :code:`np.loadtxt()` 导入文本文件，或者手动定义数组 + + .. literalinclude:: dynamic/run/run.py + :language: python + :start-after: START BUILD MODEL + :end-before: END BUILD MODEL + + +.. note:: + + 最后一行表示半空间，对应厚度值不会被使用。 \ No newline at end of file diff --git a/docs/source/Tutorial/static/run/milrow b/docs/source/Tutorial/static/run/milrow new file mode 100644 index 00000000..a99f4912 --- /dev/null +++ b/docs/source/Tutorial/static/run/milrow @@ -0,0 +1,9 @@ +0.2 3.4 1.7 2.3 9e10 9e10 +0.6 3.7 1.9 2.4 9e10 9e10 +0.5 4.2 2.1 2.4 9e10 9e10 +0.5 4.6 2.3 2.5 9e10 9e10 +0.7 4.9 2.8 2.6 9e10 9e10 +0.5 5.1 2.9 2.7 9e10 9e10 +6.0 5.9 3.3 2.7 9e10 9e10 +1. 6.9 4.0 2.8 9e10 9e10 +0.0 8.2 4.7 3.2 9e10 9e10 \ No newline at end of file diff --git a/docs/source/Tutorial/static/run/run.py b/docs/source/Tutorial/static/run/run.py new file mode 100644 index 00000000..4d769803 --- /dev/null +++ b/docs/source/Tutorial/static/run/run.py @@ -0,0 +1,39 @@ +# BEGIN GRN +import numpy as np +import pygrt + +modarr = np.loadtxt("milrow") + +pymod = pygrt.PyModel1D(modarr, depsrc=2.0, deprcv=0.0) + +xarr = np.linspace(-3, 3, 41) +yarr = np.linspace(-2.5, 2.5, 33) +static_grn = pymod.compute_static_grn(xarr, yarr) +print(static_grn.keys()) +# dict_keys(['_xarr', '_yarr', '_src_va', '_src_vb', '_src_rho', '_rcv_va', '_rcv_vb', '_rcv_rho', 'EXZ', 'VFZ', 'DDZ', 'HFZ', 'DSZ', 'SSZ', 'EXR', 'VFR', 'DDR', 'HFR', 'DSR', 'SSR', 'HFT', 'DST', 'SST']) +# END GRN + +# BEGIN SYN EXP +static_syn = pygrt.utils.gen_syn_from_gf_EXP(static_grn, M0=1e24, ZNE=True) +print(static_syn.keys()) +# dict_keys(['_xarr', '_yarr', '_src_va', '_src_vb', '_src_rho', '_rcv_va', '_rcv_vb', '_rcv_rho', 'EXZ', 'EXN', 'EXE']) +# END SYN EXP + + +# BEGIN SYN SF +static_syn = pygrt.utils.gen_syn_from_gf_SF(static_grn, S=1e16, fN=1, fE=-0.5, fZ=2, ZNE=True) +print(static_syn.keys()) +# dict_keys(['_xarr', '_yarr', '_src_va', '_src_vb', '_src_rho', '_rcv_va', '_rcv_vb', '_rcv_rho', 'SFZ', 'SFN', 'SFE']) +# END SYN SF + +# BEGIN SYN DC +static_syn = pygrt.utils.gen_syn_from_gf_DC(static_grn, M0=1e24, strike=33, dip=50, rake=120, ZNE=True) +print(static_syn.keys()) +# dict_keys(['_xarr', '_yarr', '_src_va', '_src_vb', '_src_rho', '_rcv_va', '_rcv_vb', '_rcv_rho', 'DCZ', 'DCN', 'DCE']) +# END SYN DC + +# BEGIN SYN MT +static_syn = pygrt.utils.gen_syn_from_gf_MT(static_grn, M0=1e24, MT=[0.1,-0.2,1.0,0.3,-0.5,-2.0], ZNE=True) +print(static_syn.keys()) +# dict_keys(['_xarr', '_yarr', '_src_va', '_src_vb', '_src_rho', '_rcv_va', '_rcv_vb', '_rcv_rho', 'MTZ', 'MTN', 'MTE']) +# END SYN MT diff --git a/docs/source/Tutorial/static/run/run.sh b/docs/source/Tutorial/static/run/run.sh new file mode 100755 index 00000000..6362e350 --- /dev/null +++ b/docs/source/Tutorial/static/run/run.sh @@ -0,0 +1,88 @@ +#!/bin/bash + + +function gmtplot_static(){ + local syn=$1 + local S=$2 + + gmt xyz2grd $syn -GsynZ.nc -R$y1/$y2/$x1/$x2 -I$ny+n/$nx+n -i0,1,2 -: + gmt xyz2grd $syn -GsynN.nc -R$y1/$y2/$x1/$x2 -I$ny+n/$nx+n -i0,1,3 -: + gmt xyz2grd $syn -GsynE.nc -R$y1/$y2/$x1/$x2 -I$ny+n/$nx+n -i0,1,4 -: + + gmt basemap -Baf -BWSen -JX5c/5c -R$y1/$y2/$x1/$x2 + gmt grdimage synZ.nc + gmt grdvector synE.nc synN.nc -Q0.08c+e+jc+h1+gblack $S + + rm *.nc +} + + + +# BEGIN GRN +depsrc=2 +deprcv=0 + +x1=-3 +x2=3 +nx=41 + +y1=-2.5 +y2=2.5 +ny=33 +# 输出到标准输出，重定向到grn文件中 +stgrt -Mmilrow -D${depsrc}/${deprcv} -X$x1/$x2/$nx -Y$y1/$y2/$ny > grn +# END GRN + +head -n6 grn > grn_head +echo "..." >> grn_head + + +# BEGIN SYN EXP +# 从标准输入中读取格林函数 +stgrt.syn -S1e24 -N < grn > syn_exp +# END SYN EXP + +gmt begin syn_exp png E300 + gmtplot_static syn_exp -Si0.03c + gmt colorbar -Bx+l"Z (cm)" +gmt end + + +# BEGIN SYN SF +# 从标准输入中读取格林函数 +stgrt.syn -S1e16 -F1/-0.5/2 -N < grn > syn_sf +# END SYN SF + +gmt begin syn_sf png E300 + gmtplot_static syn_sf -Si6.5c + gmt colorbar -Bx+l"Z (cm)" +gmt end + +# BEGIN SYN DC +# 从标准输入中读取格林函数 +stgrt.syn -S1e24 -M33/50/120 -N < grn > syn_dc +# END SYN DC + +gmt begin syn_dc png E300 + gmtplot_static syn_dc -Si0.03c + gmt meca -Sa0.5c < syn_mt +# END SYN MT + +# X Y depth mrr mtt mff mrt mrf mtf exp +# mzz mxx myy mxz -myz -mxy +gmt begin syn_mt png E300 + gmtplot_static syn_mt -Si0.02c + gmt meca -Sm0.5c <` ，C程序为 :command:`stgrt`。 + +建议先阅读完动态格林函数的计算过程（ :ref:`gfunc_rst` ）。静态情况与动态情况采取的计算方法一致，只是推导细节会有不同，详见 :ref:`初稿 ` 。 + + +示例程序 +----------- + +假设在 :file:`milrow` 模型中，震源深度2km，接收点位于地表。北向(X/km)在[-3,3]范围内等距采样41个点，东向(Y/km)在[-2.5,2.5]范围内等距采样33个点，计算这些点上的静态格林函数。 + +.. tabs:: + + .. tab:: C + + .. literalinclude:: run/run.sh + :language: bash + :start-after: BEGIN GRN + :end-before: END GRN + + 输出的 :file:`grn` 文件内部如下，开头的 ``#`` 用于保存源点和场点所在介质层的物性参数。 + + .. literalinclude:: run/grn_head + :language: text + + + .. tab:: Python + + .. literalinclude:: run/run.py + :language: python + :start-after: BEGIN GRN + :end-before: END GRN + + 函数返回字典类型，包括一些基本参数以及格林函数（2D矩阵）。 \ No newline at end of file diff --git a/docs/source/Tutorial/static/static_strain_stress.rst b/docs/source/Tutorial/static/static_strain_stress.rst new file mode 100644 index 00000000..42e34eac --- /dev/null +++ b/docs/source/Tutorial/static/static_strain_stress.rst @@ -0,0 +1,9 @@ +计算静态应变和应力 +=================== + +:Author: Zhu Dengda +:Email: zhudengda@mail.iggcas.ac.cn + +----------------------------------------------------------- + +TODO diff --git a/docs/source/Tutorial/static/static_syn.rst b/docs/source/Tutorial/static/static_syn.rst new file mode 100644 index 00000000..a2c85325 --- /dev/null +++ b/docs/source/Tutorial/static/static_syn.rst @@ -0,0 +1,116 @@ +合成静态位移 +================= + +:Author: Zhu Dengda +:Email: zhudengda@mail.iggcas.ac.cn + +----------------------------------------------------------- + +使用上节计算的格林函数，合成静态位移。为方便画图，以下结果都使用ZNE分量。 + +不同震源 +------------- + +爆炸源 +~~~~~~~~~~~~~~~~~ +标量矩 1e24 dyne·cm。 + +.. tabs:: + + .. tab:: C + + .. literalinclude:: run/run.sh + :language: bash + :start-after: BEGIN SYN EXP + :end-before: END SYN EXP + + .. tab:: Python + + .. literalinclude:: run/run.py + :language: python + :start-after: BEGIN SYN EXP + :end-before: END SYN EXP + + +.. image:: run/syn_exp.png + :width: 300px + :align: center + + +单力源 +~~~~~~~~~~~~~~~~~ +北向力 :math:`f_N=1`，东向力 :math:`f_E=-0.5`，垂直向下的力 :math:`f_Z=2`，单位 1e16 dyne。 + +.. tabs:: + + .. tab:: C + + .. literalinclude:: run/run.sh + :language: bash + :start-after: BEGIN SYN SF + :end-before: END SYN SF + + .. tab:: Python + + .. literalinclude:: run/run.py + :language: python + :start-after: BEGIN SYN SF + :end-before: END SYN SF + + +.. image:: run/syn_sf.png + :width: 300px + :align: center + + +剪切源 +~~~~~~~~~~~~~~ +断层走向33°，倾角50°，滑动角120°，标量矩 1e24 dyne·cm。 + +.. tabs:: + + .. tab:: C + + .. literalinclude:: run/run.sh + :language: bash + :start-after: BEGIN SYN DC + :end-before: END SYN DC + + .. tab:: Python + + .. literalinclude:: run/run.py + :language: python + :start-after: BEGIN SYN DC + :end-before: END SYN DC + + +.. image:: run/syn_dc.png + :width: 300px + :align: center + + +矩张量源 +~~~~~~~~~~~~~~ +:math:`M_{xx}=0.1, M_{xy}=-0.2, M_{xz}=1.0, M_{yy}=0.3, M_{yz}=-0.5, M_{zz}=-2.0`，单位 1e24 dyne·cm， **其中X为北向，Y为东向，Z为垂直向下**。 + + +.. tabs:: + + .. tab:: C + + .. literalinclude:: run/run.sh + :language: bash + :start-after: BEGIN SYN MT + :end-before: END SYN MT + + .. tab:: Python + + .. literalinclude:: run/run.py + :language: python + :start-after: BEGIN SYN MT + :end-before: END SYN MT + + +.. image:: run/syn_mt.png + :width: 300px + :align: center diff --git a/docs/source/_static/my_custom.js b/docs/source/_static/my_custom.js new file mode 100755 index 00000000..a4bc60c4 --- /dev/null +++ b/docs/source/_static/my_custom.js @@ -0,0 +1,6 @@ +// $(document).ready(function(){ +// let div_logo = document.getElementsByClassName("wy-side-nav-search")[0]; +// let a_logo = div_logo.getElementsByTagName("a"); +// a_logo[0].href = "https://github.com/Dengda98/PyGRT"; +// a_logo[0].target = "_blank"; +// }); diff --git a/docs/source/_static/my_theme.css b/docs/source/_static/my_theme.css new file mode 100755 index 00000000..73fac79a --- /dev/null +++ b/docs/source/_static/my_theme.css @@ -0,0 +1,30 @@ +/* .wy-nav-content { + max-width: 800px; +} */ + +/* 设置代码块的淡黄色背景 */ +.highlight { + background: #edfecb; /* 经典的淡黄色 */ + border-radius: 3px; /* 可选：添加圆角 */ + padding: 0.5em; /* 可选：添加内边距 */ +} + +/* 更改左上角标题背景色 */ +.wy-side-nav-search { + background-color: #F8EF9E !important; +} + +/* 如果需要，同时修改顶部标题栏颜色 */ +/* .wy-nav-top { + background-color: #FF5722 !important; +} */ + +/* 改变左上角房子图标的颜色 */ +.wy-side-nav-search a.icon-home { + color: #000000 !important; +} + +/* 改变左上角项目名的颜色 */ +.wy-side-nav-search .project { + color: #000000 !important; +} \ No newline at end of file diff --git a/docs/source/conf.py b/docs/source/conf.py new file mode 100755 index 00000000..c68c2bd3 --- /dev/null +++ b/docs/source/conf.py @@ -0,0 +1,126 @@ +# Configuration file for the Sphinx documentation builder. +# +# For the full list of built-in configuration values, see the documentation: +# https://www.sphinx-doc.org/en/master/usage/configuration.html + +# -- Project information ----------------------------------------------------- +# https://www.sphinx-doc.org/en/master/usage/configuration.html#project-information + +import sys, os, pathlib +import subprocess +import setuptools_scm + +html_last_updated_fmt = '%b %d, %Y' + + +def setup(app): + app.add_css_file('my_theme.css') + + + + + +project = 'PyGRT' +copyright = '2025, Zhu Dengda' +author = 'Zhu Dengda' +# 引入版本信息 +version = setuptools_scm.get_version(root='../..', relative_to=__file__) +release = version + +# -- General configuration --------------------------------------------------- +# https://www.sphinx-doc.org/en/master/usage/configuration.html#general-configuration + +extensions = [ + "sphinx.ext.napoleon", + "sphinx.ext.autodoc", + 'sphinx.ext.viewcode', + "recommonmark", + "sphinx_markdown_tables", + "sphinx_copybutton", + "sphinx.ext.intersphinx", + "breathe", + 'sphinx.ext.mathjax', + "nbsphinx", + 'sphinx_tabs.tabs', + 'sphinxcontrib.mermaid', + 'sphinx.ext.imgconverter', +] + + +nbsphinx_allow_errors = True # 在构建文档时允许 Notebook 中的错误 +nbsphinx_execute = 'never' # 不执行 Notebook，只是展示内容 +nbsphinx_requirejs_path = '' # 加上该选项，否则nbsphinx插件会导致mermaid无法使用zoom + +# 使用mermaid绘制简易流程图 +# mermaid_output_format = 'raw' +# mermaid_version = 'latest' +# mermaid_d3_zoom = True + +source_suffix = { + '.rst': 'restructuredtext', + # '.txt': 'markdown', + # '.md': 'markdown', +} + +myst_enable_extensions = [ + "tasklist", + "deflist", + "dollarmath", +] + +# Breathe configuration +breathe_projects = { + "h_PyGRT": "../doxygen_h/xml", +} +breathe_default_project = "h_PyGRT" + +templates_path = ['_templates'] +exclude_patterns = [] + +language = 'zh_CN' + +# -- Options for HTML output ------------------------------------------------- +# https://www.sphinx-doc.org/en/master/usage/configuration.html#options-for-html-output + +# html_theme = 'alabaster' +master_doc = 'index' +html_theme = 'sphinx_rtd_theme' +html_theme_options = { + 'analytics_anonymize_ip': False, + 'logo_only': False, # True + 'display_version': True, + 'prev_next_buttons_location': 'bottom', + 'style_external_links': False, + 'collapse_navigation': True, + 'sticky_navigation': True, + 'navigation_depth': 4, + 'includehidden': True, + 'titles_only': False, +} + +html_logo = "../../figs/logo_transparent.png" +html_static_path = ['_static'] +html_css_files = ['my_theme.css'] +html_js_files = [ + 'my_custom.js', +] + + +autodoc_default_options = { + 'member-order': 'bysource', + 'special-members': '__init__', +} + + +intersphinx_mapping = { + 'obspy': ('https://docs.obspy.org/', None), + 'numpy': ('https://numpy.org/doc/stable/', None) +} + +# 别名 +rst_epilog = f""" +.. |Stream| replace:: :class:`~obspy.core.stream.Stream` +.. |Trace| replace:: :class:`~obspy.core.trace.Trace` +""" + + \ No newline at end of file diff --git a/docs/source/copyright.rst b/docs/source/copyright.rst new file mode 100644 index 00000000..136ddac8 --- /dev/null +++ b/docs/source/copyright.rst @@ -0,0 +1,19 @@ +版权声明 +============== + + +开源许可 +-------------- + + 本项目采用 GPL-3.0 开源许可证发布，您可自由访问、使用及修改代码，但须遵守以下约束： + + + **禁止商用** ：未经书面许可，不得将本项目代码、衍生物或服务用于任何商业用途。 + + **禁止二次分发** ：您无权以任何形式（包括但不限于代码托管平台、打包发布等）公开分发本项目源代码。唯一官方下载渠道为： `PyGRT `_ + + +修改与贡献 +-------------- + + + **开源义务** ：任何基于本项目的修改或衍生作品，必须以相同许可证完全开源。 + + **贡献建议** ：欢迎通过 GitHub Pull Request 提交改进，所有贡献需遵循本项目代码规范与许可条款。 + diff --git a/docs/source/index.rst b/docs/source/index.rst new file mode 100755 index 00000000..9c8c9af2 --- /dev/null +++ b/docs/source/index.rst @@ -0,0 +1,99 @@ +**PyGRT** 文档 +===================== + +:Author: Zhu Dengda +:Email: zhudengda@mail.iggcas.ac.cn + +----------------------------------------------------------- + +**PyGRT** Github首页： `PyGRT `_ + +**PyGRT** : 一个用于计算在半无限水平分层均匀介质模型中的理论地震图的 **Python**\程序包，使用广义反射透射系数矩阵法(**G**\eneralized **R**\eflection-**T**\ransmission coefficient matrix) 以及离散波数法等方法计算。 + +**PyGRT** 能算什么？ +-------------------------- +目前，**PyGRT** 可计算区域尺度内半无限水平分层均匀介质中的 **动态与静态的位移、应力、应变。** + + +运行平台 +--------- +**Windows, Linux, MacOS** + + +程序特点 +--------- + ++ **格林函数频谱由C程序并行计算**。由于不同频率间的计算相互独立，在计算过程中基于 `OpenMP `_ 进行了 **并行优化**，极大提升了计算效率； +\ + ++ **将Python语言的便携、可扩展性与C语言的计算高效特点结合**。C程序被编译链接成动态库 *libgrt.so* ， **PyGRT** 再基于Python的 `ctypes `_ 标准库实现对C库函数的调用。再基于第三方库 `NumPy `_ 、 `SciPy `_ 和 `ObsPy `_ 可很方便地完成对C程序结果的数据整合、Fourier变换、卷积、滤波、保存到sac文件等操作（例如FFT点数不再强制要求2次幂依然能保证高效 `scipy.fft.fft `_）。借用Python语言的特点以及丰富成熟的第三方库，可灵活地实现后续的各种数据处理。 +\ + ++ **明确的中文注释**。我在代码中对计算过程给出了明确的中文注释以及相关的公式索引， **所有公式索引（除非特别指明）均来自** :ref:`该初稿 `，在代码中基本按照公式推导的计算逻辑，在不过多损失计算效率的前提下进行适当优化，保证了代码的可读性，希望能给学习过程中的你提供一些参考，解答一些疑惑。仅是我浅薄的理解，仅作参考，欢迎指正。 +\ + ++ **特殊情况下的计算优化**。目前包括两个部分： +\ + + - **当震中距** :math:`r` **很大时**。由于结果中Bessel函数的形式为 :math:`J_m(kr)`，当震中距变大时，积分要求波数 :math:`k` 的积分间隔就越小，导致计算变慢。目前 **PyGRT** 实现了 **基于线性插值的Filon积分** 以缓解此问题 :ref:`(纪晨, 姚振兴, 1995) ` :ref:`(初稿) ` ，将每个积分小区间内的Bessel函数取渐近公式，使每个区间内的积分近似值可获得解析解，这样可加大波数 :math:`k` 的积分间隔，提高计算速度。 +\ + + - **当震源深度和台站深度很接近时（如第一类Lamb问题，台站和震源均位于地表）**。此时核函数随着波数 :math:`k` 的增加收敛非常慢，会在收敛值上下波动，导致很难到达到指定的收敛条件，需要耗费更多时间以达到更高的积分上限（不过使用多线程并行优化后，这个问题不再那么刺眼）。目前 **PyGRT** 实现了 **峰谷平均法** 以缓解此问题 :ref:`(Zhang et al., 2003) ` :ref:`(张海明, 2021) `，当波数 :math:`k` 达到一定上限时，可以统计波动的波峰波谷值，再递归取缩减序列 :math:`M_i \leftarrow 0.5\times(M_i + M_{i+1})` 得到收敛值。 +\ + ++ **积分过程可输出记录文件**。通过函数参数设置，在计算格林函数频谱时可输出不同频率的积分过程文件，记录了核函数和被积函数随波数 :math:`k` 的变化，且 **PyGRT** 包含了读取该文件的函数，可以简易地观察到相关量的变化，方便各种测试以及学习过程中的理解。 + + +-------------------------------------------------------------------------- + + +以下包括 **Python** 和 **C** 两部分的API介绍。如果只关心程序的使用，可以只查看 **Python** 部分；若对具体实现过程感兴趣，尤其是在学习 **广义反射透射系数矩阵、离散波数法、线性Filon积分法、峰谷平均法** 等内容感到迷惑难解时，可以在运行程序的同时结合阅读相关 **C代码**，我在C代码中给出了 **明确的中文注释以及相关公式索引** ，希望能帮助你在学习相关内容时解答一些疑惑。 + + + +.. toctree:: + :maxdepth: 1 + :caption: 目录 + :numbered: 3 + + + 文档首页 + install + Tutorial/index + API/api + copyright + +-------------------------------------------------------------------------- + + +主要参考 +--------- + +.. _bouchon_1981: + +.. [1] Michel Bouchon. 1981. A simple method to calculate Green's functions for elastic layered media. Bulletin of the Seismological Society of America. 71(4). 959–971. doi: `10.1785/BSSA0710040959 `_ + +.. _jichen_1995: + +.. [2] 纪晨, 姚振兴, 1995. 区域地震范围的宽频带理论地震图算法研究[J]. 地球物理学报(4): 460-468. + +.. _xie&yao_1989: + +.. [3] 谢小碧, 姚振兴, 1989. 计算分层介质中位错点源静态位移场的广义反射、透射系数矩阵和离散波数方法[J]. 地球物理学报(3): 270-280. + +.. _yao_init_manuscripts: + +.. [4] 姚振兴, 谢小碧. 2022/03. 理论地震图及其应用（初稿）. + +.. _yao&harkrider_1983: + +.. [5] Yao Z. X. and D. G. Harkrider. 1983. A generalized refelection-transmission coefficient matrix and discrete wavenumber method for synthetic seismograms. BSSA. 73(6). 1685-1699. doi: `10.1785/BSSA07306A1685 `_ + +.. _zhang_book_2021: + +.. [6] 张海明. 2021. 地震学中的Lamb问题（上）. 科学出版社. + +.. _zhang_2003: + +.. [7] Zhang, H. M., Chen, X. F., and Chang, S. 2003. An efficient numerical method for computing synthetic seismograms for a layered half-space with sources and receivers at close or same depths. Seismic motion, lithospheric structures, earthquake and volcanic sources: The Keiiti Aki volume, 467-486. doi: `10.1007/978-3-0348-8010-7_3 `_ + diff --git a/docs/source/install.rst b/docs/source/install.rst new file mode 100755 index 00000000..a9b57e80 --- /dev/null +++ b/docs/source/install.rst @@ -0,0 +1,66 @@ + +.. _install_section: + +安装 +============= + +:Author: Zhu Dengda +:Email: zhudengda@mail.iggcas.ac.cn + +----------------------------------------------------------- + + +.. |gr| replace:: `Github Releases `__ + +依赖 +------------ + ++ 如果你想使用C程序，你需要 `Seismic Analysis Code (SAC) `_ ，因为C程序以二进制的SAC格式输出波形。 ++ 如果你想以Python脚本形式使用，依赖已在 :file:`setup.py` 中写好，直接使用 :command:`pip` 安装即可。 + + +免编译安装 +------------ + +目前 **PyGRT** 已分发内置预编译二进制文件的安装包，用户仅需运行以下命令即可 + +.. code-block:: bash + + pip install pygrt-kit + +或者从 |gr| 中下载最新版本的程序压缩包（符合自己的操作系统），解压，在根目录下运行以下命令即可 + +.. code-block:: bash + + pip install . + +如果你不想使用Python，只想使用传统的命令行形式运行C程序，也可从 |gr| 中下载最新版本的程序压缩包（符合自己的操作系统），解压，其中 :rst:dir:`pygrt/C_extension/bin` 和 :rst:dir:`pygrt/C_extension/lib` 为预编译好的可执行文件目录和动态/静态库目录，按自己习惯配置环境变量 :envvar:`PATH` 即可。 + + + +环境变量配置 +------------- +如果你使用 :command:`pip` 安装后，想使用其中的C程序（ :command:`grt` , :command:`grt.syn` etc. ），需配置环境变量 :envvar:`PATH` 。运行以下命令 + +.. code-block:: bash + + python -m pygrt.print + +输出 + +.. code-block:: text + + PyGRT installation directory: + PyGRT executable file directory: + PyGRT library directory: + +将其中的“PyGRT executable file directory”路径添加到环境变量 :envvar:`PATH` 中即可。C程序的运行独立于Python，每个C程序可使用 ``-h`` 查看帮助。 + + +常见问题 +------------ +如果运行报错，提示缺少依赖（常见于MacOS），这通常是缺少 ``OpenMP`` 库。尝试安装 ``gcc`` 编译器，其中会自带 ``OpenMP``。 + + + + diff --git a/example/compare_results/plot_cps_pygrt.py b/example/compare_results/plot_cps_pygrt.py index c6eb018a..f464e605 100644 --- a/example/compare_results/plot_cps_pygrt.py +++ b/example/compare_results/plot_cps_pygrt.py @@ -18,7 +18,7 @@ zeta = 0.8 # compute Green's Functions -st_grt = pymod.gen_gf_spectra( +st_grt = pymod.compute_grn( distarr=rs, nt=nt, dt=dt, diff --git a/example/far_field/plot_compare_pygrt.py b/example/far_field/plot_compare_pygrt.py index 6144391e..3edc0f58 100644 --- a/example/far_field/plot_compare_pygrt.py +++ b/example/far_field/plot_compare_pygrt.py @@ -20,7 +20,7 @@ L2 = -5 T0 = 190 -st_grn1 = pymod.gen_gf_spectra( +st_grn1 = pymod.compute_grn( distarr=rs, nt=nt, dt=dt, @@ -29,7 +29,7 @@ delayT0=T0 )[0] -st_grn2 = pymod.gen_gf_spectra( +st_grn2 = pymod.compute_grn( distarr=rs, nt=nt, dt=dt, diff --git a/example/lamb_problem/run2.py b/example/lamb_problem/run2.py index 35aa69f5..71821086 100644 --- a/example/lamb_problem/run2.py +++ b/example/lamb_problem/run2.py @@ -26,7 +26,7 @@ pymod = pygrt.PyModel1D(modarr, depsrc, deprcv) # compute Green's Functions -st_grn = pymod.gen_gf_spectra( +st_grn = pymod.compute_grn( distarr=rs, nt=nt, dt=dt, diff --git a/example/static_disp/run_stgrt.sh b/example/static_disp/run_stgrt.sh index 0b50e7a5..e78c0a40 100755 --- a/example/static_disp/run_stgrt.sh +++ b/example/static_disp/run_stgrt.sh @@ -31,7 +31,7 @@ gmt begin disp_dc png E300 gmt grdvector synE.nc synN.nc -Q0.1c+e+jc+h1+gblack -Si0.3c gmt meca -Sa0.5c < ∂x - * 1/r*∂t ∂y + * | | uz | ur | ut | + * |----------|-----------|-----------|-----------| + * | ∂z | | | | + * | ∂r | | | | + * | 1/r*∂t | | | | + * + * + * | | uz | ux | uy | + * |----------|-----------|-----------|-----------| + * | ∂z | | | | + * | ∂x | | | | + * | ∂y | | | | + * * * * @param theta (in)r轴相对x轴的旋转弧度 diff --git a/pygrt/C_extension/include/common/quadratic.h b/pygrt/C_extension/include/common/quadratic.h index 77d106d7..c4085c40 100755 --- a/pygrt/C_extension/include/common/quadratic.h +++ b/pygrt/C_extension/include/common/quadratic.h @@ -42,7 +42,7 @@ MYCOMPLEX quad_eval(MYREAL x, MYCOMPLEX a, MYCOMPLEX b, MYCOMPLEX c); /** - * 给定x，根据a,b,c值，估计 \f$ \MYINT_{x_1}^{x_2} f(s)ds \f$ + * 给定x，根据a,b,c值，估计 \f$ \int_{x_1}^{x_2} f(s)ds \f$ * * @param x1 (in)积分下限 * @param x2 (in)积分上限 diff --git a/pygrt/C_extension/include/dynamic/grt.h b/pygrt/C_extension/include/dynamic/grt.h index 42c51bfc..02cd2745 100755 --- a/pygrt/C_extension/include/dynamic/grt.h +++ b/pygrt/C_extension/include/dynamic/grt.h @@ -34,7 +34,7 @@ void set_num_threads(int num_threads); * @param freqs (in)所有频点的频率值（包括未计算的） * @param nr (in)震中距数量 * @param rs (in)震中距数组 - * @param wI (in)虚频率, \f$ \tilde{\omega} =\omega - i\wI \f$ + * @param wI (in)虚频率, \f$ \tilde{\omega} =\omega - i \omega_I \f$ * @param vmin_ref (in)参考最小速度，用于定义波数积分的上限 * @param keps (in)波数积分的收敛条件，要求在某震中距下所有格林函数都收敛，为负数代表不提前判断收敛，按照波数积分上限进行积分 * @param ampk (in)影响波数k积分上限的系数，见下方 @@ -113,7 +113,7 @@ void integ_grn_spec_in_C( * @param freqs (in)所有频点的频率值（包括未计算的） * @param nr (in)震中距数量 * @param rs (in)震中距数组 - * @param wI (in)虚频率, \f$ \tilde{\omega} =\omega - i\wI \f$ + * @param wI (in)虚频率, \f$ \tilde{\omega} =\omega - i \omega_I \f$ * @param vmin_ref (in)参考最小速度，用于定义波数积分的上限 * @param keps (in)波数积分的收敛条件，要求在某震中距下所有格林函数都收敛，为负数代表不提前判断收敛，按照波数积分上限进行积分 * @param ampk (in)影响波数k积分上限的系数，见下方 diff --git a/pygrt/C_extension/include/dynamic/propagate.h b/pygrt/C_extension/include/dynamic/propagate.h index c98c3624..e0ee4a39 100755 --- a/pygrt/C_extension/include/dynamic/propagate.h +++ b/pygrt/C_extension/include/dynamic/propagate.h @@ -71,7 +71,7 @@ * 对应的我们使用ZR表示(z1+, zR+)的效应，FR和ZR也满足类似的递推关系。 * @note 从公式推导上，例如RD_RS，描述的是(zR+, zS-)的效应，但由于我们假定 * 震源位于介质层内，则z=zS并不是介质的物理分界面，此时 \f$ D_{j-1}^{-1} * D_j = I \f$， - * 故在程序可更方便的编写。 + * 故在程序可更方便的编写。（这个在静态情况下不成立，不能以此优化） * @note 接收点位于自由表面的情况 不再单独考虑，合并在接受点浅于震源的情况 * * diff --git a/pygrt/C_extension/include/dynamic/signals.h b/pygrt/C_extension/include/dynamic/signals.h index 4f2ae649..402aa368 100644 --- a/pygrt/C_extension/include/dynamic/signals.h +++ b/pygrt/C_extension/include/dynamic/signals.h @@ -119,6 +119,7 @@ float * get_parabola_wave(float dt, float *Tlen, int *Nt); /** * 生成梯形波或三角波 * + * @verbatim * ^ * | * | @@ -132,6 +133,8 @@ float * get_parabola_wave(float dt, float *Tlen, int *Nt); * |------+------------------+------+----------------> * O T1 T2 T3 T * + * @endverbatim + * * * @param dt (in)采样间隔 * @param T1 (inout)上坡截止时刻 @@ -148,7 +151,7 @@ float * get_trap_wave(float dt, float *T1, float *T2, float *T3, int *Nt); /** * 生成雷克子波 * - * \f[ f(t)=(1-2 \pi^2 \f_0^2 (t-t_0)^2 ) \exp{ - \pi^2 \f_0^2 (t-t_0)^2} \f] + * \f[ f(t)=(1-2 \pi^2 f_0^2 (t-t_0)^2 ) e^{ - \pi^2 f_0^2 (t-t_0)^2} \f] * * @param dt (in)采样间隔 * @param f0 (in)主频 diff --git a/pygrt/C_extension/include/static/static_layer.h b/pygrt/C_extension/include/static/static_layer.h index 4397c5ec..3e6f4ee4 100644 --- a/pygrt/C_extension/include/static/static_layer.h +++ b/pygrt/C_extension/include/static/static_layer.h @@ -20,7 +20,7 @@ /** * 计算自由表面的静态反射系数，公式(6.3.12) * - * @param delta1 (in)表层的 \f$ \Delta \f$ + * @param delta1 (in)表层的 \f$ \Delta = \frac{\lambda + \mu}{\lambda + 3\mu} \f$ * @param R_tilt[2][2] (out)P-SV系数矩阵，SH系数为1 * */ diff --git a/pygrt/C_extension/include/static/static_source.h b/pygrt/C_extension/include/static/static_source.h index 6aa98ef3..1cfc4b68 100644 --- a/pygrt/C_extension/include/static/static_source.h +++ b/pygrt/C_extension/include/static/static_source.h @@ -3,7 +3,7 @@ * @author Zhu Dengda (zhudengda@mail.iggcas.ac.cn) * @date 2025-02-18 * - * 以下代码实现的是 静态震源系数————剪切源， 参考： + * 以下代码实现的是 静态震源系数————爆炸源，垂直力源，水平力源，剪切源， 参考： * 1. 谢小碧, 姚振兴, 1989. 计算分层介质中位错点源静态位移场的广义反射、 * 透射系数矩阵和离散波数方法[J]. 地球物理学报(3): 270-280. * diff --git a/pygrt/C_extension/src/common/model.c b/pygrt/C_extension/src/common/model.c index 9d59f0ce..30cc605b 100755 --- a/pygrt/C_extension/src/common/model.c +++ b/pygrt/C_extension/src/common/model.c @@ -288,6 +288,9 @@ PYMODEL1D * read_pymod_from_file(const char *command, const char *modelpath, dou while(fgets(line, sizeof(line), fp)) { iline++; + + // 注释行 + if(line[0]=='#') continue; h = va = vb = rho = qa = qb = -9.0; if(ncols != sscanf(line, "%lf %lf %lf %lf %lf %lf\n", &h, &va, &vb, &rho, &qa, &qb)){ diff --git a/pygrt/C_extension/src/dynamic/grt_syn.c b/pygrt/C_extension/src/dynamic/grt_syn.c index e71986c2..10a3888b 100644 --- a/pygrt/C_extension/src/dynamic/grt_syn.c +++ b/pygrt/C_extension/src/dynamic/grt_syn.c @@ -162,6 +162,7 @@ printf("\n" "\n" " -T/////\n" " Six elements of Moment Tensor. \n" +" x (North), y (East), z (Downward).\n" " Notice they will be scaled by .\n" "\n" " -F//\n" diff --git a/pygrt/C_extension/src/static/stgrt_syn.c b/pygrt/C_extension/src/static/stgrt_syn.c index 4968d905..dec43e10 100644 --- a/pygrt/C_extension/src/static/stgrt_syn.c +++ b/pygrt/C_extension/src/static/stgrt_syn.c @@ -103,6 +103,7 @@ printf("\n" "\n" " -T/////\n" " Six elements of Moment Tensor. \n" +" x (North), y (East), z (Downward).\n" " Notice they will be scaled by .\n" "\n" " -F//\n" diff --git a/pygrt/pymod.py b/pygrt/pymod.py index 49cddee4..a9af0f9b 100755 --- a/pygrt/pymod.py +++ b/pygrt/pymod.py @@ -124,6 +124,15 @@ def __init__(self, modarr0:np.ndarray, depsrc:float, deprcv:float): def compute_travt1d(self, dist:float): + r""" + 调用C程序，计算初至P波和S波的走时 + + :param dist: 震中距 + + :return: + - **travtP** - 初至P波走时(s) + - **travtS** - 初至S波走时(s) + """ travtP = C_compute_travt1d( self.c_pymod1d.Thk, self.c_pymod1d.Va, @@ -240,7 +249,7 @@ def compute_grn( :param nt: 时间点数，借助于 `SciPy`，nt不再要求是2的幂次 :param dt: 采样间隔(s) :param freqband: 频率范围(Hz)，以此确定待计算的离散频率点 - :param zeta: 定义虚频率的系数 :math:`\zeta`， 虚频率 :math:`\tilde{\omega} = \omega - j*w_I, w_I = \zeta*\pi/T, T=nt*dt` , T为时窗长度。 + :param zeta: 定义虚频率的系数 :math:`\zeta` ， 虚频率 :math:`\tilde{\omega} = \omega - j*w_I, w_I = \zeta*\pi/T, T=nt*dt` , T为时窗长度。 使用离散波数积分时为了避开附加源以及奇点的影响， :ref:`(Bouchon, 1981) ` 在频率上添加微小虚部， 更多测试见 :ref:`(张海明, 2021) ` :param vmin_ref: 最小参考速度，默认vmin=max(minimum velocity, 0.1)，用于定义波数积分上限，小于0则在达到积分上限后使用峰谷平均法 @@ -248,8 +257,8 @@ def compute_grn( :param keps: 波数k积分收敛条件，见 :ref:`(Yao and Harkrider, 1983) ` :ref:`(初稿) `， 为负数代表不提前判断收敛，按照波数积分上限进行积分 :param ampk: 影响波数k积分上限的系数，见下方 - :param iwk0: k0是否取随频率变化的线性关系，即 :math:` k_{0} = k_{0} * f/f_{max}` - :param k0: 波数k积分的上限 :math:`\tilde{k_{max}}=\sqrt{(k_{0}*\pi/hs)^2 + (ampk*w/vmin_{ref})^2} ` , 波数k积分循环必须退出, hs=max(震源和台站深度差,1.0) + :param iwk0: k0是否取随频率变化的线性关系，即 :math:`k_{0} = k_{0} * f/f_{max}` + :param k0: 波数k积分的上限 :math:`\tilde{k_{max}}=\sqrt{(k_{0}*\pi/hs)^2 + (ampk*w/vmin_{ref})^2}` , 波数k积分循环必须退出, hs=max(震源和台站深度差,1.0) :param Length: 定义波数k积分的间隔 `dk=2\pi / (L*rmax)`, 选取要求见 :ref:`(Bouchon, 1981) ` :ref:`(张海明, 2021) `，默认自动选择；负数表示使用Filon积分 :param calc_upar: 是否计算位移u的空间导数 @@ -574,7 +583,7 @@ def compute_static_grn( :param statsfile: 波数k积分（包括Filon积分和峰谷平均法）的过程记录文件，常用于debug或者观察积分过程中 :math:`F(k,\omega)` 和 :math:`F(k,\omega)J_m(kr)k` 的变化 :return: - - **dataDct* - 字典形式的格林函数 + - **dataDct** - 字典形式的格林函数 """ depsrc = self.depsrc diff --git a/pygrt/utils.py b/pygrt/utils.py index e5adbca3..d75844a3 100755 --- a/pygrt/utils.py +++ b/pygrt/utils.py @@ -515,7 +515,7 @@ def gen_syn_from_gf_MT(st:Union[Stream,dict], M0:float, MT:ArrayLike, az:float=- :param st: 计算好的时域格林函数, :class:`obspy.Stream` 类型，或者静态格林函数（字典类型） :param M0: 标量地震矩 - :param MT: 矩张量 (M11, M12, M13, M22, M23, M33),下标1,2,3分别代表北向，东向，垂直向上 + :param MT: 矩张量 (M11, M12, M13, M22, M23, M33),下标1,2,3分别代表北向，东向，垂直向下 :param az: 台站方位角，以北顺时针为正，0<=az<=360 （静态情况不需要） :param ZNE: 是否以ZNE分量输出? :param calc_upar: 是否计算位移u的空间导数