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 zmq
import re
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
from threading import Thread

plt.ion()

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, source", file=datafile)

def signal_thread():
    global signal
    context = zmq.Context()
    receiver = context.socket(zmq.PULL)
    receiver.connect("tcp://127.0.0.1:5557")
    while True:
        buff = receiver.recv()
        signal = np.frombuffer(buff, dtype="float32")[-1]

s = Thread(target=signal_thread)
s.daemon = True
s.start()

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"))

sourcenames = ["Cassiopeia A\n(supernovarestant)", "Cygnus A\n(sterrenstelsel)", "Perseus A\n(sterrenstelsel)", "Taurus A (Krabnevel)", "Andromeda (sterrenstelsel)", "Zon"]

pattern = re.compile(r'[ \n]+\(.*\)')
def clean_name(prettyname):
    """Strip things between brackets"""
    #return re.sub(r'[ \n]+\(.*\)', '', prettyname)
    return pattern.sub('', prettyname)

sources = {}
for sourcename in sourcenames:
    if sourcename == "Zon":
        continue
    sources[sourcename] = SkyCoord.from_name(clean_name(sourcename))

# create plot
plt.ion()
fig, ax = plt.subplots()
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=None)

#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')

def update_legend(ax, sourcename):
    somethingchanged = False
    scatterplot = scatterplots[sourcename]
    if scatterplot.get_label()[0] == "_":
        num_points = scatterplot.get_offsets().data.shape[0]
        if num_points > 5:
            scatterplot.set_label(sourcename)
            somethingchanged = True
            print("Adding", sourcename.replace("\n", " "), "to legend")
    if somethingchanged:
        legend = ax.get_legend()
        if legend is not None:
            legend.remove()
        ax.legend()

plt.show()
plt.pause(0.001)

prev_points = None
separation = None
maxsignal = 0.

while not stoprequest:
    time.sleep(0.2)
    if signal == prev_signal:
        continue
    prev_signal = signal

    pos = dt.radec
    az = dt.az
    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())

    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:
        try:
            prev_points.remove()
        except:
            pass

    if separation > 6 * u.deg:
        continue

    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:
        if separation_alt < 0*u.deg:
            separation = -separation

    scatterplot = scatterplots[closest_source]

    print(Time.now(), separation, az.to(u.deg).value, alt.to(u.deg).value, signal, clean_name(closest_source), file=datafile, flush=True, sep=',')

    if signal > maxsignal:
        maxsignal = signal

    scatter_data = scatterplot.get_offsets()
    scatter_data = np.vstack((scatter_data, [separation.deg, signal]))
    scatterplot.set_offsets(scatter_data)

    prev_points = ax.scatter([separation.deg], [signal], s=50, color='red')
    ax.set_ylim((2800, 1.1 * maxsignal))

    update_legend(ax, closest_source)

    # show the plot
    plt.show()

    plt.pause(0.0001)