eclipswaarneming.py

#!/usr/bin/env python3

"""
Receive numbers from gnuradio (over UDP) at about 2 per second, plot them
against time. Special color if not pointing to sun

Tammo Jan Dijkema, 29 June 2020
"""


import time
import matplotlib.pyplot as plt
import datetime
import sys
import zmq
import struct
import pandas as pd
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
import matplotlib.dates as mdates
from threading import Thread

plt.ion()

import logging
logger = logging.getLogger()
logger.setLevel(level=logging.INFO)

signal = 0.
prev_signal = 0.

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')
dt.getAzEl(waitForUpdate=True)

dwl = EarthLocation(lon=Angle("6:23:46.21 degrees"),
                    lat=Angle("52:48:43.27 degrees"))

# create plot
plt.ion()
fig, ax = plt.subplots()

#ax.set_ylim((0, None))
ax.set_xlabel("Lokale tijd (CEST)")
ax.set_ylabel("Signaalsterkte (ongekalibreerd)")
fig.canvas.set_window_title('Live radiowaarneming gedeeltelijke zonsverduistering')
ax.set_title('Live radiowaarneming gedeeltelijke zonsverduistering')

prev_points = None
separation = None
maxsignal = 0.

file_data_times = []
file_data = []
file_cal_times = []
file_caldata = []
starttime = Time.now() + 2 * u.hour

if len(sys.argv) > 1:
    df = pd.read_csv(sys.argv[1], parse_dates=[0])
    is_data = (df["tracker"] == 'suntracker') & \
              (df["offset az (deg)"] < 0.01) & \
              (df["offset el (deg)"] < 0.01)
    file_data_times = (Time(df[is_data]["time (UTC)"]) + 2*u.hour).datetime
    file_data = df[is_data]["signal"]
    file_cal_times = (Time(df[~is_data]["time (UTC)"]) + 2*u.hour).datetime
    file_caldata = df[~is_data]["signal"]
    starttime = Time(df["time (UTC)"][0]) + 2 * u.hour

scatterplot = ax.scatter(file_data_times, file_data, s=3, color='black')
calplot = ax.scatter(file_cal_times, file_caldata, s=3, color='lightgrey')
ax.xaxis.set_major_formatter(mdates.DateFormatter('%H:%M'))

eclipse_starttime = Time("2021-06-10T11:21") + 2*u.hour
eclipse_endtime = Time("2021-06-10T13:34") + 2*u.hour
if starttime < eclipse_starttime < starttime + 3*u.hour:
    ax.axvline(eclipse_starttime.datetime, linestyle='--', linewidth=1, color='lightgrey')
    ax.text(eclipse_starttime.datetime,0.1,'Start: '+eclipse_starttime.isot[11:16],rotation=90, color='black', transform=ax.get_xaxis_transform())
if starttime < eclipse_endtime < starttime + 3*u.hour:
    ax.axvline(eclipse_endtime.datetime, linestyle='--', linewidth=1, color='lightgrey')
    ax.text(eclipse_endtime.datetime,0.1,'Einde: '+eclipse_endtime.isot[11:16],rotation=90, color='black', transform=ax.get_xaxis_transform())

if len(sys.argv) == 1:
    datafile = open(f"data{int(Time.now().unix)}.txt", "w")
    print("time (UTC),tracker,offset az (deg),offset el (deg),azimuth from north (deg),elevation (deg),signal", file=datafile)
else:
    datafile = open(f"{sys.argv[1]}", "a")

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

    now = Time.now()

    tracker = dt.tracker.split("/")[-1]

    print(now, tracker, dt.offset_az.to(u.deg).value, dt.offset_el.to(u.deg).value, dt.az.to(u.deg).value, dt.el.to(u.deg).value, signal, flush=True, sep=',', file=datafile)

    if signal > maxsignal:
        maxsignal = signal

    if tracker == "suntracker" and dt.dist_az < 0.01 * u.deg and dt.dist_el < 0.01 * u.deg:
        plot_to_append_to = scatterplot
    else:
        plot_to_append_to = calplot

    scatter_data = plot_to_append_to.get_offsets()
    scatter_data = np.vstack((scatter_data, [mdates.date2num((now + 2*u.hour).datetime), signal]))
    plot_to_append_to.set_offsets(scatter_data)

    if prev_points is not None:
        prev_points.remove()

    prev_points = ax.scatter([(now+2*u.hour).datetime], [signal], s=50, color='red')
    ax.set_ylim((0.5 * maxsignal, 1.1 * maxsignal))
    ax.set_xlim(starttime.datetime, (now+2.5*u.hour).datetime) # HIERO

    # show the plot
    plt.show()
    plt.pause(0.0001)