Update stats file (kernel functions value in integration ) format, python reading functions, add grt.k2a command

.github/workflows/testbuild.yml

@@ -3,6 +3,8 @@ name: Just test
 # 触发机制为“非主分支的提交”
 on:
   push:
+    paths:
+      - 'pygrt/**'  # 源码文件夹有提交才触发测试
     branches-ignore:
       - main    # 排除主分支
       - master  # 如果有 master 分支，也可以排除

pygrt/C_extension/include/common/const.h

@@ -122,6 +122,8 @@ typedef int MYINT;  ///< 整数
 #define MIN_DEPTH_GAP_SRC_RCV  1.0  ///< 震源和台站的最小深度差（不做绝对限制，仅用于参考波数积分上限）
 #define GCC_ALWAYS_INLINE __attribute__((always_inline))  ///< gcc编译器不改动内联函数
 
-#define GRT_STATIC_STRING_FMT "%18s"  ///< stgrt系列程序字符串输出格式
-#define GRT_STATIC_REAL_FMT "%18.8e"  ///< stgrt系列程序浮点数输出格式
+#define GRT_STRING_FMT "%18s"  ///< 字符串输出格式
+#define GRT_REAL_FMT "%18.8e"  ///< 浮点数输出格式
+#define GRT_CMPLX_FMT "%18.8e%-+14.8eJ"   ///< 复数输出格式
+#define GRT_STR_CMPLX_FMT "%34s"    ///< 与复数格式同长度的字符串输出格式
 

pygrt/C_extension/include/common/dwm.h

@@ -46,7 +46,7 @@
  * @param  sum_HF_uir_J[nr][3][4]   (out)水平力源
  * @param  sum_DC_uir_J[nr][3][4]   (out)剪切源
  * 
- * @param  fstats[nr])          (out)不同震中距的格林函数积分过程文件
+ * @param  fstats               (out)文件指针，保存不同k值的格林函数积分核函数
  * @param  kerfunc              (in)计算核函数的函数指针
  * 
  * @return  k        积分截至时的波数
@@ -61,4 +61,4 @@ MYREAL discrete_integ(
     MYCOMPLEX sum_HF_uiz_J[nr][3][4],  MYCOMPLEX sum_DC_uiz_J[nr][3][4],  
     MYCOMPLEX sum_EXP_uir_J[nr][3][4], MYCOMPLEX sum_VF_uir_J[nr][3][4],  
     MYCOMPLEX sum_HF_uir_J[nr][3][4],  MYCOMPLEX sum_DC_uir_J[nr][3][4],  
-    FILE *(fstats[nr]), KernelFunc kerfunc);
+    FILE *fstats, KernelFunc kerfunc);

pygrt/C_extension/include/common/fim.h

@@ -52,7 +52,7 @@
  * @param  sum_HF_uir_J[nr][3][4]   (out)水平力源
  * @param  sum_DC_uir_J[nr][3][4]   (out)剪切源
  * 
- * @param  fstats[nr]           (out)不同震中距的格林函数积分过程文件
+ * @param  fstats               (out)文件指针，保存不同k值的格林函数积分核函数
  * @param  kerfunc              (in)计算核函数的函数指针
  * 
  * @return  k        积分截至时的波数
@@ -67,6 +67,6 @@ MYREAL linear_filon_integ(
     MYCOMPLEX sum_HF_uiz_J[nr][3][4],  MYCOMPLEX sum_DC_uiz_J[nr][3][4],  
     MYCOMPLEX sum_EXP_uir_J[nr][3][4], MYCOMPLEX sum_VF_uir_J[nr][3][4],  
     MYCOMPLEX sum_HF_uir_J[nr][3][4],  MYCOMPLEX sum_DC_uir_J[nr][3][4],  
-    FILE *(fstats[nr]), KernelFunc kerfunc);
+    FILE *fstats, KernelFunc kerfunc);
 
 

pygrt/C_extension/include/common/iostats.h

@@ -26,10 +26,6 @@
  * @param    VF_qwv[3][3]     (in)垂直力源核函数
  * @param    HF_qwv[3][3]     (in)水平力源核函数
  * @param    DC_qwv[3][3]     (in)剪切源核函数
- * @param    EXP_J[3][4]      (in)爆炸源，该dk区间内的积分值，下同
- * @param    VF_J[3][4]       (in)垂直力源
- * @param    HF_J[3][4]       (in)水平力源
- * @param    DC_J[3][4]       (in)剪切源
  * 
  * 
  * @note     文件记录的值均为波数积分的中间结果，与最终的结果还差一系列的系数，
