<jupyter_start><jupyter_text>Read MeerKAT config file with antenna positions<jupyter_code>class MKAT(object):
    def __init__(self, mkat_config):
        self.ref_position = None
        arr_desc = {'names': ('name', 'east', 'north', 'up', 'location'),
                    'formats': ('S5', float, float, float, EarthLocation)}
        self.antennas = recarray((0,), dtype=arr_desc)
        self.build_array(mkat_config)
        
        self.frequency = 1420e6 / u.s
        self.wavelength = phys.c / self.frequency
        self.nr_antennas = len(self.antennas)
        self.nr_baselines = (self.nr_antennas * (self.nr_antennas - 1)) // 2
        
        self.observer = self.ephem_observer()
        
    def build_array(self, filename):
        """Read array default .yaml file."""
        with open(filename, "r") as stream:
            data = yaml.safe_load(stream)
        self.get_reference(data['reference'])
        self.get_antennas(data['antennas'])
 
    def get_reference(self, ref_position):
        """General Earth position for the Telescope"""
        ref_LAT = Latitude(ref_position['latitude'],
                           u.degree,
                           copy=False)
        ref_LON = Longitude(ref_position['longitude'],
                            u.degree,
                            wrap_angle=180. * u.degree,
                            copy=False)
        ref_ALT = u.Quantity(ref_position['altitude'],
                             u.m,
                             copy=False)
        self.ref_position = EarthLocation(lat=ref_LAT,
                                          lon=ref_LON,
                                          height=ref_ALT)
 
    def get_antennas(self, antennas):
        """ENU coordinates per antenna"""
        if self.ref_position is not None:
            [x, y, z] = self.ref_position.to_geocentric()

        self.antennas.resize(len(antennas))
        for cnt, antenna in enumerate(antennas):
            ant_ = [item.strip() for item in antenna.split(",")]
            for item_ in ant_:
                for prefix_ in ('name', 'east', 'north', 'up'):
                    val_ =  None
                    if item_.startswith(prefix_):
                        val_ =  item_.split('=')[-1].strip()
                    if val_ is not None:
                        self.antennas[cnt][prefix_] = val_
            if self.ref_position is not None:
                ant_North = x.value + self.antennas[cnt]['north']
                ant_East = y.value + self.antennas[cnt]['east']
                ant_Up = z.value + self.antennas[cnt]['up']
                self.antennas[cnt]['location'] = EarthLocation(x=ant_North * u.m,
                                                               y=ant_East * u.m,
                                                               z=ant_Up * u.m)

    def ephem_observer(self):
        lon = self.ref_position.geodetic[0].degree
        lat = self.ref_position.geodetic[1].degree
        alt = self.ref_position.geodetic[2].value
        observer = ephem.Observer()
        observer.lon = str(lon)
        observer.lat = str(lat)
        observer.elevation = alt
        return observer
    
    def baselines(self):
        """The function evaluates baselines lenghts
           and the angle between antennas
        """
        P = np.array([self.antennas['north'],
                      self.antennas['east']]).T
        # antenna position in wavelength units
        P = P / self.wavelength.value  # baseline

        bl_length = np.zeros((self.nr_baselines, ))
        bl_az_angle = np.zeros((self.nr_baselines, ))
        cnt = 0
        for idx0 in range(self.nr_antennas):
            for idx1 in range(idx0 + 1, self.nr_antennas):
                bl_len_p0 = (P[idx0, 0] - P[idx1, 0])**2
                bl_len_p1 = (P[idx0, 1] - P[idx1, 1])**2
                bl_length[cnt] = np.sqrt(bl_len_p0 + bl_len_p1)
                bl_az_angle[cnt] = np.arctan2((P[idx0, 1] - P[idx1, 1]),
                                              (P[idx0, 0] - P[idx1, 0]))
                cnt += 1

        return bl_length, bl_az_angle
mkat = MKAT('mkat_antennas.yml')

