eclipswaarneming.py 2.53 KB
Newer Older
Tammo Jan Dijkema's avatar
Tammo Jan Dijkema committed
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
#!/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 zmq
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
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.

datafile = open(f"data{int(Time.now().unix)}.txt", "w")
print("Time (UTC),separation (deg),azimuth from north (deg), elevation (deg), signal", 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')
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.

scatterplot_70cm = ax.scatter([], [], s=3, color='black')
scatterplot_23cm = ax.scatter([], [], s=3, color='blue')

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

    tracker = dt.tracker
    print("Tracker is ", tracker)
    now = Time.now()

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

    if signal > maxsignal:
        maxsignal = signal

    scatter_data = scatterplot_70cm.get_offsets()
    print(scatter_data)
    scatter_data = np.vstack((scatter_data, [mdates.date2num(now.datetime), signal]))
    print(scatter_data)
    scatterplot_70cm.set_offsets(scatter_data)

    if prev_points is not None:
        prev_points.remove()

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

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