@@ -38,9 +34,10 @@
 void write_stats(
     FILE *f0, MYREAL k, 
     const MYCOMPLEX EXP_qwv[3][3], const MYCOMPLEX VF_qwv[3][3], 
-    const MYCOMPLEX HF_qwv[3][3],  const MYCOMPLEX DC_qwv[3][3], 
-    const MYCOMPLEX EXP_J[3][4], const MYCOMPLEX VF_J[3][4], 
-    const MYCOMPLEX HF_J[3][4],  const MYCOMPLEX DC_J[3][4]);
+    const MYCOMPLEX HF_qwv[3][3],  const MYCOMPLEX DC_qwv[3][3]
+    // const MYCOMPLEX EXP_J[3][4], const MYCOMPLEX VF_J[3][4], 
+    // const MYCOMPLEX HF_J[3][4],  const MYCOMPLEX DC_J[3][4]
+);
 
 
 

pygrt/C_extension/include/common/ptam.h

@@ -34,8 +34,6 @@
  * @param    k0        (in)先前的积分已经进行到了波数k0
  * @param    predk     (in)先前的积分使用的积分间隔dk，因为峰谷平均法使用的
  *                     积分间隔会和之前的不一致，这里传入该系数以做预先调整
- * @param    rmin      (in)最小震中距 
- * @param    rmax      (in)最大震中距 
  * @param    omega     (in)复数频率 
  * @param    nr        (in)震中距数量
  * @param    rs        (in)震中距数组  
@@ -54,23 +52,22 @@
  * @param  sum_HF_uir_J0[nr][3][4]   (out)水平力源
  * @param  sum_DC_uir_J0[nr][3][4]   (out)剪切源
  * 
- * @param    fstats[nr]          (out)波数积分过程文件指针
- * @param    ptam_fstats[nr]     (out)峰谷平均法过程文件指针
+ * @param    ptam_fstatsnr        (out)峰谷平均法过程文件指针数组
  * @param    kerfunc              (in)计算核函数的函数指针
  * 
  * 
  */
 void PTA_method(
-    const MODEL1D *mod1d, MYREAL k0, MYREAL predk, MYREAL rmin, MYREAL rmax, MYCOMPLEX omega, 
-    MYINT nr, MYREAL *rs, 
+    const MODEL1D *mod1d, MYREAL k0, MYREAL predk, MYCOMPLEX omega, 
+    MYINT nr, MYREAL *rs,
     MYCOMPLEX sum_EXP_J0[nr][3][4], MYCOMPLEX sum_VF_J0[nr][3][4],  
     MYCOMPLEX sum_HF_J0[nr][3][4],  MYCOMPLEX sum_DC_J0[nr][3][4],  
     bool calc_upar,
     MYCOMPLEX sum_EXP_uiz_J0[nr][3][4], MYCOMPLEX sum_VF_uiz_J0[nr][3][4],  
     MYCOMPLEX sum_HF_uiz_J0[nr][3][4],  MYCOMPLEX sum_DC_uiz_J0[nr][3][4],  
     MYCOMPLEX sum_EXP_uir_J0[nr][3][4], MYCOMPLEX sum_VF_uir_J0[nr][3][4],  
     MYCOMPLEX sum_HF_uir_J0[nr][3][4],  MYCOMPLEX sum_DC_uir_J0[nr][3][4], 
-    FILE *(fstats[nr]), FILE *(ptam_fstats[nr]), KernelFunc kerfunc);
+    FILE *ptam_fstatsnr[nr][2], KernelFunc kerfunc);
 
 
 

pygrt/C_extension/src/common/Makefile

@@ -3,14 +3,20 @@ CC := gcc
 CFLAGS := -Wall -g -fPIC -I../../include -lm
 
 BUILD_DIR = ../../build/common
+BIN_DIR = ../../bin
 SRCS := $(wildcard *.c)
 OBJS := $(patsubst %.c, $(BUILD_DIR)/%.o, $(SRCS))
 DEPS := $(OBJS:.o=.d)
 
-progs: 
-	@echo "No rule to make target 'progs'"
+OBJS := $(filter-out \
+	$(BUILD_DIR)/grt_k2a.o \
+	, $(OBJS))
 
-all: objs
+PROGS := $(BIN_DIR)/grt.k2a
+
+progs: $(OBJS) $(BIN_DIR) $(PROGS)
+
+all: objs progs
 
 objs: $(BUILD_DIR) $(OBJS)
 
@@ -21,6 +27,13 @@ $(BUILD_DIR):
 $(BUILD_DIR)/%.o: %.c
 	$(CC) -o $@  -c $< $(CFLAGS) 
 
+# ------------------------Executable files------------------------
+$(BIN_DIR):
+	@mkdir -p $(BIN_DIR)
+
+$(BIN_DIR)/grt.k2a: grt_k2a.c $(OBJS)
+	$(CC) -o $@ $^  $(CFLAGS)
+
 
 # ----------------------- Dependency generation -----------------------
 -include $(DEPS)
@@ -33,7 +46,8 @@ $(BUILD_DIR)/%.d: %.c
 
 # ---------------------------------------------------------------------
 