def decode(text_):
    if isinstance(text_, str):
        text = text_
    else:
        text = text_.decode("utf-8")
    return text

# Get array
arr_latitude = []
arr_longitude = []
arr_names = []
for antenna in mkat.antennas:
    arr_names.append(decode(antenna['name']))
    arr_latitude.append(antenna['location'].lat.value)
    arr_longitude.append(antenna['location'].lon.value)<jupyter_output><empty_output><jupyter_text>Plotting functionality<jupyter_code>class UVplot(object):
    def __init__(self,
                 latitude,
                 declination,
                 elevation):
        self.lat = latitude
        self.dec = declination
        self.elev = elevation

    def __baseline_to_xyz__(self,
                            baseline_lengths,
                            baseline_azimuth):
        """The following function transform baseline to x,y,z coordinates
            elevation: Elevation angle in radian
            The result is a vector of (x,y,z) with unit the same as baseline length
            Interferometry and Synthesis in Radio Astronomy, Chapter 4, Equation 4.4
        """
        x0 = np.cos(self.lat) * np.sin(baseline_azimuth)
        x1 = np.sin(self.lat) * np.cos(baseline_azimuth) * np.cos(self.elev)
        x = x0 - x1
        y = np.cos(baseline_azimuth) * np.sin(self.elev)
        z0 = np.sin(self.lat) * np.sin(baseline_azimuth)
        z1 = np.cos(self.lat) * np.cos(baseline_azimuth) * np.cos(self.elev)
        z = z0 + z1
        xyz = np.array([(x, y, z)])
        return baseline_lengths * xyz.T

    def __xyz_to_baseline__(self,
                            hour_angle):
        """Transform x, y, z to u, v, w components
           hour_angle: Source hour angle in radians
           dec: Source declination in radian
           Interferometry and Synthesis in Radio Astronomy, Chapter 4, Equation 4.1
        """
        a1 = np.sin(hour_angle)
        a2 = np.cos(hour_angle)
        a3 = 0.

        b1 = -1 * np.sin(self.dec) * np.cos(hour_angle)
        b2 = np.sin(self.dec) * np.sin(hour_angle)
        b3 = np.cos(self.dec)

        c1 = np.cos(self.dec) * np.cos(hour_angle)
        c2 = -1 * np.cos(self.dec) * np.sin(hour_angle)
        c3 = np.sin(self.dec)

        return np.array([(a1, a2, a3),
                         (b1, b2, b3),
                         (c1, c2, c3)])

    def track_uv(self,
                 ha_range,
                 bl_length,
                 bl_azimuth,
                 timeslots):
        """Evaluate the track of a single antenna pair
           The function return a set of u,v,w components
        """
        uvw = np.zeros((timeslots, 3), dtype=float)
        for i in range(timeslots):
            uvw[i, :] = np.dot(self.__xyz_to_baseline__(ha_range[i]),
                               self.__baseline_to_xyz__(bl_length,
                                                        bl_azimuth)).T
        return uvw

    def uvMask(self,
               ha_range,
               bl_length,
               bl_azimuth,
               timeslots,
               maxsize,
               uvscaling):
        maskmat = np.zeros((maxsize, maxsize))
        uvw = self.track_uv(ha_range,
                            bl_length,
                            bl_azimuth,
                            timeslots)
        sctrl = maxsize // 2 + 1
        for i in range(timeslots):
            maskmat[sctrl + int(np.ceil(uvw[i, 0] * uvscaling)),
                    sctrl + int(np.ceil(uvw[i, 1] * uvscaling))] = 1.
            maskmat[sctrl - int(np.ceil(uvw[i, 0] * uvscaling)),
                    sctrl - int(np.ceil(uvw[i, 1] * uvscaling))] = 1.
        return maskmat

    def plot_uv(self,
                ha_range,
                bl_length,
                bl_azimuth,
                timeslots,
                comment=''):
        """UV coverage graph"""
        nr_baselines = len(bl_length)
        fig, ax = plt.subplots(nrows=1, ncols=1,
                               figsize=(20, 13),  # W x H
                               facecolor='white')
        for i in range(nr_baselines):
            uv = self.track_uv(ha_range,
                               bl_length[i],
                               bl_azimuth[i],
                               timeslots)
            ax.plot(uv[:, 0], uv[:, 1], 'b.', markersize=2)
            ax.plot(-uv[:, 0], -uv[:, 1], 'r.', markersize=2)
        ax.set_xlabel('u')
        ax.set_ylabel('v')
        ax.set_title('Declination {} deg: UV coverage'
                     .format(comment))
        return fig, ax

    def plot_mask(self,
                  mask,
                  psf,
                  comment=''):
        """Display simulated beamshape"""
        fig, ax = plt.subplots(nrows=1, ncols=2,
                               figsize=(20, 13),  # W x H
                               facecolor='white')
        ax[0].imshow(mask)
        [nr, nc] = np.shape(mask)
        ax[0].set_xlim(0.5 * nr - 0.25 * nr,
                       0.5 * nr + 0.25 * nr)
        ax[0].set_ylim(0.5 * nc - 0.25 * nc,
                       0.5 * nc + 0.25 * nc)
        ax[0].set_xlabel('u')
        ax[0].set_ylabel('v')
        ax[0].set_title('uv coverage, "%s"' % comment)

        ax[1].imshow(psf,
                     interpolation='gaussian',
                     origin='lower',
                     vmin=psf.min() / 1000.,
                     vmax=psf.max())
        [nr, nc] = np.shape(psf)
        ax[1].set_xlim(0.5 * nr - 0.25 * nr,
                       0.5 * nr + 0.25 * nr)
        ax[1].set_ylim(0.5 * nc - 0.25 * nc,
                       0.5 * nc + 0.25 * nc)
        ax[1].set_xlabel('u')
        ax[1].set_ylabel('v')
        ax[1].set_title('psf, "%s"' % comment)

        return fig, ax
