Initial commit

zonnewaarneming.py

+#!/usr/bin/env python3
+
+"""
+Receive numbers from gnuradio (over UDP) at about 2 per second, plot them
+against distance to to Sun or Cas A
+
+Tammo Jan Dijkema, 29 June 2020
+"""
+
+
+import time
+import matplotlib.pyplot as plt
+import pandas as pd
+import socket
+import struct
+import numpy as np
+from telescope import telescope
+from astropy.coordinates import get_sun, SkyCoord, AltAz, EarthLocation, Angle
+from astropy.time import Time
+import astropy.units as u
+from matplotlib import cm
+
+colormap = cm.Greys
+
+dt = telescope(consoleHost='console')
+
+cas_a = SkyCoord.from_name("Cassiopeia A")
+dwl = EarthLocation(lon=Angle("6:23:46.21 degrees"),
+                    lat=Angle("52:48:43.27 degrees"))
+
+UDP_IP = "127.0.0.1"
+UDP_PORT = 7805
+
+sock = socket.socket(socket.AF_INET, # Internet
+                     socket.SOCK_DGRAM) # UDP
+sock.bind((UDP_IP, UDP_PORT))
+
+# create an empty dataframe that will store streaming data
+df = pd.DataFrame(np.zeros((0, 2)), columns=["Separation", "Signal"])
+
+# create plot
+plt.ion()
+fig, ax = plt.subplots()
+ax.set_xlim((-10, 10))
+ax.set_xlabel("Afstand tot zon (graden)")
+ax.set_ylabel("Signaalsterkte (ongekalibreerd)")
+fig.canvas.set_window_title('Live waarneming Dwingeloo radiotelescoop op 1330 MHz')
+ax.set_title('Live waarneming Dwingeloo radiotelescoop op 1330 MHz')
+
+SEPARATION_AZ = True
+
+while True:
+    data, addr = sock.recvfrom(1024)
+    if len(data) != 4:
+        continue
+
+    signal = struct.unpack('f', data)[0]
+
+    # append new row to dataframe
+    pos = dt.radec
+    az = dt.az
+    if pos is None or az is None:
+        continue
+
+    az = Angle(az + 180 * u.deg)
+
+    if SEPARATION_AZ:
+        dwl_altaz = AltAz(location=dwl, obstime=Time.now())
+        az_sun = get_sun(Time.now()).transform_to(dwl_altaz).az
+        az_cas = cas_a.transform_to(dwl_altaz).az
+        separation_sun = (az - az_sun).wrap_at(180*u.deg).degree
+        separation_cas = (az - az_cas).wrap_at(180*u.deg).degree
+        if abs(separation_sun) < abs(separation_cas):
+            separation = separation_sun
+        else:
+            separation = separation_cas
+    else:
+        separation_sun = get_sun(Time.now()).separation(pos).deg
+        separation_cas = cas_a.separation(pos).deg
+        separation = min(separation_sun, separation_cas)
+
+    if abs(separation) > 10 and len(df) > 10:
+        plt.pause(2)
+        continue
+
+    df.loc[len(df)] = [separation, signal]
+
+    # plot all data
+    #ax.plot(df["Separation"], df["Signal"], 'k.')
+    ax.scatter(df["Separation"], df["Signal"], s=3, c=colormap(np.arange(len(df))**2))
+
+    # show the plot
+    plt.show()
+    plt.pause(0.0001) # <-- sets the current plot until refreshed
+
+    # be nice to the cpu :)
+    time.sleep(.1)