From 7e220aa7c9da3511fe46d7611be6b2813c696e5c Mon Sep 17 00:00:00 2001
From: Aleksandar Shulevski <shulevski@astron.nl>
Date: Thu, 1 Aug 2024 08:49:05 +0000
Subject: [PATCH] Remove import of the beam computation function.
---
plot_sidelobe_tracks.py | 47 ++++++++++++++++++++++++++++-------------
1 file changed, 32 insertions(+), 15 deletions(-)
diff --git a/plot_sidelobe_tracks.py b/plot_sidelobe_tracks.py
index 1fe3825..8661411 100644
--- a/plot_sidelobe_tracks.py
+++ b/plot_sidelobe_tracks.py
@@ -10,14 +10,25 @@ warnings.filterwarnings("ignore")
import lofarantpos.db
from lofarantpos.geo import localnorth_to_etrs, geographic_array_from_xyz, xyz_from_geographic, geographic_from_xyz
-from lofar_beams import voltage_beam
-
import ephem
from photutils.detection import DAOStarFinder
from astropy.stats import mad_std
from astropy import log
+def voltage_beam(lgrid, mgrid, lbeam, mbeam, c, p, q, r, freq, weights):
+ """Compute the beam response."""
+ # n-coordinates
+ ngrid = np.sqrt(1 - lgrid**2 - mgrid**2)
+ nbeam = np.sqrt(1 - lbeam**2 - mbeam**2)
+
+ # lmn-uvw product
+ u, v, w = np.array([p, q, r]) * freq / c
+ prod = 2j * np.pi * (u[:, np.newaxis, np.newaxis] * (lgrid - lbeam) +
+ v[:, np.newaxis, np.newaxis] * (mgrid - mbeam) +
+ w[:, np.newaxis, np.newaxis] * (ngrid - nbeam))
+ return np.sum(weights * np.exp(prod), axis=0)
+
def createFixedBody(ra, dec, epoch):
"""Fixed body adapted to us."""
fixedBody = ephem.FixedBody()
@@ -65,7 +76,7 @@ def compute_fixed_sources_altaz(station, end_of_day, source):
time.append(station.date.tuple()[3])
return az, el, time
-def get_dynspec_for_source(station, end_of_day, pointing, source, freq_arr, p, q, r):
+def get_dynspec_for_source(station, end_of_day, pointing, source, freq_arr, c, p, q, r):
"""Construct the dynamic spectrum for a station.
Parameters
@@ -78,7 +89,10 @@ def get_dynspec_for_source(station, end_of_day, pointing, source, freq_arr, p, q
RA, DEC of the pointing direction
source : str
Source name, 'CasA' etc.
- freq_arr : array of frequencies for which we compute the dynamic spectrum
+ freq_arr : array
+ Array of frequencies for which we compute the dynamic spectrum
+ c : float
+ Speed of light, m/s
p, q, r : floats
Antenna position parameters
@@ -105,7 +119,7 @@ def get_dynspec_for_source(station, end_of_day, pointing, source, freq_arr, p, q
l_src, m_src = np.sin(az_src * np.pi / 180) * np.cos(el_src * np.pi / 180), np.cos(az_src * np.pi / 180) * np.cos(el_src * np.pi / 180)
# Compute beam
log.debug('Compute beam for {} Hz'.format(fr))
- beam = voltage_beam(lgrid, mgrid, ltile, mtile, p, q, r, fr, 1.0)
+ beam = voltage_beam(lgrid, mgrid, ltile, mtile, c, p, q, r, fr, 1.0)
masked_beam = ma.masked_invalid(np.abs(beam))
# Select beam values around source position - is delta of 0.01 right?
# TODO: Normalize beam?
@@ -114,7 +128,7 @@ def get_dynspec_for_source(station, end_of_day, pointing, source, freq_arr, p, q
dynspec.append(time)
return dynspec
-def compute_lobes(station, end_of_day, pointing, min_elevation, freq, p, q, r):
+def compute_lobes(station, end_of_day, pointing, min_elevation, freq, c, p, q, r):
"""Compute the position of the grating lobes.
Parameters
@@ -123,11 +137,14 @@ def compute_lobes(station, end_of_day, pointing, min_elevation, freq, p, q, r):
'e.g. CS001'
end_of_day : date
End of the day specified
- pointing : Pointing direction
- RA, DEC
+ pointing : Object
+ Pointing direction
min_elevation : float
Minimum elevation below which we do not calculate, in degrees
- freq: frequency for which we are performing the calculation
+ freq: float
+ Frequency for which we are performing the calculation
+ c : float
+ Speed of light, m/s
p, q, r : floats
Antenna position parameters
@@ -154,7 +171,7 @@ def compute_lobes(station, end_of_day, pointing, min_elevation, freq, p, q, r):
y = (mtile + 1) * int(npix / 2)
# Compute beam
- beam = voltage_beam(lgrid, mgrid, ltile, mtile, p, q, r, freq, 1.0)
+ beam = voltage_beam(lgrid, mgrid, ltile, mtile, c, p, q, r, freq, 1.0)
masked_beam = ma.masked_invalid(np.abs(beam))
bkg_sigma = mad_std(masked_beam)
# Extract grating lobe maxima
@@ -474,13 +491,13 @@ if __name__ == "__main__":
ax[i].annotate(str(time[0])+'h', (az[0], el[0]), textcoords="offset points", xytext=(0,10), ha='center')
ax[i].annotate(str(time[-2])+'h', (az[-2], el[-2]), textcoords="offset points", xytext=(0,10), ha='center')
if dynspec_source == at.name:
- dynspec = get_dynspec_for_source(station, end_of_day, pointing_direction, at, freq_arr, p, q, r)
+ dynspec = get_dynspec_for_source(station, end_of_day, pointing_direction, at, freq_arr, c, p, q, r)
ax1[i].imshow(dynspec, origin='lower', extent=[time[0], time[-2], freq_arr[0] / 1.e6, freq_arr[-1] / 1.e6], aspect='equal')
else:
ax[i].scatter(ephem.degrees(at.body.ra), at.body.dec * (180. / np.pi), label=at.name)
if dynspec_source == at.name:
log.info("Calculating dynspec plot for {}".format(dynspec_source))
- dynspec = get_dynspec_for_source(station, end_of_day, pointing_direction, at, freq_arr, p, q, r)
+ dynspec = get_dynspec_for_source(station, end_of_day, pointing_direction, at, freq_arr, c, p, q, r)
ax1[i].imshow(dynspec, origin='lower', extent=[time[0], time[-2], freq_arr[0] / 1.e6, freq_arr[-1] / 1.e6], aspect='equal')
ax1[i].set_title("Dynspec for station " + station_name, pad=10)
ax1[i].set_xlabel("Time [UT]")
@@ -497,12 +514,12 @@ if __name__ == "__main__":
ax[i].annotate(str(time[0])+'h', (az[0], el[0]), textcoords="offset points", xytext=(0,10), ha='center')
ax[i].annotate(str(time[-2])+'h', (az[-2], el[-2]), textcoords="offset points", xytext=(0,10), ha='center')
if dynspec_source == sol.name:
- dynspec = get_dynspec_for_source(station, end_of_day, pointing_direction, sol, freq_arr, p, q, r)
+ dynspec = get_dynspec_for_source(station, end_of_day, pointing_direction, sol, freq_arr, c, p, q, r)
ax1[i].imshow(dynspec, origin='lower', extent=[time[0], time[-2], freq_arr[0] / 1.e6, freq_arr[-1] / 1.e6], aspect='equal')
else:
ax[i].scatter(ephem.degrees(sol.body.ra), sol.body.dec * (180. /np.pi), label=sol.name)
if dynspec_source == sol.name:
- dynspec = get_dynspec_for_source(station, end_of_day, pointing_direction, sol, freq_arr, p, q, r)
+ dynspec = get_dynspec_for_source(station, end_of_day, pointing_direction, sol, freq_arr, c, p, q, r)
ax1[i].imshow(dynspec, origin='lower', extent=[time[0], time[-2], freq_arr[0] / 1.e6, freq_arr[-1] / 1.e6], aspect='equal')
ax1[i].set_title("Dynspec for station " + station_name, pad=10)
ax1[i].set_xlabel("Time [UT]")
@@ -513,7 +530,7 @@ if __name__ == "__main__":
station.date = start_date
end_of_day = station.date
end_of_day += 1.
- long, lat, time = compute_lobes(station, end_of_day, pointing_direction, min_elevation, freq * 1.e6, p, q, r)
+ long, lat, time = compute_lobes(station, end_of_day, pointing_direction, min_elevation, freq * 1.e6, c, p, q, r)
if len(long) != 0 and len(lat) != 0 and len(time) != 0:
if altaz:
ax[i].scatter(long, lat, marker='x', s=40, label='lobe: ' + str(freq) + ' MHz')
--
GitLab