marcel/pylot
3
0
Fork 0
Code Issues 7 Pull Requests 1 Releases Wiki Activity

calcMoMw: modified for calculating Mo and Mw using metric units. calcsourcespec: modified correction for attenutation using an exponential including Q(f).

Browse Source
This commit is contained in:
Ludger Küperkoch 2015-12-11 15:53:35 +01:00
parent f29d285910
commit f6930618f2
1 changed files with 10 additions and 5 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

15
pylot/core/analysis/magnitude.py
View File

@ -240,8 +240,8 @@ class M0Mw(Magnitude):
# call subfunction to estimate source spectrum # call subfunction to estimate source spectrum
# and to derive w0 and fc # and to derive w0 and fc
[w0, fc] = calcsourcespec(selwf, picks[key]['P']['mpp'], \ [w0, fc] = calcsourcespec(selwf, picks[key]['P']['mpp'], \
self.getinvdir(), az, inc, self.getQp(), \ self.getinvdir(), self.getvp(), delta, az, \
self.getiplot()) inc, self.getQp(), self.getiplot())
if w0 is not None: if w0 is not None:
# call subfunction to calculate Mo and Mw # call subfunction to calculate Mo and Mw
@ -268,6 +268,7 @@ def calcMoMw(wfstream, w0, rho, vp, delta, inv):
''' '''
tr = wfstream[0] tr = wfstream[0]
delta = delta * 1000 # hypocentral distance in [m]
print("calcMoMw: Calculating seismic moment Mo and moment magnitude Mw for station %s ..." \ print("calcMoMw: Calculating seismic moment Mo and moment magnitude Mw for station %s ..." \
% tr.stats.station) % tr.stats.station)
@ -278,7 +279,8 @@ def calcMoMw(wfstream, w0, rho, vp, delta, inv):
Mo = w0 * 4 * np.pi * rho * np.power(vp, 3) * delta / (rP * freesurf) Mo = w0 * 4 * np.pi * rho * np.power(vp, 3) * delta / (rP * freesurf)
Mw = np.log10(Mo * 1e07) * 2 / 3 - 10.7 #after Hanks & Kanamori (1979), defined for [dyn*cm]! #Mw = np.log10(Mo * 1e07) * 2 / 3 - 10.7 # after Hanks & Kanamori (1979), defined for [dyn*cm]!
Mw = np.log10(Mo) * 2 / 3 - 6.7 # for metric units
print("calcMoMw: Calculated seismic moment Mo = %e Nm => Mw = %3.1f " % (Mo, Mw)) print("calcMoMw: Calculated seismic moment Mo = %e Nm => Mw = %3.1f " % (Mo, Mw))
@ -286,7 +288,7 @@ def calcMoMw(wfstream, w0, rho, vp, delta, inv):
def calcsourcespec(wfstream, onset, inventory, azimuth, incidence, Qp, iplot): def calcsourcespec(wfstream, onset, inventory, vp, delta, azimuth, incidence, Qp, iplot):
''' '''
Subfunction to calculate the source spectrum and to derive from that the plateau Subfunction to calculate the source spectrum and to derive from that the plateau
(usually called omega0) and the corner frequency assuming Aki's omega-square (usually called omega0) and the corner frequency assuming Aki's omega-square
@ -302,6 +304,7 @@ def calcsourcespec(wfstream, onset, inventory, azimuth, incidence, Qp, iplot):
Q = int(qu[0]) Q = int(qu[0])
# A, i.e. power of frequency # A, i.e. power of frequency
A = float(qu[1]) A = float(qu[1])
delta = delta * 1000 # hypocentral distance in [m]
fc = None fc = None
w0 = None w0 = None
@ -408,7 +411,9 @@ def calcsourcespec(wfstream, onset, inventory, azimuth, incidence, Qp, iplot):
F = f[fi] F = f[fi]
YY = Y[fi] YY = Y[fi]
# correction for attenuation # correction for attenuation
YYcor = YY.real*Q**A wa = 2 * np.pi * F #angular frequency
D = np.exp((wa * delta) / (2 * vp * Q*F**A))
YYcor = YY.real*D
# get plateau (DC value) and corner frequency # get plateau (DC value) and corner frequency
# initial guess of plateau # initial guess of plateau
w0in = np.mean(YYcor[0:100]) w0in = np.mean(YYcor[0:100])