Support stress and strain calculation

pygrt/C_extension/include/common/attenuation.h

@@ -16,6 +16,7 @@
  * \f[
  *  c(\omega) = c(2\pi)\times (1 + \frac{1}{\pi Q} \log(\frac{\omega}{2\pi}) + \frac{i}{2Q}) 
  * \f] 
+ * 其中虚数部分的正负号和书中不同，是因为书中使用的傅里叶变换的e指数符号和我们通常使用的相反。
  * 
  * @param Qinv     (in) 1/Q
  * @param omega    (in)复数频率\f$ \tilde{\omega} =\omega - i\zeta \f$ 

pygrt/C_extension/include/common/sacio2.h

@@ -0,0 +1,34 @@
+/**
+ * @file   sacio2.h
+ * @author Zhu Dengda (zhudengda@mail.iggcas.ac.cn)
+ * @date   2025-03-31
+ * 
+ *    在已有的sacio基础上进行部分函数的封装
+ * 
+ */
+
+#pragma once 
+
+#include "sacio.h"
+
+/**
+ * 读取SAC头段变量
+ * 
+ * @param      command        (in)当前程序命令名称
+ * @param      name           (in)SAC文件路径
+ * @param      hd             (out)SAC头段变量结构体
+ */
+void read_SAC_HEAD(const char *command, const char *name, SACHEAD *hd);
+
+
+/**
+ * 读取SAC文件
+ * 
+ * @param      command        (in)当前程序命令名称
+ * @param      name           (in)SAC文件路径
+ * @param      hd             (out)SAC头段变量结构体
+ * @param      arrout         (inout)预分配内存，不需要则设为NULL
+ * 
+ * @return     浮点数指针
+ */
+float * read_SAC(const char *command, const char *name, SACHEAD *hd, float *arrout);

pygrt/C_extension/src/common/sacio2.c

@@ -0,0 +1,41 @@
+/**
+ * @file   sacio2.c
+ * @author Zhu Dengda (zhudengda@mail.iggcas.ac.cn)
+ * @date   2025-03-31
+ * 
+ */
+
+
+
+#include <stdio.h>
+#include <stdlib.h>
+
+#include "common/sacio2.h"
+#include "common/sacio.h"
+#include "common/colorstr.h"
+
+
+void read_SAC_HEAD(const char *command, const char *name, SACHEAD *hd){
+    int lswap = read_sac_head(name, hd);
+    if(lswap == -1){
+        fprintf(stderr, "[%s] " BOLD_RED "read %s head failed.\n" DEFAULT_RESTORE, command, name);
+        exit(EXIT_FAILURE);
+    }
+}
+
+
+float * read_SAC(const char *command, const char *name, SACHEAD *hd, float *arrout){
+    float *arrin=NULL;
+    if((arrin = read_sac(name, hd)) == NULL){
+        fprintf(stderr, "[%s] " BOLD_RED "read %s failed.\n" DEFAULT_RESTORE, command, name);
+        exit(EXIT_FAILURE);
+    }
+
+    if(arrout!=NULL){
+        for(int i=0; i<hd->npts; ++i)  arrout[i] = arrin[i];
+        free(arrin);
+        arrin = arrout;
+    }
+
+    return arrin;
+}

pygrt/C_extension/src/dynamic/Makefile

@@ -16,11 +16,15 @@ OBJS := $(filter-out \
 	$(BUILD_DIR)/grt_main.o \
 	$(BUILD_DIR)/grt_syn.o \
 	$(BUILD_DIR)/grt_b2a.o \
+	$(BUILD_DIR)/grt_strain.o \
+	$(BUILD_DIR)/grt_stress.o \
 	, $(OBJS))
 
 PROGS := $(BIN_DIR)/grt \
          $(BIN_DIR)/grt.syn \
-		 $(BIN_DIR)/grt.b2a
+		 $(BIN_DIR)/grt.b2a \
+		 $(BIN_DIR)/grt.strain \
+		 $(BIN_DIR)/grt.stress
 
 all: objs progs
 
@@ -40,14 +44,20 @@ $(BUILD_DIR)/%.o: %.c
 $(BIN_DIR):
 	@mkdir -p $(BIN_DIR)
 
-$(BIN_DIR)/grt: grt_main.c $(OBJS)
-	$(CC) -o $@ $^  $(COMMON_OBJS) $(TRAVT_OBJS)  $(CFLAGS)
+$(BIN_DIR)/grt: grt_main.c $(OBJS) $(COMMON_OBJS) $(TRAVT_OBJS)
+	$(CC) -o $@ $^ $(CFLAGS)
 
