-
Notifications
You must be signed in to change notification settings - Fork 1
/
syn_seismogram.py
33 lines (27 loc) · 909 Bytes
/
syn_seismogram.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
from database import AxiBasicDB
import numpy as np
import os
basedir = '../SOLVER/ak135_th/'
# create dir
os.makedirs("SEISMOGRAMS",exist_ok=True)
db:AxiBasicDB = AxiBasicDB()
db.read_basic(basedir + "/MZZ/Data/axisem_output.nc4")
db.set_iodata(basedir)
cmtfile = "CMTSOLUTION"
stacords = np.loadtxt(basedir + "/MZZ/STATIONS",dtype = str)
t = np.arange(db.nt) * db.dtsamp + db.t0
nsta = stacords.shape[0]
for i in range(nsta):
print(i+1,nsta)
stla,stlo = np.float32(stacords[i,2:4])
ue,un,uz = db.syn_seismo(stla,stlo,0.,'enz',basedir + cmtfile)
name = stacords[i,1] + "." + stacords[i,0]
newname = "SEISMOGRAMS/" + name
dataout = np.zeros((len(un),2))
dataout[:,0] = t
dataout[:,1] = ue
np.savetxt(f'{newname}.BXE.dat',dataout)
dataout[:,1] = un
np.savetxt(f'{newname}.BXN.dat',dataout)
dataout[:,1] = uz
np.savetxt(f'{newname}.BXZ.dat',dataout)