# hour angle range in hours given nr of timeslots
# At meridian the hour angle is 0.000 degree,
# with the time before meridian crossing expressed as negative degrees,
# and the local time after expressed as positive degrees.
def get_ha_range(telescope,
                 target,
                 start_time,
                 duration,
                 ntimeslots=50):

    start_time = datetime.strptime(start_time, '%Y-%m-%d %H:%M:%S')
    stop_time = start_time + timedelta(seconds=duration)

    telescope.observer.date = ephem.Date(start_time)
    star = ephem.FixedBody()
    star._ra = target.ra.degree
    star._dec = target.dec.degree
    rise_time = telescope.observer.previous_rising(star)
    set_time = telescope.observer.next_setting(star,
                                               start=rise_time).datetime()
    if start_time > set_time:
        rise_time = telescope.observer.next_rising(star)
        set_time = telescope.observer.next_setting(star,
                                                   start=rise_time).datetime()
    if start_time < rise_time.datetime():
        start_time = rise_time.datetime()

    transit_time = telescope.observer.next_transit(star,
                                                   start=rise_time).datetime()

    if stop_time > set_time:
        stop_time = set_time
        duration = (stop_time - start_time).total_seconds()

    time_step = float(duration) / float(ntimeslots)
    dtime = np.arange(0, duration + time_step, time_step)
    ha_range = np.zeros((dtime.size), dtype=float)
    for cnt, dt in enumerate(dtime):
        time_ = start_time + timedelta(seconds=dt)
        if time_ < transit_time:
            dt_transit = -(transit_time - time_).total_seconds()
        else:
            dt_transit = (time_ - transit_time).total_seconds()
        ha_range[cnt] = dt_transit / 3600.

    # 1hr = 15 deg = pi/12 rad
    return ha_range * np.pi / 12.