-$(BIN_DIR)/grt.syn: grt_syn.c $(OBJS)
-	$(CC) -o $@ $^  $(COMMON_OBJS) $(TRAVT_OBJS)  $(CFLAGS)
+$(BIN_DIR)/grt.syn: grt_syn.c $(OBJS) $(COMMON_OBJS) $(TRAVT_OBJS)
+	$(CC) -o $@ $^  $(CFLAGS)
 
-$(BIN_DIR)/grt.b2a: grt_b2a.c
-	$(CC) -o $@ $^  $(COMMON_OBJS) $(CFLAGS)
+$(BIN_DIR)/grt.b2a: grt_b2a.c $(COMMON_OBJS)
+	$(CC) -o $@ $^  $(CFLAGS)
+
+$(BIN_DIR)/grt.strain: grt_strain.c $(COMMON_OBJS)
+	$(CC) -o $@ $^  $(CFLAGS)
+
+$(BIN_DIR)/grt.stress: grt_stress.c $(COMMON_OBJS)
+	$(CC) -o $@ $^  $(CFLAGS)
 
 # ----------------------- Dependency generation -----------------------
 -include $(DEPS)

pygrt/C_extension/src/dynamic/grt_b2a.c

@@ -14,7 +14,7 @@
 #include <stdlib.h>
 #include <string.h>
 
-#include "common/sacio.h"
+#include "common/sacio2.h"
 #include "common/logo.h"
 #include "common/colorstr.h"
 
