Remove import of the beam computation function.

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