From a0fa9cd560c322c915eb440178dea1a7b610d3da Mon Sep 17 00:00:00 2001
From: Cees Bassa <cgbsat@gmail.com>
Date: Sat, 13 Apr 2019 16:59:32 +0200
Subject: [PATCH] Misc changes

---
 python/rfplot.py    | 12 +++++-------
 python/strf/cmap.py |  2 +-
 python/strf/rfio.py | 23 +++++++++++++++++++++--
 3 files changed, 27 insertions(+), 10 deletions(-)

diff --git a/python/rfplot.py b/python/rfplot.py
index 29a066f..1d40805 100755
--- a/python/rfplot.py
+++ b/python/rfplot.py
@@ -1,4 +1,4 @@
-#!/usr/bin/env python3
+#!/usr/bin/env python
 import numpy as np
 import ppgplot
 from strf.rfio import Spectrogram
@@ -10,17 +10,15 @@ if __name__ == "__main__":
     prefix = "2018-12-31T09:34:11"
     ifile = 48000
     nsub = 900
-    path = "/data2/satobs/radio/20190207"
-    prefix = "2019-02-07T08:26:59"
-    ifile = 0
-    nsub = 21600
 
     # Read spectrogram
     s = Spectrogram.read_bin_files(path, prefix, ifile, nsub, 4171)
+
+    s.correct_spectrogram(143050000.0, 200.0)
+
     
     zmin = 0.0
     zmax = 0.00008
-    zmax = 0.02
 
 
     # Start ppgplot
@@ -73,7 +71,7 @@ if __name__ == "__main__":
             fmin -= fcen
             fmax -= fcen
             ppgplot.pgswin(xmin, xmax, fmin, fmax)
-            ppgplot.pgbox("", 0., 0, "BCTSN", 0., 0)
+            ppgplot.pgbox("", 0., 0, "BTSN", 0., 0)
             ppgplot.pglab("UT Date", "Frequency - {:.3f} MHz".format(fcen), "")
 
             # Back to pixels
diff --git a/python/strf/cmap.py b/python/strf/cmap.py
index aff9c76..f363209 100644
--- a/python/strf/cmap.py
+++ b/python/strf/cmap.py
@@ -1,4 +1,4 @@
-#!/usr/bin/env python3
+#!/usr/bin/env python
 import numpy as np
 
 cm_cool = np.array([[-0.5, 0.0, 0.17, 0.33, 0.50, 0.67, 0.83, 1.0, 1.7], 
diff --git a/python/strf/rfio.py b/python/strf/rfio.py
index 76dbfd4..bc43256 100644
--- a/python/strf/rfio.py
+++ b/python/strf/rfio.py
@@ -1,9 +1,10 @@
-#!/usr/bin/env python3
+#!/usr/bin/env python
 import numpy as np
 from datetime import datetime
 import re
 import os
 from astropy.time import Time
+import scipy.optimize
 
 class Spectrogram:
     """Spectrogram class"""
@@ -34,7 +35,7 @@ class Spectrogram:
         while isub<nsub:
             # File name of file
             fname = os.path.join(path, "{:s}_{:06d}.bin".format(prefix, ifile))
-            print(fname)
+#            print(fname)
             with open(fname, "rb") as fp:
                 next_header = fp.read(256)
                 while next_header:
@@ -48,6 +49,17 @@ class Spectrogram:
             
         return Spectrogram(np.transpose(np.vstack(zs)), np.array(mjds), freq, siteid)
 
+    def correct_spectrogram(self, fcen, bw):
+        c = np.abs(self.freq-fcen)<bw
+        x = self.freq[c]-fcen
+        y = np.median(self.z, axis=1)[c]
+
+        p = [x[np.argmax(y)], 5.0, 1.0, 0.0]
+
+        q, cov_q = scipy.optimize.leastsq(gaussian_residual, p, args=(x, y))
+
+        self.freq -= q[0]
+    
 def parse_header(header_b):
     header_s = header_b.decode('ASCII').strip('\x00')
     regex = r"^HEADER\nUTC_START    (.*)\nFREQ         (.*) Hz\nBW           (.*) Hz\nLENGTH       (.*) s\nNCHAN        (.*)\nNSUB         (.*)\nEND\n$"
@@ -62,3 +74,10 @@ def parse_header(header_b):
             'nchan': int(match.group(5)),
             'nsub': int(match.group(6))}
     
+def gaussian(x, a):
+    arg = -0.5*((x-a[0])/a[1])**2
+    return a[2]*np.exp(arg)+a[3]
+
+def gaussian_residual(a, x, y):
+    ym = gaussian(x, a)
+    return y-ym