@@ -71,11 +71,7 @@ int main(int argc, char **argv){
 
     // 读入SAC文件
     SACHEAD hd;
-    float *arr=NULL;
-    if((arr = read_sac(filepath, &hd)) == NULL){
-        fprintf(stderr, "[%s] " BOLD_RED "read %s failed.\n" DEFAULT_RESTORE, command, filepath);
-        exit(EXIT_FAILURE);
-    }
+    float *arr = read_SAC(command, filepath, &hd, NULL);
 
     // 将波形写入标准输出，第一列时间，第二列振幅
     float begt = hd.b;

pygrt/C_extension/src/dynamic/grt_main.c

@@ -726,9 +726,6 @@ static void ifft_one_trace(
         float_arr[i] = out[i];
     }
 
-    // 写入虚频率
-    hd->user0 = wI;
-
     write_sac(outpath, *hd, float_arr);
     // FILE *fp = fopen(outpath, "wb");
     // fwrite(out, sizeof(float), nt, fp);
@@ -1143,6 +1140,14 @@ int main(int argc, char **argv) {
     // 记录震源和台站深度
     hd.evdp = depsrc; // km
     hd.stel = (-1.0)*deprcv*1e3; // m
+    // 写入虚频率
+    hd.user0 = wI;
+    // 写入接受点的Vp,Vs,rho
+    hd.user1 = pymod->Va[pymod->ircv];
+    hd.user2 = pymod->Vb[pymod->ircv];
+    hd.user3 = pymod->Rho[pymod->ircv];
+    hd.user4 = RONE/pymod->Qa[pymod->ircv];
+    hd.user5 = RONE/pymod->Qb[pymod->ircv];
 
 
     // 做反傅里叶变换，保存SAC文件

pygrt/C_extension/src/dynamic/grt_strain.c

@@ -0,0 +1,148 @@
+/**
+ * @file   grt_strain.c
+ * @author Zhu Dengda (zhudengda@mail.iggcas.ac.cn)
+ * @date   2025-03-28
+ * 
+ *    根据预先合成的位移空间导数，组合成应变
+ * 
+ */
+
+
+#include <stdio.h>
+#include <unistd.h>
+#include <stdlib.h>
+#include <dirent.h>
+#include <ctype.h>
+#include <string.h>
+
+#include "common/sacio2.h"
+#include "common/const.h"
+#include "common/logo.h"
+#include "common/colorstr.h"
+
+
+
+/**
+ * 打印使用说明
+ */
+static void print_help(){
+print_logo();
+printf("\n"
+"[grt.strain]\n\n"
+"    Conbine spatial derivatives of displacements into strain.\n"
+"\n\n"
+"Usage:\n"
+"----------------------------------------------------------------\n"
+"    grt.strain <syn_dir>/<name>\n"
+"\n\n\n"
+);
+}
+
+
+
+int main(int argc, char **argv){
+    const char *command = argv[0];
+
+    // 输入不够
+    if(argc < 2){
+        fprintf(stderr, "[%s] " BOLD_RED "Error! Need set an input. Use '-h' for help.\n" DEFAULT_RESTORE, command);
+        exit(EXIT_FAILURE);
+    }
+
+    // 输入过多
+    if(argc > 2){
+        fprintf(stderr, "[%s] " BOLD_RED "Error! You should set only one input. Use '-h' for help.\n" DEFAULT_RESTORE, command);
+        exit(EXIT_FAILURE);
+    }
+
+    // 使用-h查看帮助
+    if(strcmp(argv[1], "-h") == 0){
+        print_help();
+        exit(EXIT_SUCCESS);
+    }
+
+
+    // 合成地震图目录路径
+    char *s_synpath = (char*)malloc(sizeof(char)*(strlen(argv[1])+1));
+    // 保存文件前缀 
+    char *s_prefix = (char*)malloc(sizeof(char)*(strlen(argv[1])+1));
+    if(2 != sscanf(argv[1], "%[^/]/%s", s_synpath, s_prefix)){
+        fprintf(stderr, "[%s] " BOLD_RED "Error format in \"%s\".\n" DEFAULT_RESTORE, command, argv[1]);
+        exit(EXIT_FAILURE);
+    }
+
+    // 检查是否存在该目录
+    DIR *dir = opendir(s_synpath);
+    if (dir == NULL) {
+        fprintf(stderr, "[%s] " BOLD_RED "Error! Directory \"%s\" not exists.\n" DEFAULT_RESTORE, command, s_synpath);
+        exit(EXIT_FAILURE);
+    } 
+
+
+
+    // ----------------------------------------------------------------------------------
+    // 开始读取计算，输出6个量
+    float *arrin = NULL;
+    char c1, c2;
+    char *s_filepath = (char*)malloc(sizeof(char) * (strlen(s_synpath)+strlen(s_prefix)+100));
+    const char chs[3] = {'R', 'T', 'Z'};
+
+    // 读取一个头段变量，获得基本参数，分配数组内存
+    SACHEAD hd;
+    sprintf(s_filepath, "%s/r%sR.sac", s_synpath, s_prefix);
+    read_SAC_HEAD(command, s_filepath, &hd);
+    int npts=hd.npts;
+    float dist=hd.dist;
+    float *arrout = (float*)calloc(npts, sizeof(float));
+
+    // ----------------------------------------------------------------------------------
+    // 循环6个分量
+    for(int i1=0; i1<3; ++i1){
+        c1 = chs[i1];
+        for(int i2=i1; i2<3; ++i2){
+            c2 = chs[i2];
+
+            // 读取数据 u_{i,j}
+            sprintf(s_filepath, "%s/%c%s%c.sac", s_synpath, tolower(c1), s_prefix, c2);
+            arrin = read_SAC(command, s_filepath, &hd, arrin);
+
+            // 累加
+            for(int i=0; i<npts; ++i)  arrout[i] += arrin[i];
+
+            // 读取数据 u_{j,i}
+            sprintf(s_filepath, "%s/%c%s%c.sac", s_synpath, tolower(c2), s_prefix, c1);
+            arrin = read_SAC(command, s_filepath, &hd, arrin);
+
+            // 累加
+            for(int i=0; i<npts; ++i)  arrout[i] = (arrout[i] + arrin[i]) * 0.5f;
+
+            // 特殊情况需加上协变导数
+            if(c1=='R' && c2=='T'){
+                // 读取数据 u_T
+                sprintf(s_filepath, "%s/%sT.sac", s_synpath, s_prefix);
+                arrin = read_SAC(command, s_filepath, &hd, arrin);
+                for(int i=0; i<npts; ++i)  arrout[i] -= 0.5f * arrin[i] / dist;
+            }
+            else if(c1=='T' && c2=='T'){
+                // 读取数据 u_R
+                sprintf(s_filepath, "%s/%sR.sac", s_synpath, s_prefix);
+                arrin = read_SAC(command, s_filepath, &hd, arrin);
+                for(int i=0; i<npts; ++i)  arrout[i] += arrin[i] / dist;
+            }
+
+            // 保存到SAC
+            sprintf(hd.kcmpnm, "%c%c", c1, c2);
+            sprintf(s_filepath, "%s/%s.strain.%c%c.sac", s_synpath, s_prefix, c1, c2);
+            write_sac(s_filepath, hd, arrout);
+
+            // 置零
+            for(int i=0; i<npts; ++i)  arrout[i] = 0.0f;
+        }
+    }
+
+    if(arrin)   free(arrin);
+    if(arrout)  free(arrout);
+    free(s_filepath);
+    free(s_synpath);
+    free(s_prefix);
+}