Add many sources to zonnewaarneming

zonnewaarneming.py

@@ -14,7 +14,7 @@ import pandas as pd
 import zmq
 import struct
 import numpy as np
-from telescope import telescope
+from telescope import Telescope
 from astropy.coordinates import get_sun, SkyCoord, AltAz, EarthLocation, Angle
 from astropy.time import Time
 import astropy.units as u
@@ -27,11 +27,17 @@ import logging
 logger = logging.getLogger()
 logger.setLevel(level=logging.INFO)
 
+stoprequest = False
+
+def set_stoprequest(event):
+    global stoprequest
+    stoprequest = True
+
 signal = 0.
 prev_signal = 0.
 
 datafile = open(f"data{int(Time.now().unix)}.txt", "w")
-print("Time (UTC),separation (deg),azimuth from north (deg), elevation (deg), signal", file=datafile)
+print("Time (UTC),separation (deg),azimuth from north (deg), elevation (deg), signal, source", file=datafile)
 
 def signal_thread():
     global signal
@@ -46,38 +52,43 @@ s = Thread(target=signal_thread)
 s.daemon = True
 s.start()
 
-dt = telescope(consoleHost='console')
+dt = Telescope(consoleHost='localhost')
 
 cas_a = SkyCoord.from_name("Cassiopeia A")
 dwl = EarthLocation(lon=Angle("6:23:46.21 degrees"),
                     lat=Angle("52:48:43.27 degrees"))
 
-SEPARATION_AZ = True
+sourcenames = ["Cassiopeia A", "Cygnus A", "Perseus A", "Taurus A", "Zon"]
+sources = {}
+for sourcename in sourcenames:
+    if sourcename == "Zon":
+        continue
+    sources[sourcename] = SkyCoord.from_name(sourcename)
 
 # create plot
 plt.ion()
 fig, ax = plt.subplots()
-if SEPARATION_AZ:
-    ax.set_xlim((-6, 6))
-else:
-    ax.set_xlim((0, 6))
+ax.set_xlim((-6, 6))
+fig.canvas.mpl_connect('close_event', set_stoprequest)
 
 ax.axvline(x=0, linestyle='--', color='grey', linewidth=.5)
 
+scatterplots = {}
+for sourcename in sourcenames:
+    scatterplots[sourcename] = ax.scatter([], [], s=3, label=sourcename)
+
 #ax.set_ylim((0, None))
 ax.set_xlabel("Afstand tot middelpunt (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')
+ax.legend()
 
 prev_points = None
 separation = None
 maxsignal = 0.
 
-scatterplot_sun = ax.scatter([], [], s=3, color='black')
-scatterplot_cas = ax.scatter([], [], s=3, color='blue')
-
-while True:
+while not stoprequest:
     time.sleep(0.2)
     if signal == prev_signal:
         continue
@@ -85,51 +96,54 @@ while True:
 
     pos = dt.radec
     az = dt.az
-    el = dt.el
+    alt = dt.el
     if pos is None or az is None:
         print("No azimuth received!")
         continue
 
     az = Angle(az + 180 * u.deg)
 
+    sources["Zon"] = get_sun(Time.now())
+
     color = 'black'
     dwl_altaz = AltAz(location=dwl, obstime=Time.now())
-    altaz_sun = get_sun(Time.now()).transform_to(dwl_altaz)
-    az_sun, alt_sun = altaz_sun.az, altaz_sun.alt
-    altaz_cas = cas_a.transform_to(dwl_altaz)
-    az_cas, alt_cas = altaz_cas.az, altaz_cas.alt
-    separation_sun = (az - az_sun).wrap_at(180*u.deg).degree * np.cos(el)
-    separation_cas = (az - az_cas).wrap_at(180*u.deg).degree * np.cos(el)
-
-    if min(abs((alt_sun - el).degree), abs((alt_cas - el).degree)) > 0.2:
-        print("Elevation off, skipping")
+
+    separation = 180*u.deg
+    for sourcename in sources:
+        source_separation = sources[sourcename].separation(pos)
+        if source_separation < separation:
+            closest_source = sourcename
+            separation = source_separation
+
+    if prev_points is not None:
+        prev_points.remove()
+
+    if separation > 6 * u.deg:
+        print("Far from everything, plotting nothing")
         continue
 
-    if abs(separation_sun) < abs(separation_cas):
-        separation = separation_sun
-        scatterplot = scatterplot_sun
+    source_altaz = sources[closest_source].transform_to(dwl_altaz)
+    separation_az = (az - source_altaz.az).wrap_at(180*u.deg)
+    separation_alt = (alt - source_altaz.alt)
+    if np.abs(separation_az) > np.abs(separation_alt):
+        if separation_az < 0*u.deg:
+            separation = -separation
     else:
-        separation = separation_cas
-        scatterplot = scatterplot_cas
+        if separation_alt < 0*u.deg:
+            separation = -separation
 
-    print(Time.now(), separation, az.to(u.deg).value, el.to(u.deg).value, signal, file=datafile, flush=True, sep=',')
+    scatterplot = scatterplots[closest_source]
+
+    print(Time.now(), separation, az.to(u.deg).value, alt.to(u.deg).value, signal, closest_source, file=datafile, flush=True, sep=',')
 
     if signal > maxsignal:
         maxsignal = signal
 
-    if abs(separation) > 10:
-        plt.pause(2)
-        print(f"Out of range, separation is {separation:.1f}")
-        continue
-
     scatter_data = scatterplot.get_offsets()
-    scatter_data = np.vstack((scatter_data, [separation, signal]))
+    scatter_data = np.vstack((scatter_data, [separation.deg, signal]))
     scatterplot.set_offsets(scatter_data)
 
-    if prev_points is not None:
-        prev_points.remove()
-
-    prev_points = ax.scatter([separation], [signal], s=50, color='red')
+    prev_points = ax.scatter([separation.deg], [signal], s=50, color='red')
     ax.set_ylim((2800, 1.1 * maxsignal))
 
     # show the plot