Add some features and fix some bugs about source signals

.github/workflows/testbuild.yml

@@ -274,6 +274,13 @@ jobs:
       #     python plot_compare_pygrt.py
       #   continue-on-error: true  # 即使失败，仍然标记为成功
 
+      - name: Test signal convolution
+        working-directory: ${{ env.EXAMPLE_COPY_PATH }}/convolve_signals
+        run: |
+          chmod +x *.sh
+          ./run_grt.sh
+          python run_pygrt.py
+
       - name: Remove the test files
         run: |
           rm -rf ${{ env.EXAMPLE_COPY_PATH }}

example/convolve_signals/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
+28.  6.9  4.0  2.8  9e10  9e10
+-0.1   8.2  4.7  3.2  9e10  9e10

example/convolve_signals/plot.sh

@@ -0,0 +1,27 @@
+#!/bin/bash
+
+SAC_DISPLAY_COPYRIGHT=0
+
+out="compare_pdf"
+rm -rf $out
+mkdir -p $out 
+
+
+for synpath in syn*; do 
+for cout in $synpath/cout*; do
+fname=$(basename $cout)
+echo $cout $(dirname $cout)/p${fname:1:}
+sac << EOF
+r $cout $(dirname $cout)/p${fname:1}
+qdp off
+color red inc list blue red 
+p1
+saveimg $out/compare_${synpath}_${fname}.pdf
+q
+EOF
+
+done
+done
+
+
+pdftk $out/* cat output out.pdf

example/convolve_signals/run_grt.sh

@@ -0,0 +1,40 @@
+#!/bin/bash
+
+
+dist=10
+depsrc=2
+deprcv=0.5
+
+nt=1024
+dt=0.01
+
+modname="milrow"
+out="GRN"
+
+grt -M${modname} -O${out} -N${nt}/${dt} -D${depsrc}/${deprcv} -R${dist}
+
+# convolve different signals
+G=$out/${modname}_${depsrc}_${deprcv}_${dist}
+P=cout
+S=1e24
+az=39.2
+grt.syn -G$G -Osyn_exp -P$P -A$az -S$S -Dt/0.2/0.2/0.4
+
+fn=2
+fe=-1
+fz=4
+grt.syn -G$G -Osyn_sf -P$P -A$az -S$S -F$fn/$fe/$fz -Dt/0.1/0.3/0.6
+
+stk=77
+dip=88
+rak=99
+grt.syn -G$G -Osyn_dc -P$P -A$az -S$S -M$stk/$dip/$rak -Dp/0.6
+
+M11=1
+M12=-2
+M13=-5
+M22=0.5
+M23=3
+M33=1.2
+grt.syn -G$G -Osyn_mt -P$P -A$az -S$S -T$M11/$M12/$M13/$M22/$M23/$M33 -Dr/3
+

example/convolve_signals/run_pygrt.py

@@ -0,0 +1,41 @@
+import numpy as np
+import pygrt 
+
+dist=10
+depsrc=2
+deprcv=0.5
+
+nt=1024
+dt=0.01
+
+modname="milrow"
+
+modarr = np.loadtxt(modname)
+
+pymod = pygrt.PyModel1D(modarr, depsrc, deprcv)
+
+# compute green functions
+st_grn = pymod.gen_gf_spectra(dist, nt, dt)[0]
+
+# synthetic
+S=1e24
+az=39.2
+st = pygrt.utils.gen_syn_from_gf_EXP(st_grn, S, az)
+sigs = pygrt.sigs.gen_triangle_wave(0.4, dt)
+pygrt.utils.stream_convolve(st, sigs)
+pygrt.utils.stream_write_sac(st, "syn_exp/pout")
+
+st = pygrt.utils.gen_syn_from_gf_SF(st_grn, S, 2, -1, 4, az)
+sigs = pygrt.sigs.gen_trap_wave(0.1, 0.3, 0.6, dt)
+pygrt.utils.stream_convolve(st, sigs)
+pygrt.utils.stream_write_sac(st, "syn_sf/pout")
+
+st = pygrt.utils.gen_syn_from_gf_DC(st_grn, S, 77, 88, 99, az)
+sigs = pygrt.sigs.gen_parabola_wave(0.6, dt)
+pygrt.utils.stream_convolve(st, sigs)
+pygrt.utils.stream_write_sac(st, "syn_dc/pout")
+
+st = pygrt.utils.gen_syn_from_gf_MT(st_grn, S, [1,-2,-5,0.5,3,1.2], az)
+sigs = pygrt.sigs.gen_ricker_wave(3, dt)
+pygrt.utils.stream_convolve(st, sigs)
+pygrt.utils.stream_write_sac(st, "syn_mt/pout")

pygrt/C_extension/include/dynamic/signals.h

@@ -167,4 +167,11 @@ float * get_ricker_wave(float dt, float f0, int *Nt);
  * 
  * @return   float指针
  */
-float * get_custom_wave(int *Nt, const char *tfparams);
+float * get_custom_wave(int *Nt, const char *tfparams);
+
+/**
+ * 专用于在Python端释放C中申请的内存
+ * 
+ * @param     pt    指针
+ */
+void free1d(void *pt);

pygrt/C_extension/src/dynamic/signals.c

@@ -348,3 +348,8 @@ float * get_custom_wave(int *Nt, const char *tfparams){
     *Nt = nt;
     return tfarr;
 }
+
+
+void free1d(void *pt){
+    free(pt);
+}

pygrt/c_interfaces.py

@@ -13,6 +13,9 @@
 
 from .c_structures import * 
 
+FPOINTER = POINTER(c_float)
+IPOINTER = POINTER(c_int)
+
 
 c_PGRN = POINTER(c_GRN)
 
@@ -81,4 +84,29 @@ def set_num_threads(n):
 C_compute_travt1d.argtypes = [
     PREAL, PREAL, c_int, 
     c_int, c_int, REAL
-]
+]
+
+
+
+# -------------------------------------------------------------------
+#                      C函数定义的时间函数
+# -------------------------------------------------------------------
+C_free = libgrt.free1d
+"""释放在C中申请的内存"""
+C_free.restype = None
+C_free.argtypes = [c_void_p]
+
+C_get_trap_wave = libgrt.get_trap_wave
+"""梯形波"""
+C_get_trap_wave.restype = FPOINTER
+C_get_trap_wave.argtypes = [c_float, FPOINTER, FPOINTER, FPOINTER, IPOINTER]
+
+C_get_parabola_wave = libgrt.get_parabola_wave
+"""抛物波"""
+C_get_parabola_wave.restype = FPOINTER
+C_get_parabola_wave.argtypes = [c_float, FPOINTER, IPOINTER]
+
+C_get_ricker_wave = libgrt.get_ricker_wave
+"""雷克子波"""
+C_get_ricker_wave.restype = FPOINTER
+C_get_ricker_wave.argtypes = [c_float, c_float, IPOINTER]

pygrt/pygrn.py

@@ -153,7 +153,7 @@ def freq2time(self, T0:float, travtP:float, travtS:float, mult:float=1.0):
         sac = self.SACTrace
         nt = sac.npts
         dt = sac.delta
-        wI = self.wI
+        wI = sac.user0
 
         T = nt*dt
         if not np.isclose(T*df, 1.0):