From 4f76560752b3102c7c7281a6104a5f72fabc21d2 Mon Sep 17 00:00:00 2001 From: Dengda98 Date: Fri, 28 Mar 2025 15:02:39 +0800 Subject: [PATCH 1/8] FEAT: write Vp, Vs and Rho of receiver site to SAC header --- pygrt/C_extension/src/dynamic/grt_main.c | 9 ++++++--- 1 file changed, 6 insertions(+), 3 deletions(-) diff --git a/pygrt/C_extension/src/dynamic/grt_main.c b/pygrt/C_extension/src/dynamic/grt_main.c index 6e128ce2..3d40da89 100644 --- a/pygrt/C_extension/src/dynamic/grt_main.c +++ 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,12 @@ 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]; // 做反傅里叶变换，保存SAC文件 From 7f9e7ccefbe77e40bb3fde1bba964d3f5394da31 Mon Sep 17 00:00:00 2001 From: Dengda98 Date: Fri, 28 Mar 2025 15:17:57 +0800 Subject: [PATCH 2/8] FEAT: change -O in `get_syn.c` to be mandatory --- pygrt/C_extension/src/dynamic/grt_syn.c | 24 +++++++++++------------- 1 file changed, 11 insertions(+), 13 deletions(-) diff --git a/pygrt/C_extension/src/dynamic/grt_syn.c b/pygrt/C_extension/src/dynamic/grt_syn.c index 346c9768..85b21b8f 100644 --- a/pygrt/C_extension/src/dynamic/grt_syn.c +++ b/pygrt/C_extension/src/dynamic/grt_syn.c @@ -128,10 +128,10 @@ printf("\n" "\n\n" "Usage:\n" "----------------------------------------------------------------\n" -" grt.syn -G -A -S \n" +" grt.syn -G -A -S -O \n" " [-M//]\n" " [-T/////]\n" -" [-F//] [-O] \n" +" [-F//] \n" " [-D/] [-I] [-J]\n" " [-P] [-e] [-s]\n" "\n" @@ -452,6 +452,10 @@ static void getopt_from_command(int argc, char **argv){ fprintf(stderr, "[%s] " BOLD_RED "Error! Need set -S. Use '-h' for help.\n" DEFAULT_RESTORE, command); exit(EXIT_FAILURE); } + if(O_flag == 0){ + fprintf(stderr, "[%s] " BOLD_RED "Error! Need set -O. Use '-h' for help.\n" DEFAULT_RESTORE, command); + exit(EXIT_FAILURE); + } // 只能使用一种震源 if(M_flag + F_flag + T_flag > 1){ @@ -504,18 +508,12 @@ static void getopt_from_command(int argc, char **argv){ } - if(O_flag == 1){ - // 建立保存目录 - if(mkdir(s_output_dir, 0777) != 0){ - if(errno != EEXIST){ - fprintf(stderr, "[%s] " BOLD_RED "Error! Unable to create folder %s. Error code: %d\n" DEFAULT_RESTORE, command, s_output_dir, errno); - exit(EXIT_FAILURE); - } + // 建立保存目录 + if(mkdir(s_output_dir, 0777) != 0){ + if(errno != EEXIST){ + fprintf(stderr, "[%s] " BOLD_RED "Error! Unable to create folder %s. Error code: %d\n" DEFAULT_RESTORE, command, s_output_dir, errno); + exit(EXIT_FAILURE); } - } else { - // 使用当前目录 - s_output_dir = (char*)malloc(sizeof(char)*100); - strcpy(s_output_dir, "."); } if(P_flag == 0){ From 21a3503f7af7d334e87d3a612d7b9941f1792238 Mon Sep 17 00:00:00 2001 From: Dengda98 Date: Mon, 31 Mar 2025 22:56:09 +0800 Subject: [PATCH 3/8] DOC: add comment about the sign of i/2Q --- pygrt/C_extension/include/common/attenuation.h | 1 + 1 file changed, 1 insertion(+) diff --git a/pygrt/C_extension/include/common/attenuation.h b/pygrt/C_extension/include/common/attenuation.h index 9a066947..9e3d270c 100755 --- a/pygrt/C_extension/include/common/attenuation.h +++ b/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$ From 5129d6f927043b1b03ed3d924d5f8941f61162d2 Mon Sep 17 00:00:00 2001 From: Dengda98 Date: Mon, 31 Mar 2025 22:58:19 +0800 Subject: [PATCH 4/8] FEAT: Add SAC file reading functions to encapsulate related functions --- pygrt/C_extension/include/common/sacio2.h | 34 +++++++++++++++++++ pygrt/C_extension/src/common/sacio2.c | 41 +++++++++++++++++++++++ 2 files changed, 75 insertions(+) create mode 100644 pygrt/C_extension/include/common/sacio2.h create mode 100644 pygrt/C_extension/src/common/sacio2.c diff --git a/pygrt/C_extension/include/common/sacio2.h b/pygrt/C_extension/include/common/sacio2.h new file mode 100644 index 00000000..a5789b8e --- /dev/null +++ b/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); \ No newline at end of file diff --git a/pygrt/C_extension/src/common/sacio2.c b/pygrt/C_extension/src/common/sacio2.c new file mode 100644 index 00000000..9b5e64ed --- /dev/null +++ b/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 +#include + +#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; inpts; ++i) arrout[i] = arrin[i]; + free(arrin); + arrin = arrout; + } + + return arrin; +} \ No newline at end of file From 1ea1c96d4e8550aac286b10bb3280cf79697f9e5 Mon Sep 17 00:00:00 2001 From: Dengda98 Date: Mon, 31 Mar 2025 22:59:27 +0800 Subject: [PATCH 5/8] FEAT: use new `read_SAC` function --- pygrt/C_extension/src/dynamic/grt_b2a.c | 8 ++------ pygrt/C_extension/src/dynamic/grt_syn.c | 9 +++------ 2 files changed, 5 insertions(+), 12 deletions(-) diff --git a/pygrt/C_extension/src/dynamic/grt_b2a.c b/pygrt/C_extension/src/dynamic/grt_b2a.c index 24d5d003..f7c6aa90 100644 --- a/pygrt/C_extension/src/dynamic/grt_b2a.c +++ b/pygrt/C_extension/src/dynamic/grt_b2a.c @@ -14,7 +14,7 @@ #include #include -#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; diff --git a/pygrt/C_extension/src/dynamic/grt_syn.c b/pygrt/C_extension/src/dynamic/grt_syn.c index 85b21b8f..399c13e5 100644 --- a/pygrt/C_extension/src/dynamic/grt_syn.c +++ b/pygrt/C_extension/src/dynamic/grt_syn.c @@ -21,7 +21,7 @@ #include #include "dynamic/signals.h" -#include "common/sacio.h" +#include "common/sacio2.h" #include "common/const.h" #include "common/logo.h" #include "common/colorstr.h" @@ -698,11 +698,8 @@ int main(int argc, char **argv){ if(coef == 0.0) continue; sprintf(buffer, "%s/%s%s%c.sac", s_grnpath, sacin_prefixes[ityp], srcName[k], ch); - if((arr = read_sac(buffer, &hd)) == NULL){ - fprintf(stderr, "[%s] " BOLD_RED "read %s failed.\n" DEFAULT_RESTORE, command, buffer); - exit(EXIT_FAILURE); - } - + arr = read_SAC(command, buffer, &hd, NULL); + nt = hd.npts; dt = hd.delta; dist = hd.dist; From 46180b4723c7a273cc63797f02b0ba63dedae235 Mon Sep 17 00:00:00 2001 From: Dengda98 Date: Mon, 31 Mar 2025 23:00:30 +0800 Subject: [PATCH 6/8] FEAT: write Qa and Qb of receiver site in SAC header --- pygrt/C_extension/src/dynamic/grt_main.c | 2 ++ 1 file changed, 2 insertions(+) diff --git a/pygrt/C_extension/src/dynamic/grt_main.c b/pygrt/C_extension/src/dynamic/grt_main.c index 3d40da89..806ced25 100644 --- a/pygrt/C_extension/src/dynamic/grt_main.c +++ b/pygrt/C_extension/src/dynamic/grt_main.c @@ -1146,6 +1146,8 @@ int main(int argc, char **argv) { 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文件 From 3b2fe3ee1f781eb47f29df0aa078b838493a6548 Mon Sep 17 00:00:00 2001 From: Dengda98 Date: Mon, 31 Mar 2025 23:01:31 +0800 Subject: [PATCH 7/8] FEAT: add `grt.stress` and `grt.strain` to calculate stress and strain --- pygrt/C_extension/src/dynamic/grt_strain.c | 148 +++++++++++++ pygrt/C_extension/src/dynamic/grt_stress.c | 241 +++++++++++++++++++++ 2 files changed, 389 insertions(+) create mode 100644 pygrt/C_extension/src/dynamic/grt_strain.c create mode 100644 pygrt/C_extension/src/dynamic/grt_stress.c diff --git a/pygrt/C_extension/src/dynamic/grt_strain.c b/pygrt/C_extension/src/dynamic/grt_strain.c new file mode 100644 index 00000000..ecd0bc26 --- /dev/null +++ b/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 +#include +#include +#include +#include +#include + +#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 /\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 +#include +#include +#include +#include +#include + +#include +#include + +#include "common/attenuation.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.stress]\n\n" +" Conbine spatial derivatives of displacements into stress.\n" +"\n\n" +"Usage:\n" +"----------------------------------------------------------------\n" +" grt.stress /\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个量 + 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 dt=hd.delta; + float dist=hd.dist; + float df=1.0/(npts*dt); + int nf=npts/2 + 1; + float va=hd.user1; + float vb=hd.user2; + float rho=hd.user3; + float Qainv=hd.user4; + float Qbinv=hd.user5; + if(va < 0.0 || vb < 0.0 || rho < 0.0){ + fprintf(stderr, "[%s] " BOLD_RED "Error! read necessary header value from \"%s\" error.\n" DEFAULT_RESTORE, command, s_filepath); + exit(EXIT_FAILURE); + } + // 申请内存 + // lamda * 体积应变 + fftwf_complex *lam_ukk = (fftwf_complex*)fftwf_malloc(sizeof(fftwf_complex)*nf); + // 不同频率的lambda和mu + fftwf_complex *lams = (fftwf_complex*)fftwf_malloc(sizeof(fftwf_complex)*nf); + fftwf_complex *mus = (fftwf_complex*)fftwf_malloc(sizeof(fftwf_complex)*nf); + // 分配FFTW数组 + fftwf_complex *carrin = (fftwf_complex*)fftwf_malloc(sizeof(fftwf_complex)*nf); + fftwf_complex *carrout = (fftwf_complex*)fftwf_malloc(sizeof(fftwf_complex)*nf); + float *arrin = (float*)malloc(sizeof(float)*npts); + float *arrout = (float*)malloc(sizeof(float)*npts); + fftwf_plan plan = fftwf_plan_dft_r2c_1d(npts, arrin, carrin, FFTW_ESTIMATE); + fftwf_plan plan_inv = fftwf_plan_dft_c2r_1d(npts, carrout, arrout, FFTW_ESTIMATE); + // 初始化 + for(int i=0; i Date: Mon, 31 Mar 2025 23:01:51 +0800 Subject: [PATCH 8/8] FEAT: update Makefile --- pygrt/C_extension/src/dynamic/Makefile | 24 +++++++++++++++++------- pygrt/C_extension/src/travt/Makefile | 4 ++-- 2 files changed, 19 insertions(+), 9 deletions(-) diff --git a/pygrt/C_extension/src/dynamic/Makefile b/pygrt/C_extension/src/dynamic/Makefile index 81022052..742679ba 100644 --- a/pygrt/C_extension/src/dynamic/Makefile +++ b/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) diff --git a/pygrt/C_extension/src/travt/Makefile b/pygrt/C_extension/src/travt/Makefile index 2bedc7f0..6ff88d85 100644 --- a/pygrt/C_extension/src/travt/Makefile +++ b/pygrt/C_extension/src/travt/Makefile @@ -35,8 +35,8 @@ $(BUILD_DIR)/%.o: %.c $(BIN_DIR): @mkdir -p $(BIN_DIR) -$(BIN_DIR)/grt.travt: grt_travt.c $(OBJS) - $(CC) -o $@ $^ $(COMMON_OBJS) $(CFLAGS) +$(BIN_DIR)/grt.travt: grt_travt.c $(OBJS) $(COMMON_OBJS) + $(CC) -o $@ $^ $(CFLAGS) # ----------------------- Dependency generation ----------------------- -include $(DEPS)