From 5f0e59d95a1a590dc8719487b069ee7b6082c8d8 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Ludger=20K=C3=BCperkoch?= <kueperkoch@bestec-for-nature.com>
Date: Mon, 14 Dec 2015 09:34:56 +0100
Subject: [PATCH] Additional comments to make the code clearer.

---
 pylot/core/analysis/magnitude.py | 63 ++++++++++++++++++++++++++++----
 1 file changed, 56 insertions(+), 7 deletions(-)

diff --git a/pylot/core/analysis/magnitude.py b/pylot/core/analysis/magnitude.py
index 4353b1e5..827d6055 100644
--- a/pylot/core/analysis/magnitude.py
+++ b/pylot/core/analysis/magnitude.py
@@ -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/m³]
+    :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))