zonnewaarneming.py 4.85 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
#!/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
Tammo Jan Dijkema's avatar
Tammo Jan Dijkema committed
14
import zmq
Tammo Jan Dijkema's avatar
Tammo Jan Dijkema committed
15
import re
Tammo Jan Dijkema's avatar
Tammo Jan Dijkema committed
16
17
import struct
import numpy as np
18
from telescope import Telescope
Tammo Jan Dijkema's avatar
Tammo Jan Dijkema committed
19
20
21
22
from astropy.coordinates import get_sun, SkyCoord, AltAz, EarthLocation, Angle
from astropy.time import Time
import astropy.units as u
from matplotlib import cm
Tammo Jan Dijkema's avatar
Tammo Jan Dijkema committed
23
from threading import Thread
Tammo Jan Dijkema's avatar
Tammo Jan Dijkema committed
24

Tammo Jan Dijkema's avatar
Tammo Jan Dijkema committed
25
plt.ion()
Tammo Jan Dijkema's avatar
Tammo Jan Dijkema committed
26

Tammo Jan Dijkema's avatar
Tammo Jan Dijkema committed
27
import logging
Tammo Jan Dijkema's avatar
Tammo Jan Dijkema committed
28
29
30
logger = logging.getLogger()
logger.setLevel(level=logging.INFO)

31
32
33
34
35
36
stoprequest = False

def set_stoprequest(event):
    global stoprequest
    stoprequest = True

Tammo Jan Dijkema's avatar
Tammo Jan Dijkema committed
37
38
signal = 0.
prev_signal = 0.
Tammo Jan Dijkema's avatar
Tammo Jan Dijkema committed
39

Tammo Jan Dijkema's avatar
Tammo Jan Dijkema committed
40
datafile = open(f"data{int(Time.now().unix)}.txt", "w")
41
print("Time (UTC),separation (deg),azimuth from north (deg), elevation (deg), signal, source", file=datafile)
Tammo Jan Dijkema's avatar
Tammo Jan Dijkema committed
42

Tammo Jan Dijkema's avatar
Tammo Jan Dijkema committed
43
44
45
46
47
48
49
50
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]
Tammo Jan Dijkema's avatar
Tammo Jan Dijkema committed
51

Tammo Jan Dijkema's avatar
Tammo Jan Dijkema committed
52
53
54
55
s = Thread(target=signal_thread)
s.daemon = True
s.start()

56
dt = Telescope(consoleHost='localhost')
Tammo Jan Dijkema's avatar
Tammo Jan Dijkema committed
57
58
59
60
61

cas_a = SkyCoord.from_name("Cassiopeia A")
dwl = EarthLocation(lon=Angle("6:23:46.21 degrees"),
                    lat=Angle("52:48:43.27 degrees"))

Tammo Jan Dijkema's avatar
Tammo Jan Dijkema committed
62
63
64
65
66
67
68
69
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)

70
71
72
73
sources = {}
for sourcename in sourcenames:
    if sourcename == "Zon":
        continue
Tammo Jan Dijkema's avatar
Tammo Jan Dijkema committed
74
    sources[sourcename] = SkyCoord.from_name(clean_name(sourcename))
Tammo Jan Dijkema's avatar
Tammo Jan Dijkema committed
75

Tammo Jan Dijkema's avatar
Tammo Jan Dijkema committed
76
77
78
# create plot
plt.ion()
fig, ax = plt.subplots()
79
80
ax.set_xlim((-6, 6))
fig.canvas.mpl_connect('close_event', set_stoprequest)
Tammo Jan Dijkema's avatar
Tammo Jan Dijkema committed
81

Tammo Jan Dijkema's avatar
Tammo Jan Dijkema committed
82
83
ax.axvline(x=0, linestyle='--', color='grey', linewidth=.5)

84
85
scatterplots = {}
for sourcename in sourcenames:
Tammo Jan Dijkema's avatar
Tammo Jan Dijkema committed
86
    scatterplots[sourcename] = ax.scatter([], [], s=3, label=None)
87

Tammo Jan Dijkema's avatar
Tammo Jan Dijkema committed
88
#ax.set_ylim((0, None))
Tammo Jan Dijkema's avatar
Tammo Jan Dijkema committed
89
ax.set_xlabel("Afstand tot middelpunt (graden)")
Tammo Jan Dijkema's avatar
Tammo Jan Dijkema committed
90
91
92
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')
Tammo Jan Dijkema's avatar
Tammo Jan Dijkema committed
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110

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)
Tammo Jan Dijkema's avatar
Tammo Jan Dijkema committed
111

Tammo Jan Dijkema's avatar
Tammo Jan Dijkema committed
112
113
114
115
prev_points = None
separation = None
maxsignal = 0.

116
while not stoprequest:
Tammo Jan Dijkema's avatar
Tammo Jan Dijkema committed
117
118
119
120
    time.sleep(0.2)
    if signal == prev_signal:
        continue
    prev_signal = signal
Tammo Jan Dijkema's avatar
Tammo Jan Dijkema committed
121
122
123

    pos = dt.radec
    az = dt.az
124
    alt = dt.el
Tammo Jan Dijkema's avatar
Tammo Jan Dijkema committed
125
    if pos is None or az is None:
Tammo Jan Dijkema's avatar
Tammo Jan Dijkema committed
126
        print("No azimuth received!")
Tammo Jan Dijkema's avatar
Tammo Jan Dijkema committed
127
128
129
130
        continue

    az = Angle(az + 180 * u.deg)

131
132
    sources["Zon"] = get_sun(Time.now())

Tammo Jan Dijkema's avatar
Tammo Jan Dijkema committed
133
    color = 'black'
134
    dwl_altaz = AltAz(location=dwl, obstime=Time.now())
135
136
137
138
139
140
141
142
143

    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:
Tammo Jan Dijkema's avatar
Tammo Jan Dijkema committed
144
145
146
147
        try:
            prev_points.remove()
        except:
            pass
148
149

    if separation > 6 * u.deg:
150
151
        continue

152
153
154
155
156
157
    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
Tammo Jan Dijkema's avatar
Tammo Jan Dijkema committed
158
    else:
159
160
        if separation_alt < 0*u.deg:
            separation = -separation
Tammo Jan Dijkema's avatar
Tammo Jan Dijkema committed
161

162
163
    scatterplot = scatterplots[closest_source]

Tammo Jan Dijkema's avatar
Tammo Jan Dijkema committed
164
    print(Time.now(), separation, az.to(u.deg).value, alt.to(u.deg).value, signal, clean_name(closest_source), file=datafile, flush=True, sep=',')
Tammo Jan Dijkema's avatar
Tammo Jan Dijkema committed
165
166
167
168

    if signal > maxsignal:
        maxsignal = signal

169
    scatter_data = scatterplot.get_offsets()
170
    scatter_data = np.vstack((scatter_data, [separation.deg, signal]))
171
    scatterplot.set_offsets(scatter_data)
Tammo Jan Dijkema's avatar
Tammo Jan Dijkema committed
172

173
    prev_points = ax.scatter([separation.deg], [signal], s=50, color='red')
174
    ax.set_ylim((2800, 1.1 * maxsignal))
Tammo Jan Dijkema's avatar
Tammo Jan Dijkema committed
175

Tammo Jan Dijkema's avatar
Tammo Jan Dijkema committed
176
177
    update_legend(ax, closest_source)

Tammo Jan Dijkema's avatar
Tammo Jan Dijkema committed
178
179
    # show the plot
    plt.show()
Tammo Jan Dijkema's avatar
Tammo Jan Dijkema committed
180

181
    plt.pause(0.0001)