-cleanbuild: clean
+cleanbuild:
+	rm -rf $(BUILD_DIR)
 
-clean:
-	rm -rf $(BUILD_DIR)
+clean: cleanbuild
+	rm -f $(PROGS)

pygrt/C_extension/src/common/dwm.c

@@ -33,7 +33,7 @@ MYREAL discrete_integ(
     MYCOMPLEX sum_HF_uiz_J[nr][3][4],  MYCOMPLEX sum_DC_uiz_J[nr][3][4],  
     MYCOMPLEX sum_EXP_uir_J[nr][3][4], MYCOMPLEX sum_VF_uir_J[nr][3][4],  
     MYCOMPLEX sum_HF_uir_J[nr][3][4],  MYCOMPLEX sum_DC_uir_J[nr][3][4],  
-    FILE *(fstats[nr]), KernelFunc kerfunc)
+    FILE *fstats, KernelFunc kerfunc)
 {
     MYCOMPLEX EXP_J[3][4], VF_J[3][4], HF_J[3][4],  DC_J[3][4];
 
@@ -67,15 +67,21 @@ MYREAL discrete_integ(
 
     // 波数k循环 (5.9.2)
     while(true){
-        k += dk; 
-
+
         if(k > kmax)  break;
+        k += dk; 
 
         // printf("w=%15.5e, ik=%d\n", CREAL(omega), ik);
         // 计算核函数 F(k, w)
         kerfunc(mod1d, omega, k, pEXP_qwv, pVF_qwv, pHF_qwv, pDC_qwv, 
                 calc_upar, pEXP_uiz_qwv, pVF_uiz_qwv, pHF_uiz_qwv, pDC_uiz_qwv); 
-
+
+        // 记录积分核函数
+        if(fstats!=NULL){
+            write_stats(
+                fstats, k, 
+                EXP_qwv, VF_qwv, HF_qwv, DC_qwv);
+        }
 
         // 震中距rs循环
         iendk = true;
@@ -92,15 +98,6 @@ MYREAL discrete_integ(
             int_Pk(k, rs[ir], 
                    pEXP_qwv, pVF_qwv, pHF_qwv, pDC_qwv, false,
                    EXP_J, VF_J, HF_J, DC_J);
-
-            // 记录积分结果
-            if(fstats[ir]!=NULL){
-                write_stats(
-                    fstats[ir], k, 
-                    EXP_qwv, VF_qwv, HF_qwv, DC_qwv,
-                    EXP_J, VF_J, HF_J, DC_J);
-            }
-
 
             iendk0 = true;
             for(MYINT m=0; m<3; ++m){

pygrt/C_extension/src/common/fim.c

@@ -32,7 +32,7 @@ MYREAL linear_filon_integ(
     MYCOMPLEX sum_HF_uiz_J[nr][3][4],  MYCOMPLEX sum_DC_uiz_J[nr][3][4],  
     MYCOMPLEX sum_EXP_uir_J[nr][3][4], MYCOMPLEX sum_VF_uir_J[nr][3][4],  
     MYCOMPLEX sum_HF_uir_J[nr][3][4],  MYCOMPLEX sum_DC_uir_J[nr][3][4],  
-    FILE *(fstats[nr]), KernelFunc kerfunc)
+    FILE *fstats, KernelFunc kerfunc)
 {   
     for(MYINT ir=0; ir<nr; ++ir){
         for(MYINT m=0; m<3; ++m){
@@ -102,14 +102,21 @@ MYREAL linear_filon_integ(
     // k循环 
     ik = 0;
     while(true){
-        k += dk; 
-
+
         if(k > kmax) break;
+        k += dk; 
 
         // 计算核函数 F(k, w)
         kerfunc(mod1d, omega, k, pEXP_qwv, pVF_qwv, pHF_qwv, pDC_qwv,
                 calc_upar, pEXP_uiz_qwv, pVF_uiz_qwv, pHF_uiz_qwv, pDC_uiz_qwv); 
 
+        // 记录积分结果
+        if(fstats!=NULL){
+            write_stats(
+                fstats, k, 
+                pEXP_qwv, pVF_qwv, pHF_qwv, pDC_qwv);
+        }
+
         // 震中距rs循环
         iendk = true;
         for(MYINT ir=0; ir<nr; ++ir){
@@ -121,14 +128,6 @@ MYREAL linear_filon_integ(
                 pEXP_qwv, pVF_qwv, pHF_qwv, pDC_qwv, false,
                 EXP_J, VF_J, HF_J, DC_J);
 
-            // 记录积分结果
-            if(fstats[ir]!=NULL){
-                write_stats(
-                    fstats[ir], k, 
-                    pEXP_qwv, pVF_qwv, pHF_qwv, pDC_qwv,
-                    EXP_J, VF_J, HF_J, DC_J);
-            }
-
 
             iendk0 = true;
             for(MYINT m=0; m<3; ++m){