From 8035903fa5025955c822595f2c67423eaa93a77a Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Ludger=20K=C3=BCperkoch?= <kueperkoch@bestec-for-nature.com>
Date: Tue, 29 Sep 2015 11:15:51 +0200
Subject: [PATCH] Zero crossings are calculated to derive only P pulse for
 calculating source spectrum.

---
 pylot/core/pick/autopick.py | 43 +++++++++++++++++++++++--------------
 1 file changed, 27 insertions(+), 16 deletions(-)

diff --git a/pylot/core/pick/autopick.py b/pylot/core/pick/autopick.py
index 10b4f7d5..e3771d7a 100755
--- a/pylot/core/pick/autopick.py
+++ b/pylot/core/pick/autopick.py
@@ -313,7 +313,6 @@ def autopickstation(wfstream, pickparam):
                 ##############################################################
                 # get DC value (w0) and corner frequency (fc) of source spectrum
                 # from P pulse
-                # restitute streams
                 # initialize Data object
                 data = Data()
                 [corzdat, restflag] = data.restituteWFData(invdir, zdat)
@@ -323,33 +322,45 @@ def autopickstation(wfstream, pickparam):
                         # class needs stream object => build it
                         z_copy = zdat.copy()
                         z_copy[0].data = corintzdat
-                	# calculate source spectrum and get w0 and fc
-                        calcwin = 1 / bpz2[0] # largest detectable period == window length 
-                                              # around P pulse for calculating source spectrum
-                	specpara = DCfc(z_copy, mpickP, calcwin, iplot)
-                        w0 = specpara.getw0()
-                        fc = specpara.getfc()
+                        # largest detectable period == window length
+                        # after P pulse for calculating source spectrum
+                        winzc = (1 / bpz2[0]) * z_copy[0].stats.sampling_rate  
+                        impickP = mpickP * z_copy[0].stats.sampling_rate
+                        wfzc = z_copy[0].data[impickP : impickP + winzc]
+                        # calculate spectrum using only first cycles of
+                        # waveform after P onset!
+                        zc = crossings_nonzero_all(wfzc)
+                        if np.size(zc) == 0:
+                            print ("Something is wrong with the waveform, " \
+                                   "no zero crossings derived!")
+                            print ("Cannot calculate source spectrum!")
+                        else:
+                            calcwin = (zc[3] - zc[0]) * z_copy[0].stats.delta
+                	    # calculate source spectrum and get w0 and fc
+                	    specpara = DCfc(z_copy, mpickP, calcwin, iplot)
+                            w0 = specpara.getw0()
+                            fc = specpara.getfc()
 
-                print 'autopickstation: P-weight: %d, SNR: %f, SNR[dB]: %f, ' \
-                      'Polarity: %s' % (Pweight, SNRP, SNRPdB, FM)
+                print ("autopickstation: P-weight: %d, SNR: %f, SNR[dB]: %f, " \
+                       "Polarity: %s" % (Pweight, SNRP, SNRPdB, FM))
                 Sflag = 1
 
         else:
-            print 'Bad initial (AIC) P-pick, skipping this onset!'
+            print ("Bad initial (AIC) P-pick, skipping this onset!")
             print 'AIC-SNR=', aicpick.getSNR(), 'AIC-Slope=', aicpick.getSlope(), 'counts/s'
             print '(min. AIC-SNR=', minAICPSNR, ', min. AIC-Slope=', minAICPslope, 'counts/s)'
             Sflag = 0
 
     else:
-        print 'autopickstation: No vertical component data available!, ' \
-              'Skipping station!'
+        print ("autopickstation: No vertical component data available!, " \
+               "Skipping station!")
 
     if edat is not None and ndat is not None and len(edat) > 0 and len(
             ndat) > 0 and Pweight < 4:
-        print 'Go on picking S onset ...'
-        print '##################################################'
-        print 'Working on S onset of station %s' % edat[0].stats.station
-        print 'Filtering horizontal traces ...'
+        print ("Go on picking S onset ...")
+        print ("##################################################")
+        print ("Working on S onset of station %s" % edat[0].stats.station)
+        print ("Filtering horizontal traces ...")
 
         # determine time window for calculating CF after P onset
         cuttimesh = [round(max([mpickP + sstart, 0])),