def coverage(telescope,
             target,
             start_time=None,
             duration=300,  # sec
             ntimeslots=10):

    # declination convert in radian
    dec = np.radians(target.dec.degree)
    print(dec)
    # telescope latitude
    latitude = telescope.ref_position.geodetic[1].radian
    print(latitude)
    # number of time slots
    ntimeslots = ntimeslots
    print(ntimeslots)
    # hour angle range from target sky track
    ha_range = get_ha_range(telescope,
                            target,
                            start_time,
                            duration,
                            ntimeslots=ntimeslots)
    print(ha_range)

    # baseline lengths and azimuth angles
    [bl_length,
     bl_az_angle] = telescope.baselines()

    # Plot the uv-Coverage
    uvplot = UVplot(latitude=latitude,
                    declination=dec,
                    elevation=0.)

    fig, ax = uvplot.plot_uv(ha_range,
                             bl_length,
                             bl_az_angle,
                             ntimeslots,
                             comment='{:.0f}'.format(dec))
    uv = uvplot.track_uv(ha_range=ha_range,
                         bl_length=bl_length[-1],
                         bl_azimuth=bl_az_angle[-1],
                         timeslots=ntimeslots)
    mb = 5 * np.sqrt((uv**2).sum(1)).max()
    ax.set_xlim(-mb, mb)
    ax.set_ylim(-mb, mb)
    plt.axis('equal')
    plt.gca().invert_xaxis()

    # Show UV mask and synthesize beam, natural weighting without a taper
    # number of pixels
    npix = 2**10  # should be a power of two
    # kernel matrix
    mask = np.zeros((npix, npix))
    # uv-grid scale factor to fit tracks into matrix
    uvscale = npix / 2 / mb * 0.95 * 0.5
    for i in range(telescope.nr_baselines):
        mask = mask + uvplot.uvMask(ha_range,
                                    bl_length[i],
                                    bl_az_angle[i],
                                    ntimeslots,
                                    npix,
                                    uvscale)

    psf = np.fft.ifftshift(np.fft.ifft2(mask.T)).real
    uvplot.plot_mask(mask.T,
                     psf,
                     comment="natural weighting")<jupyter_output><empty_output><jupyter_text>Select antennas for array<jupyter_code>class antennas():
    def __init__(self):
        
        self.antennas = {}
        self.construct_array()
        
        # build widget grid
        for ant_name in self.antennas.keys():
            self.antennas[ant_name] = self.ant_checkbox(ant_name)
        
    def construct_array(self):
        for cntr in range(0, 64):
            ant_name = 'm{:03d}'.format(cntr)
            self.antennas[ant_name] = None        
        
    def ant_checkbox(self, ant_name):
        return widgets.Checkbox(value=True,
                                description=ant_name,
                                disabled=False,
                                indent=False)<jupyter_output><empty_output><jupyter_text>Target coordinates and observation information<jupyter_code>class RaDec(object):
    def __init__(self):
        radec_header = widgets.Button(description='Enter celestial target coordinates',
                                      disabled=True, button_style='', tooltip='', icon='',
                                      layout={'width':'608px'})        
        header_1_0 = self._add_header_(description='J1939-6342')
        header_1_1 = self._add_header_(description='19h39m25.0264s')        
        header_1_2 = self._add_header_(description='-63d42m45.624s')  
        
        self.cname = self._add_text_(description='Name', value='Test')
        # The right ascension for this object
        self.ra = self._add_text_(description='RA', value='00h42m44.2992s')
        # The declination for this object
        self.dec = self._add_text_(description='Dec', value='+41d16m09.012s')
    
        headers_1 = widgets.HBox([header_1_0, header_1_1, header_1_2])
        celestial_target = widgets.HBox([self.cname, self.ra, self.dec])
        
        header_2 = self._add_header_(description='YYYY-MM-DD HH:mm:SS.f')
        self.starttime = self._add_text_(description='Start time', value='2013-6-21 21:00:00')
        self.duration = self._add_text_(description='Duration [sec]', value='600')
        self.timeslots = self._add_text_(description='# timeslots', value='50')
        start = widgets.VBox([header_2, widgets.HBox([self.starttime, self.duration, self.timeslots])])
        
        self.ui = widgets.VBox([radec_header,
                                headers_1,
                                celestial_target,
                                start
                               ])
    
    def _add_header_(self, description):
        header_layout = {'width':'200px', 'min_width':'200px',
                         'height':'28px', 'min_height':'28px'}
        return widgets.Button(description=description,
                              disabled=True,
                              button_style='',
                              tooltip='',
                              icon='',
                              layout=header_layout)
    
    def _add_text_(self, description, value=''):
        text_style = {'description_width': 'initial'}
        text_layout = {'width':'200px', 'min_width':'200px'}
        return widgets.Text(description=description,
                            value=value,
                            layout=text_layout,
                            style=text_style)<jupyter_output><empty_output><jupyter_text>User input<jupyter_code>def select_antennas(ant_array):
    array_ = np.empty([4, 16], dtype=widgets.Checkbox)
    for cnt, ant in enumerate(ant_array.antennas.keys()):
        array_[cnt//16, cnt%16] = ant_array.antennas[ant]
    antennas_ = widgets.VBox([widgets.HBox(array_[0,:].tolist()),
                              widgets.HBox(array_[1,:].tolist()),
                              widgets.HBox(array_[2,:].tolist()),
                              widgets.HBox(array_[3,:].tolist())])
    display(antennas_)
ant_array = antennas()
select_antennas(ant_array)
target = RaDec()
display(target.ui)<jupyter_output><empty_output><jupyter_text>Run to show UV coverage for selected antennasSelect the appropriate antennas from those listed above and and run the cells generate the coverage plots<jupyter_code>radec_coord = SkyCoord(ra=target.ra.value, dec=target.dec.value,
                       unit=(u.hourangle, u.deg),
                       frame=ICRS)
coverage(telescope=mkat,
         target=radec_coord,
         start_time=target.starttime.value,
         duration=float(target.duration.value),
         ntimeslots=int(target.timeslots.value))<jupyter_output>Test 00h42m44.2992s +41d16m09.012s
<SkyCoord (ICRS): (ra, dec) in deg
    (10.68458, 41.26917)>
2013-6-21 21:00:00 600 50
0.7202828960652683
-0.5360090362073406
50
[-1.85318907 -1.85231641 -1.85144374 -1.85057108 -1.84969841 -1.84882575
 -1.84795308 -1.84708042 -1.84620776 -1.84533509 -1.84446243 -1.84358976
 -1.8427171  -1.84184443 -1.84097177 -1.8400991  -1.83922644 -1.83835377
 -1.83748111 -1.83660844 -1.83573578 -1.83486311 -1.83399045 -1.83311779
 -1.83224512 -1.83137246 -1.83049979 -1.82962713 -1.82875446 -1.8278818
 -1.82700913 -1.82613647 -1.8252638  -1.82439114 -1.82351847 -1.82264581
 -1.82177315 -1.82090048 -1.82002782 -1.81915515 -1.81828249 -1.81740982
 -1.81653716 -1.81566449 -1.81479183 -1.81391916 -1.8130465  -1.81217383
 -1.81130117 -1.81042851 -1.80955584]
{'fraction_non_alphanumeric': 0.11369495720325874, 'fraction_numerical': 0.04439517376508199, 'mean_word_length': 1.6532147742818057, 'pattern_counts': {'":': 0, '<': 27, '<?xml version=': 0, '>': 27, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 25, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '6', 'dir': 'jupyter-structured-clean-dedup', 'id': '453431', 'license': 'BSD-2-Clause', 'n_tokens_mistral': 6624, 'n_tokens_neox': 5935, 'n_words': 1166, 'path': 'planning/MeerKAT Array UV coverage and PSF.ipynb', 'repo_name': 'BongaNjamela001/MeerKAT-Cookbook', 'size': '19394'}