Additional comments to make the code clearer.

This commit is contained in:
Ludger Küperkoch 2015-12-14 09:34:56 +01:00
parent f6930618f2
commit 5f0e59d95a

View File

@ -1,7 +1,7 @@
#!/usr/bin/env python
# -*- coding: utf-8 -*-
"""
Created August/September 2015.
Created autumn/winter 2015.
:author: Ludger Küperkoch / MAGS2 EP3 working group
"""
@ -18,7 +18,7 @@ from pylot.core.read.data import Data
class Magnitude(object):
'''
Superclass for calculating Wood-Anderson peak-to-peak
amplitudes, local magnitudes, seismic moments
amplitudes, local magnitudes, source spectra, seismic moments
and moment magnitudes.
'''
@ -52,7 +52,7 @@ class Magnitude(object):
:param: vp [m/s], P-velocity
:param: integer
:param: invdir, path to inventory or dataless-SEED file
:param: invdir, name and path to inventory or dataless-SEED file
:type: string
'''
@ -265,6 +265,21 @@ def calcMoMw(wfstream, w0, rho, vp, delta, inv):
'''
Subfunction of run_calcMoMw to calculate individual
seismic moments and corresponding moment magnitudes.
:param: wfstream
:type: `~obspy.core.stream.Stream`
:param: w0, height of plateau of source spectrum
:type: float
:param: rho, rock density [kg/]
:type: integer
:param: delta, hypocentral distance [km]
:type: integer
:param: inv, name/path of inventory or dataless-SEED file
:type: string
'''
tr = wfstream[0]
@ -293,8 +308,35 @@ def calcsourcespec(wfstream, onset, inventory, vp, delta, azimuth, incidence, Qp
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
source model. Has to be derived from instrument corrected displacement traces,
thus restitution and integration necessary! Integrated traces have to be rotated
into ray-coordinate system ZNE => LQT!
thus restitution and integration necessary! Integrated traces are rotated
into ray-coordinate system ZNE => LQT using Obspy's rotate modul!
:param: wfstream
:type: `~obspy.core.stream.Stream`
:param: onset, P-phase onset time
:type: float
:param: inventory, path/name of inventory or dataless-SEED file
:type: string
:param: vp, Vp-wave velocity
:type: float
:param: delta, hypocentral distance [km]
:type: integer
:param: azimuth
:type: integer
:param: incidence
:type: integer
:param: Qp, quality factor for P-waves
:type: integer
:param: iplot, show results (iplot>1) or not (iplot(<2)
:type: integer
'''
print ("Calculating source spectrum ....")
@ -346,7 +388,11 @@ def calcsourcespec(wfstream, onset, inventory, vp, delta, azimuth, incidence, Qp
LQT=cordat_copy.rotate('ZNE->LQT',azimuth, incidence)
ldat = LQT.select(component="L")
if len(ldat) == 0:
# yet Obspy's rotate can not handle channels 3/2/1
# if horizontal channels are 2 and 3
# no azimuth information is available and thus no
# rotation is possible!
print("calcsourcespec: Azimuth information is missing, "
"no rotation of components possible!")
ldat = LQT.select(component="Z")
# integrate to displacement
@ -410,10 +456,12 @@ def calcsourcespec(wfstream, onset, inventory, vp, delta, azimuth, incidence, Qp
fi = np.where((f >= 1) & (f < 100))
F = f[fi]
YY = Y[fi]
# correction for attenuation
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
# initial guess of plateau
w0in = np.mean(YYcor[0:100])
@ -433,7 +481,8 @@ def calcsourcespec(wfstream, onset, inventory, vp, delta, azimuth, incidence, Qp
# use of conventional fitting
[w02, fc2] = fitSourceModel(F, YYcor, Fcin, iplot)
# get w0 and fc as median
# get w0 and fc as median of both
# source spectrum fits
w0 = np.median([w01, w02])
fc = np.median([fc1, fc2])
print("calcsourcespec: Using w0-value = %e m/Hz and fc = %f Hz" % (w0, fc))