diff --git a/example/compare_results_liquid/compare_cps_grt.pdf b/example/compare_results_liquid/compare_cps_grt.pdf new file mode 100644 index 00000000..db28951e Binary files /dev/null and b/example/compare_results_liquid/compare_cps_grt.pdf differ diff --git a/example/compare_results_liquid/compare_cps_pygrt.pdf b/example/compare_results_liquid/compare_cps_pygrt.pdf new file mode 100644 index 00000000..ba9f0fba Binary files /dev/null and b/example/compare_results_liquid/compare_cps_pygrt.pdf differ diff --git a/example/compare_results_liquid/liquid b/example/compare_results_liquid/liquid new file mode 100644 index 00000000..e15f7611 --- /dev/null +++ b/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 \ No newline at end of file diff --git a/example/compare_results_liquid/liquid.mod b/example/compare_results_liquid/liquid.mod new file mode 100644 index 00000000..4b2d3443 --- /dev/null +++ b/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 \ No newline at end of file diff --git a/example/compare_results_liquid/plot_cps_grt.py b/example/compare_results_liquid/plot_cps_grt.py new file mode 100644 index 00000000..3bdf9593 --- /dev/null +++ b/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") \ No newline at end of file diff --git a/example/compare_results_liquid/plot_cps_pygrt.py b/example/compare_results_liquid/plot_cps_pygrt.py new file mode 100644 index 00000000..d2644613 --- /dev/null +++ b/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)) diff --git a/example/compare_results_liquid/run_cps330.sh b/example/compare_results_liquid/run_cps330.sh new file mode 100644 index 00000000..62b2e01f --- /dev/null +++ b/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/ diff --git a/example/compare_results_liquid/run_grt.sh b/example/compare_results_liquid/run_grt.sh new file mode 100644 index 00000000..b4e7bac1 --- /dev/null +++ b/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 +