TEST: add compare_results_liquid example

example/compare_results_liquid/liquid

@@ -0,0 +1,5 @@
+4.0     6.0        3.0      2.5       9e10  9e10
+4.0     6.2        0.0      2.6        9e10  9e10
+7.0     6.4        3.5      2.7        9e10  9e10
+26.0     6.6        3.6      2.8       9e10  9e10
+0.0     8.0        4.7      3.3       9e10  9e10

example/compare_results_liquid/liquid.mod

@@ -0,0 +1,17 @@
+MODEL.01
+LIQUID
+ISOTROPIC
+KGS
+FLAT EARTH
+1-D
+CONSTANT VELOCITY
+LINE08
+LINE09
+LINE10
+LINE11
+      H(KM)   VP(KM/S)   VS(KM/S) RHO(GM/CC)         QP         QS       ETAP       ETAS      FREFP      FREFS
+      4.0     6.0        3.0      2.5         0 0 0 0 1 1  
+      4.0     6.2        0.0      2.6         0 0 0 0 1 1  
+      7.0     6.4        3.5      2.7         0 0 0 0 1 1  
+     26.0     6.6        3.6      2.8         0 0 0 0 1 1  
+      0.0     8.0        4.7      3.3         0 0 0 0 1 1  

example/compare_results_liquid/plot_cps_grt.py

@@ -0,0 +1,74 @@
+import numpy as np 
+import matplotlib.pyplot as plt 
+from obspy import * 
+
+def plot(st_grt:Stream, st_cps:Stream, out:str):
+    nt = st_grt[0].stats.npts
+    dt = st_grt[0].stats.delta
+    t = np.arange(0, nt)*dt
+
+    fig, axs = plt.subplots(6, 3, figsize=(12, 15), gridspec_kw=dict(hspace=0.0, wspace=0.1)) # 
+    axs = axs.ravel()
+
+    itr = 0
+    for src, src2 in zip(['EX', 'VF', 'HF', 'DD', 'DS', 'SS'],
+                        ['Explosion', 'Vertical\nForce', 'Horizontal\nForce', '45°\nDip Slip', '90°\nDip Slip', 'Strike Slip']):
+        for comp, comp2 in zip(['Z', 'R', 'T'], ['Vertical', 'Radial', 'Transverse']):
+            ax = axs[itr]
+
+            for (_,spine) in ax.spines.items():
+                spine.set_visible(False)
+            ax.set_xticks([])
+            ax.set_yticks([])
+
+            if itr < 3:
+                ax.set_title(comp2, fontsize=18, pad=10)
+
+            itr += 1
+            try:
+                tr = st_grt.select(channel=f'{src}{comp}')[0]
+                tr_cps = st_cps.select(channel=f'{comp}{src}')[0]
+                data1 = tr.data.copy()
+                data2 = tr_cps.data.copy()
+            except:
+                # data1 = data2 = np.zeros_like(t)
+                continue
+
+
+            norm = np.linalg.norm(data1, np.inf)+1e-10
+
+            # adjust sign
+            if comp == 'T':
+                data2 *= -1
+            if src == 'DS':
+                data2 *= -1
+
+
+            ax.plot(t, data2/norm, c='k', ls='-', lw=1.8)
+            ax.plot(t, data1/norm, c='r', ls='-', lw=0.4)
+
+            ax.set_xlim([10, 60])
+            ax.set_ylim([-1.1,1.1])
+
+            ylen = ax.get_ylim()[1] - ax.get_ylim()[0]
+
+            if comp == 'Z':
+                ax.text(0.8, 0, src2, fontsize=18, clip_on=False, ha='center', va='center', )
+
+
+    for i in range(1, 4):        
+        ax = axs[-i]
+        ax.spines['bottom'].set_visible(True)
+        ax.set_xticks(np.arange(10, 70, 10))
+        ax.set_xlabel('Time (s)', fontsize=12)
+
+    fig.savefig(out)
+
+
+
+
+if __name__=='__main__':
+    st_grt = read("GRN/liquid_20_8.2_100/*")
+    st_cps = read("liquid_sdep20_rdep8.2/GRN/*")
+
+    plot(st_grt, st_cps, "compare_cps_grt.pdf")

example/compare_results_liquid/plot_cps_pygrt.py

@@ -0,0 +1,36 @@
+import pygrt
+import numpy as np 
+import matplotlib.pyplot as plt 
+from obspy import * 
+
+
+# load model 
+modarr = np.loadtxt("./liquid")
+
+depsrc = 20.0   
+deprcv = 8.2    
+pymod = pygrt.PyModel1D(modarr, depsrc, deprcv) 
+
+rs = np.array([100]) 
+
+nt = 512   
+dt = 0.2   
+zeta = 0.8  
+
+# compute Green's Functions
+st_grt = pymod.compute_grn(
+    distarr=rs, 
+    nt=nt, 
+    dt=dt, 
+    zeta=zeta, 
+    Length=20,
+)[0]
+
+
+try:
+    st_cps = read("liquid_sdep20_rdep8.2/GRN/*")
+
+    from plot_cps_grt import plot
+    plot(st_grt, st_cps, "compare_cps_pygrt.pdf")
+except Exception as e:
+    print(str(e))

example/compare_results_liquid/run_cps330.sh

@@ -0,0 +1,36 @@
+#!/bin/bash
+
+dist=100
+depsrc=20
+deprcv=8.2
+
+nt=512
+dt=0.2
+
+XF=1.0
+XL=$(echo $dist*20 | bc)
+
+modname="liquid"
+mod="${modname}.mod"
+
+maindir=${modname}_sdep${depsrc}_rdep${deprcv}
+mkdir -p ${maindir}
+cd ${maindir}
+# mkdir -p ${modname}
+# cd ${modname}
+
+cat > dfile << EOF
+$dist $dt $nt 0 0
+EOF
+
+
+# compute Green's Functions
+hprep96 -M ../$mod -d dfile -HS $depsrc -HR $deprcv -ALL -XL $XL -XF $XF
+rspec96 
+
+# to ascii
+hpulse96 -D -i -IMP > hpulse96.out
+
+mkdir -p GRN
+f96tosac -B < hpulse96.out
+mv *.sac GRN/

example/compare_results_liquid/run_grt.sh

@@ -0,0 +1,14 @@
+#!/bin/bash
+
+dist=100
+depsrc=20
+deprcv=8.2
+
+nt=512
+dt=0.2
+
+modname="liquid"
+out="GRN"
+
+grt -M${modname} -O${out} -N${nt}/${dt} -D${depsrc}/${deprcv} -R${dist} -L20
+