FEAT: exact closed-form solution for the Lamb problem of the first kind (#117)

Author: Dengda98

docs/Makefile

@@ -13,8 +13,8 @@ help:
 	@$(SPHINXBUILD) -M help "$(SOURCEDIR)" "$(BUILDDIR)" $(SPHINXOPTS) $(O)
 
 clean:
-	rm $(BUILDDIR)/* -rf
-	rm doxygen_h/* -rf
+	rm $(BUILDDIR)/ -rf
+	rm doxygen_h/ -rf
 
 .PHONY: help Makefile clean
 

pygrt/C_extension/include/grt.h

@@ -41,9 +41,10 @@
     X(static_strain)       \
     X(static_stress)       \
     /* other */            \
-    X(ker2asc)                 \
-    X(sac2asc)                 \
+    X(ker2asc)             \
+    X(sac2asc)             \
     X(travt)               \
+    X(lamb1)               \
 // ------------------------------------------------------
 
 /** 子模块的函数指针类型 */

pygrt/C_extension/include/grt/common/checkerror.h

@@ -13,6 +13,7 @@
 #include <sys/stat.h>
 #include <dirent.h>
 #include <errno.h>
+#include <stdlib.h>
 
 #include "grt/common/colorstr.h"
 

pygrt/C_extension/include/grt/lamb/elliptic.h

@@ -0,0 +1,54 @@
+/**
+ * @file   elliptic.h
+ * @author Zhu Dengda ([email protected])
+ * @date   2025-11
+ * 
+ *    使用渐近解和 Carlson 算法计算三类完全椭圆积分，参考：
+ *    
+ *        1. Abramowitz, Milton, and Irene A. Stegun. 1964. Handbook of Mathematical Functions with Formulas, Graphs, and Mathematical Tables. Vol. 55. US Government printing office.
+ *        2. Carlson, B. C. 1995. “Numerical Computation of Real or Complex Elliptic Integrals.” Numerical Algorithms 10 (1): 13–26. https://doi.org/10.1007/BF02198293.
+ *        3. 张海明, 冯禧 著. 2024. 地震学中的 Lamb 问题（下）. 科学出版社
+ */
+
+#pragma once
+
+#include "grt/common/const.h"
+
+/**
+ * 使用 Abramowitz and Stegun (1964, P591-592) 提供的渐近表达式计算第一类完全椭圆积分， \f$ m \in (0,1) \f$
+ * \f[
+ *  K(m) = \int_0^{\pi/2} \dfrac{1}{\sqrt{1 - m \sin^2\theta}}\ d\theta
+ * \f]  
+ * 
+ * @param[in]     m     参数 m
+ * @return    积分值
+ * 
+ */
+MYREAL grt_ellipticK(const MYREAL m);
+
+
+/**
+ * 使用 Abramowitz and Stegun (1964, P591-592) 提供的渐近表达式计算第二类完全椭圆积分， \f$ m \in (0,1) \f$
+ * \f[
+ *  E(m) = \int_0^{\pi/2} \sqrt{1 - m \sin^2\theta}\ d\theta
+ * \f]  
+ * 
+ * @param[in]     m     参数 m
+ * @return    积分值
+ * 
+ */
+MYREAL grt_ellipticE(const MYREAL m);
+
+
+/**
+ * 使用 Carlson (1995) 算法计算第三类完全椭圆积分， \f$ m \in (0,1) \f$
+ * \f[
+ *  {\it \Pi}(m) = \int_0^{\pi/2} \dfrac{1}{(1-nx^2)\sqrt{1-x^2}\sqrt{1 - m \sin^2\theta}}\ d\theta
+ * \f]  
+ * 
+ * @param[in]     n     参数 n
+ * @param[in]     m     参数 m
+ * @return    积分值
+ * 
+ */
+MYCOMPLEX grt_ellipticPi(const MYCOMPLEX n, const MYREAL m);

pygrt/C_extension/include/grt/lamb/lamb1.h

@@ -0,0 +1,29 @@
+/**
+ * @file   lamb1.h
+ * @author Zhu Dengda ([email protected])
+ * @date   2025-11
+ * 
+ *    使用广义闭合解求解第一类 Lamb 问题，参考：
+ * 
+ *        张海明, 冯禧 著. 2024. 地震学中的 Lamb 问题（下）. 科学出版社
+ */
+
+#pragma once
+
+#include "grt/common/const.h"
+#include "grt/lamb/lamb_util.h"
+
+
+/**
+ * 使用广义闭合解求解第一类 Lamb 问题
+ * 
+ * @param[in]    alpha     P 波速度
+ * @param[in]    beta      S 波速度
+ * @param[in]    ts        归一化时间序列
+ * @param[in]    nt        时间序列点数
+ * @param[in]    azimuth   方位角，单位度
+ * @param[out]   u         记录结果的指针，如果为NULL则输出到标准输出
+ * 
+ */
+void grt_solve_lamb1(
+    const MYREAL alpha, const MYREAL beta, const MYREAL *ts, const int nt, const MYREAL azimuth, MYREAL (*u)[3][3]);

pygrt/C_extension/include/grt/lamb/lamb_util.h

@@ -0,0 +1,44 @@
+/**
+ * @file   lamb_util.h
+ * @author Zhu Dengda ([email protected])
+ * @date   2025-11
+ * 
+ *    一些使用广义闭合解求解 Lamb 问题过程中可能用到的辅助函数
+ */
+
+#pragma once
+
+#include "grt/common/const.h"
+
+/**
+ * 求解一元三次形式的 Rayleigh 方程的根,  \f$ x^3 + ax^2 + bx + c = 0 \f$
+ * 
+ * @param[in]      a     系数 a
+ * @param[in]      b     系数 b
+ * @param[in]      c     系数 c
+ * @param[out]     y3    三个根，其中 y3[2] 为正根
+ */
+void grt_roots3(const MYREAL a, const MYREAL b, const MYREAL c, MYCOMPLEX y3[3]);
+
+/**
+ * 做如下多项式求值， \f$ \sum_{m=0}^n C_{2m+o} y^m \f$
+ * 
+ * @param[in]    C       数组 C
+ * @param[in]    n       最高幂次 n
+ * @param[in]    y       自变量 y
+ * @param[in]    o       偏移量
+ * 
+ * @return    多项式结果
+ * 
+ */
+MYCOMPLEX grt_evalpoly2(const MYCOMPLEX *C, const int n, const MYCOMPLEX y, const int offset);
+
+
+/**
+ * 根据 Vs/Vp 比值，计算泊松比，  \f$ \nu = \dfrac{1}{2} \dfrac{1 - 2k^2}{1 - k^2},  k=\dfrac{V_S}{V_P} \f$
+ * 
+ * @param[in]    k    Vs/Vp
+ * @return   泊松比
+ */
+
+MYREAL grt_possion_ratio(const MYREAL k);