From 30ee81a39d1d5238fb43dfa982ccf7ac3bb2578b Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Ludger=20K=C3=BCperkoch?= <kueperkoch@bestec-for-nature.com>
Date: Wed, 23 Sep 2015 16:31:48 +0200
Subject: [PATCH] New class DCfc of object Magnitude for calculating source
 spectrum and to derive DC value and corner frequency.

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

diff --git a/pylot/core/analysis/magnitude.py b/pylot/core/analysis/magnitude.py
index 072984eb..2af92b37 100644
--- a/pylot/core/analysis/magnitude.py
+++ b/pylot/core/analysis/magnitude.py
@@ -25,7 +25,8 @@ class Magnitude(object):
         :type:  float
 
         :param: pwin, pick window [To To+pwin] to get maximum
-                peak-to-peak amplitude
+                peak-to-peak amplitude (WApp) or to calculate
+                source spectrum (DCfc)
         :type:  float
 
         :param: iplot, no. of figure window for plotting interims results
@@ -40,6 +41,7 @@ class Magnitude(object):
         self.setpwin(pwin)
         self.setiplot(iplot)
         self.calcwapp()
+        self.calcsourcespec()
 
 
     def getwfstream(self):
@@ -68,18 +70,23 @@ class Magnitude(object):
 
     def getwapp(self):
         return self.wapp
+  
+    def getw0(self):
+        return self.w0
+
+    def getfc(self):
+        return self.fc
 
     def calcwapp(self):
         self.wapp = None
 
-
     def calcsourcespec(self):
         self.sourcespek = None
 
 class WApp(Magnitude):
     '''
     Method to derive peak-to-peak amplitude as seen on a Wood-Anderson-
-    seismograph. Has to be derived from corrected traces!
+    seismograph. Has to be derived from instrument corrected traces!
     '''
 
     def calcwapp(self):
@@ -110,6 +117,7 @@ class WApp(Magnitude):
         iwin = getsignalwin(th, self.getTo(), self.getpwin())
         self.wapp = np.max(sqH[iwin])
         print ("Determined Wood-Anderson peak-to-peak amplitude: %f mm") % self.wapp
+
         if self.getiplot() > 1:
             stream.plot()
             f = plt.figure(2)
@@ -128,10 +136,48 @@ class WApp(Magnitude):
 class DCfc(Magnitude):
     '''
     Method to calculate the source spectrum and to derive from that the plateau 
-    (the so-called DC-value) and the corner frequency assuming Aki's omega-square 
-    source model. Has to be derived from corrected traces!
+    (so-called DC-value) and the corner frequency assuming Aki's omega-square 
+    source model. Has to be derived from instrument corrected displacement traces!
     '''
 
     def calcsourcespec(self):
     	print ("Calculating source spectrum ....")
 
+        self.w0 = None # DC-value
+        self.fc = None # corner frequency 
+
+        stream = self.getwfstream()
+        tr = stream[0]
+
+        # get time array
+        t = np.arange(0, len(tr) * tr.stats.delta, tr.stats.delta)
+        iwin = getsignalwin(t, self.getTo(), self.getpwin())
+        xdat = tr.data[iwin]
+
+        # fft 
+        fny = tr.stats.sampling_rate / 2
+        N = 1024
+        y = tr.stats.delta * np.fft.fft(xdat, N)
+        Y = abs(y[: N/2])
+        L = (N - 1) / tr.stats.sampling_rate
+        f = np.arange(0, fny, 1 / L)
+        
+        #if self.getiplot() > 1:
+        iplot=2
+        if iplot > 1:
+            f1 = plt.figure(1)
+            plt.subplot(2,1,1)
+            plt.plot(t, np.multiply(tr, 1000), 'k') # show displacement in mm
+            plt.plot(t[iwin], np.multiply(xdat, 1000), 'g') # show displacement in mm
+            plt.title('Seismogram and P pulse, station %s' % tr.stats.station)
+            plt.xlabel('Time since %s' % tr.stats.starttime)
+            plt.ylabel('Displacement [mm]')
+
+            plt.subplot(2,1,2)
+            plt.semilogy(f, Y.real)
+            plt.title('Source Spectrum from P Pulse')
+            plt.xlabel('Frequency [Hz]')
+            plt.ylabel('Amplitude [m/Hz]')
+            plt.show()
+            raw_input()
+            plt.close(f1)