Some proposed visualizations of Cross-coupling beam models

Jacobs and B. Gehlot
14 Dec 2020 - First draft
2 Feb 2021 - Add more description and make cross power plots same as auto power.

Introduction

Cross-coupling is expected to play a role in HERA’s systematic response. Nic Fagnoni generated a simulation with several antennas to study the impact of cross coupling in the HERA Phase 1 feeds, the PAPER cross-dipoles surrounded by a "tuna can" reflector (Fagnoni et al 2020). Nic simulated a 19 element Hex array (map shown below). From this simulation we have beams for three antennas indicated by the green boxes. In this memo we take a closer look at these simulation results from the perspective identifying a physical intuition about cross-coupling and connecting insights to observations.

Assumptions and Predictions

A simplistic physical model for antenna to antenna cross coupling is the addition of multiple paths caused by reflections. This will have the effect of adding copies of the signal each delayed according to the length of the return path. A beam model encodes these simulations in the beam model. In the HERA delay spectrum context reflections emerge additional direction dependent response which has an additional delay/phase slope corresponding to the multipath length and an amplitude which is net of the gain at each multipath reflection and the path loss. Beam maps can not represent multi-path lengths longer than c/(channel width). Auto correlations will include reflections undergoing a round trip out and back. Since two reflections are required, we might expect them to be lower by 2x in dB.

Execution note

Execute this notebook in the HERA-Beams repo synced to the trunk branch. This repo uses git lfs. Make sure to git lfs checkout both NF_CrossCouplingBeams and NicolasFagnoniBeams to force the files to download.

In [1]:
from IPython.core.display import display, HTML
display(HTML("<style>.container { width:90% !important; }</style>"))
from IPython.display import clear_output
%matplotlib inline
In [2]:
import numpy as np
from time import sleep
from copy import deepcopy
import matplotlib.pyplot as plt
from pyuvdata import UVData, UVBeam
from matplotlib.colors import LogNorm
from matplotlib.animation import FuncAnimation, PillowWriter

from numpy.fft import fftshift, fftfreq, fft
from scipy import signal
from glob import glob
import os
In [3]:
import matplotlib as mpl
mpl.rcParams['xtick.major.size']  = 6
mpl.rcParams['ytick.major.size']  = 6
mpl.rcParams['xtick.minor.size']  = 3
mpl.rcParams['ytick.minor.size']  = 3
mpl.rcParams['xtick.major.width'] = 1.0
mpl.rcParams['ytick.major.width'] = 1.0
mpl.rcParams['xtick.minor.width'] = 1.0
mpl.rcParams['ytick.minor.width'] = 1.0
mpl.rcParams['axes.linewidth'] = 1.0
mpl.rcParams['font.family'] = 'serif'
mpl.rcParams['figure.facecolor'] = 'white'
In [4]:
def colorbar_plotter(fig,ax,im,label):
        
    pos = ax.get_position()
    cbarax = fig.add_axes([pos.x0 + pos.width+0.01, pos.y0, 0.008, pos.height])

    cbar = fig.colorbar(im,cax=cbarax)

    cbar.ax.tick_params(labelsize=16) 
    cbar.ax.set_ylabel(label,fontsize=16)
    
    return cbar

def get_normalized_power_beam(data):
    power = np.real_if_close(data[0,0]*np.conj(data[0,0]) + data[1,0]*np.conj(data[1,0]))
    norm = np.max(power,axis=(1,2,3))
    
    return power/norm[:,np.newaxis,np.newaxis,np.newaxis]

def cross_power_beam(beam1, beam2):
    power = beam1[0,0]*np.conj(beam2[0,0]) + beam1[1,0]*np.conj(beam2[1,0])
    norm = np.max(np.abs(power),axis=(1,2,3))
    return power/norm[:,np.newaxis,np.newaxis,np.newaxis]
        
def dB(x):
    return 10*np.log10(x)

def FFT(d,dt=0.5):
    x = d*(signal.blackmanharris(d.shape[1]))[np.newaxis,:,np.newaxis,np.newaxis]
    return fftshift(fftfreq(d.shape[1],d=dt)),fftshift(fft(x,axis=1),axes=1)

def nearest_idx(data,val):
    return np.abs(data - val).argmin()

def load_beamfitses(files):
    # input a list of beam fites, 
    # output a single concatenated beam file
    #   ASSUMES you have a single sensible axis to append along
    b = UVBeam()
    b.read_beamfits(files[0])
    if len(files)==1:return b
    for f in files[1:]:
        L = UVBeam()
        L.read_beamfits(f)
        b += L
        del(L)
    return deepcopy(b)
def load_beampols(files):
    # input a dict {pol:filename}, where pol is E or N 
    # output a single concatenated beam object
    # 
    E = UVBeam()
    E.read_beamfits(files['E'])
    E.select(feeds=['x'],inplace=True)
    N = UVBeam()
    N.read_beamfits(files['N'])
    N.select(feeds=['y'],inplace=True)
    return E+N

Load Data

Two data sets. Isolated beam which we load as "old" and cross coupling beams. Here we load the ports corresponding to the EE polarization.

In [5]:
beamfile_old = ['HERA-Beams/NicolasFagnoniBeams/NF_HERA_Dipole_efield_beam.fits']

beamfiles_0506  = {'E':'HERA-Beams/NF_CrossCouplingBeams/NF_HERA-Dipole_ccbeam_port05_E-field.fits',
                   'N':'HERA-Beams/NF_CrossCouplingBeams/NF_HERA-Dipole_ccbeam_port06_E-field.fits'}
beamfiles_2122 = {'E':'HERA-Beams/NF_CrossCouplingBeams/NF_HERA-Dipole_ccbeam_port21_E-field.fits',
                  'N':'HERA-Beams/NF_CrossCouplingBeams/NF_HERA-Dipole_ccbeam_port22_E-field.fits'}
beamfiles_2526 = {'E':'HERA-Beams/NF_CrossCouplingBeams/NF_HERA-Dipole_ccbeam_port25_E-field.fits',
                  'N':'HERA-Beams/NF_CrossCouplingBeams/NF_HERA-Dipole_ccbeam_port26_E-field.fits'}
# beamfile_0506 = 'HERA_ccbeam_port5-6.uvbeam'
# beamfile_2122 = 'HERA_ccbeam_port21-22.uvbeam'
# beamfile_2526 = 'HERA_ccbeam_port25-26.uvbeam'

beam_old = load_beamfitses(beamfile_old)
beam_0506 = load_beampols(beamfiles_0506)
beam_2122 = load_beampols(beamfiles_2122)
beam_2526 = load_beampols(beamfiles_2526)

# beam_old.read_beamfits(beamfile_old)
# beam_0506.read_beamfits(beamfile_0506)
# beam_2122.read_beamfits(beamfile_2122)
# beam_2526.read_beamfits(beamfile_2526)

Aligning

The two datasets were computed on different frequency grids. The biggest difference is that the cross coupling beam is at a higher resolution. We interpolate the isolated beam down to that denser sampling. We'll need to remember not to draw any conclusions beyond the original nyquist limit.

In [6]:
beam_old.freq_interp_kind = 'linear'
beam_old.interpolation_function = 'az_za_simple'
beam_old_interp = beam_old.interp(freq_array=beam_0506.freq_array[0],new_object=True)
In [8]:
data_old = beam_old_interp.data_array

data_0506 = beam_0506.data_array
data_2122 = beam_2122.data_array
data_2526 = beam_2526.data_array

Power beams

We form up the power beams by squaring the electric field vector components and normalizing by the maximum taken along theta, phi, and frequency but allowing a normalization for each polarization.

Here we use a custom power beam which arose as we developed a function to make cross-power beams. This could be added to uvbeam if theres more interest.

In [9]:
power_old = get_normalized_power_beam(data_old)
power_0506 = get_normalized_power_beam(data_0506)
power_2122 = get_normalized_power_beam(data_2122)
power_2526 = get_normalized_power_beam(data_2526)

Auto Beams

The autocorrelation beams are the power beams at 150MHz (channel # 200).

In [45]:
power_beams = [power_old,power_0506,power_2122,power_2526]
titles = ['Old_beam', 'Port 5-6', 'Port 21-22', 'Port 25-26']
#beam_models = [beam_old, ]
theta = beam_0506.axis2_array
phi = beam_0506.axis1_array

AZ, ZA = np.meshgrid(phi,theta)

fig1, axs = plt.subplots(nrows=2,ncols=4,figsize=(24,12),subplot_kw=dict(projection='polar'))
fig1.subplots_adjust(wspace=0.2,hspace=0.1)

for i in range(4):
    for j in range(2):
        im = axs[j,i].scatter(AZ.ravel(),ZA.ravel(),c=dB(power_beams[i][j,200,:,:]).ravel(),vmin=-70,vmax=0.)
        axs[j,i].set_ylim(0,np.pi/2.)
        axs[j,i].tick_params(axis='both',direction='in',which='both',labelsize=14)
        if i==3:
            colorbar_plotter(fig1,axs[j,i],im,r'')        
        if j==0:
            axs[j,i].set_title(titles[i],y=1.1,fontsize=16)

The auto correlation beams left to right: Isolated beam, South antenna, Middle antenna, Eastern antenna. Top row is EE and bottom NN. Azimuthal distribution of ripples is directionally consistent with reflections in the nearest hex positions. These are fringes.

Lets compare directly to the isolated antenna with a difference.

In [52]:
beam_diff = [power_0506 - power_old, power_2122 - power_old, power_2526 - power_old]
In [53]:
fig2, axs = plt.subplots(nrows=2,ncols=3,figsize=(24,16),subplot_kw=dict(projection='polar'))

for i in range(3):
    for j in range(2):       
        
        im = axs[j,i].scatter(AZ.ravel(),ZA.ravel(),c=beam_diff[i][j,200,:,:].ravel(),vmin=-0.002,vmax=0.002,cmap='bwr')

        axs[j,i].set_ylim(0,np.pi/2.)
        if i==2:
            colorbar_plotter(fig2,axs[j,i],im,r'')
        axs[j,i].tick_params(axis='both',direction='in',which='both',labelsize=14)
        if j==0:    
            axs[j,i].set_title('Port 5-6',y=1.1,fontsize=16)

Auto correlation beams of cross coupling sim, isolated beam subtracted. Looks like a variation of 0.001 in linear units on top of ~-10 to -20 dB beam level. So thats 1 to 10%.

Delay Response of Auto Beam

Sources fringing off of a nearby reflector will have a delay corresponding to their geometric path length plus the delay from the reflector. The fringe amplitude should average down with frequency causing the zero delay / spectral average to be smaller amplitude than at a single frequency.

In [12]:
freqs = beam_0506.freq_array[0,100:301]
df = freqs[1] - freqs[0]
In [13]:
F_ccbeams = []
for i in range(4):
    delays, FPbeam = FFT(power_beams[i][:,100:301],dt=df)
    FPbeam_norm  = np.abs(FPbeam)
    FPbeam_norm /= np.max(FPbeam_norm,axis=(1,2,3))[:,np.newaxis,np.newaxis,np.newaxis] #normalize by max across delays, az, za for ea pol
    F_ccbeams.append(FPbeam_norm)
In [27]:
# fig3, axs = plt.subplots(nrows=2,ncols=4,figsize=(24,12),subplot_kw=dict(projection='polar'))

# for i in range(2): #pol
#     for j in range(4): #which antenna

#         im = axs[i,j].scatter(AZ.ravel(),ZA.ravel(),c=dB(F_ccbeams[j][i,100,:,:].ravel()),vmin=-50,vmax=0)  
        
#         if j==3:
#             colorbar_plotter(fig3,axs[i,j],im,r'')        
        
#         axs[i,j].set_ylim(0,0.5*np.pi)
#         axs[i,j].tick_params(axis='both',direction='in',which='both',labelsize=14)
#         if i==0: 
#             axs[i,j].set_title(titles[j] + r' and $\tau = 0$',y=1.05,fontsize=18)
In [35]:
# Make more maps at useful delay channels
# plot delay spectra pierce points like the figure below
# the cross power beam in cross coupling model, delay space is: F_beams[j][0,100+idx[i],:,:]  #bl,pol, channel,theta,phi
# the auto power beam in cross coupling model, delay space is: F_ccbeams (same indexing, except bl becomes ant)
# sometimes we see A,B,C naming. A is south, B is center, C is East.
tau_cut = delays[100:]
idx = [0, np.abs(tau_cut - 1.e-7).argmin(), np.abs(tau_cut - 2.e-7).argmin(), np.abs(tau_cut - 4.e-7).argmin()]
print(idx)
print(tau_cut[idx[-1]]*1e9)

[0, 10, 20, 40]
398.0099502487562
In [43]:
tau_cut = delays[100:]
idx = [0, np.abs(tau_cut - 1.e-7).argmin(), np.abs(tau_cut - 2.5e-7).argmin(), np.abs(tau_cut - 4.e-7).argmin()]

fig, axs = plt.subplots(nrows=8,ncols=4,figsize=(24,48),subplot_kw=dict(projection='polar'))

for i in range(4): #delay bins
    for j in range(4):  #beams

        im1 = axs[2*i,j].scatter(AZ.ravel(),ZA.ravel(),c=dB(F_ccbeams[j][0,100+idx[i],:,:].ravel()))  #E
        im2 = axs[2*i+1,j].scatter(AZ.ravel(),ZA.ravel(),c=dB(F_ccbeams[j][1,100+idx[i],:,:].ravel())) #N
        
        if j==3:
            colorbar_plotter(fig,axs[2*i,j],im1,r'')        
        
        axs[2*i,j].set_ylim(0,0.5*np.pi)
        axs[2*i+1,j].set_ylim(0,0.5*np.pi)             
        axs[2*i,j].tick_params(axis='both',direction='in',which='both',labelsize=14)
        axs[2*i+1,j].tick_params(axis='both',direction='in',which='both',labelsize=14)
    
        if j==1:
            print(r'$\tau =$' + '%.2f ns'%(tau_cut[idx[i]]*1.e9))
            axs[2*i,j].set_title(r'$\tau =$' + '%.2f ns'%(tau_cut[idx[i]]*1.e9),y=1.05,fontsize=18)
        if j==0:
            axs[2*i,j].set_ylabel('E',labelpad=50)
            axs[2*i+1,j].set_ylabel('N',labelpad=50)

$\tau =$0.00 ns
$\tau =$99.50 ns
$\tau =$248.76 ns
$\tau =$398.01 ns

Delay response of cross-coupling beam

The cross products of the beams tell a more complete story. Lets look at the same plots as before but now $E_i * E_j^\dag$. First the delay spectra along the same peirce points and then beams at interesting delay slices.

In [24]:
colors = ['r','b','k']
labels = [r'$\phi_{az} = 0^{\circ}$',r'$\phi_{az} = 45^{\circ}$',r'$\phi_{az} = 90^{\circ}$']


arrows = [350.,-30]
dx = [0.,100.,100.]
dy = [10,0,10]

phi_in = np.array([90.,0.,45.])*np.pi/180. + np.pi #az angles requested
theta_in = np.array([30.,45.,60.])*np.pi/180. #el angles requested
phi_rot = np.array([0.,45.,90.])*np.pi/180. #baseline orientations
titles = ['Isolated', 'South',  'Center', 'East']
pols = ['E','N']

fig6, axs = plt.subplots(nrows=2,ncols=3,figsize=(24,12))

for i in range(2): # pol
    for j in range(3): #which antenna (was baseline)
        for k in range(3): #angle
                
            #phi_rel = phi_in[j] + phi_rot[k]
            phi_rel = phi_rot[k]
            if phi_rel == 2*np.pi:
                phi_rel = 0.
                    
            axs[i,j].plot(delays*1.e9,dB(F_ccbeams[j+1][i,:,nearest_idx(theta,theta_in[1]),nearest_idx(phi,phi_rel)]),lw=2,c=colors[k],label=labels[k])
        
            if phi_rot[k]==0.:
                axs[i,j].plot(delays*1.e9,dB(F_ccbeams[0][i,:,nearest_idx(theta,theta_in[1]),nearest_idx(phi,phi_rel)]),lw=2,c=colors[k],ls='--',
                              label=r'Old Beam ($\phi_{rel} = 0^{\circ}$)')
        
        axs[i,j].set_ylim(-90.,-20.)
        axs[i,j].set_xlim(-500.,500.)
        axs[i,j].tick_params(axis='both',direction='in',which='both',labelsize=16)
        
        
        #axs[i,j].arrow(arrows[0],-50.,dx[j],dy[j],lw=3,ls='--',color='grey')
        axs[i,j].set_title(titles[j+1]+' '+pols[i])
        
        if i==1:
            axs[i,j].set_xlabel(r'$\tau\,[$ns$]$',fontsize=16)
        else:
            if j==0:
                handles, labels = axs[i,j].get_legend_handles_labels()
                axs[i,j].legend(handles,labels,fontsize=16,ncol=4,loc='upper left',bbox_to_anchor=(0.8,1.2))
    
        if j==0:
            axs[i,j].set_ylabel('[dB]',fontsize=16)
In [21]:
beams_bl = [cross_power_beam(data_0506,data_2122),
            cross_power_beam(data_2122,data_2526),
            cross_power_beam(data_0506,data_2526)]
In [22]:
F_beams = []
for i in range(3):
    delays, Fbeam = FFT(beams_bl[i][:,100:301],dt=df)
    Fbeam_norm  = np.abs(Fbeam)
    Fbeam_norm /= np.max(Fbeam_norm,axis=(1,2,3))[:,np.newaxis,np.newaxis,np.newaxis]
    F_beams.append(Fbeam_norm)

Cross-power beams at interesting delays

Looking at the same delays from the autos. Now might expect some potential for asymmetry between positive and negative delays, so we add some negative bins.

Caption data: Columns are baselines: South, E, SE Rows are alternating EE and NN pols (todo indicate pol on figure)

In [23]:
tau_cut = delays#[100:]
idx = [0, np.abs(tau_cut - 1.e-7).argmin(), 
       np.abs(tau_cut - 2.e-7).argmin(), 
       np.abs(tau_cut - 4.e-7).argmin(), 
       np.abs(tau_cut + 4.e-7).argmin()]


fig5, axs = plt.subplots(nrows=6,ncols=3,figsize=(24,48),subplot_kw=dict(projection='polar'))

for i in range(3):
    for j in range(3):

        im1 = axs[2*i,j].scatter(AZ.ravel(),ZA.ravel(),c=dB(F_beams[j][0,idx[i],:,:].ravel()))  
        im2 = axs[2*i+1,j].scatter(AZ.ravel(),ZA.ravel(),c=dB(F_beams[j][1,idx[i],:,:].ravel()))
        
        if j==2:
            colorbar_plotter(fig5,axs[2*i,j],im1,r'')        
        
        axs[2*i,j].set_ylim(0,0.5*np.pi)
        axs[2*i+1,j].set_ylim(0,0.5*np.pi)             
        axs[2*i,j].tick_params(axis='both',direction='in',which='both',labelsize=14)
        axs[2*i+1,j].tick_params(axis='both',direction='in',which='both',labelsize=14)
    
        if j==1:
            axs[2*i,j].set_title(r'$\tau =$' + '%.2f ns'%(tau_cut[idx[i]]*1.e9),y=1.05,fontsize=18)
In [72]:
phi_in = np.array([90.,0.,45.])*np.pi/180. + np.pi
theta_in = np.array([30.,45.,60.])*np.pi/180.
phi_rot = np.array([0.,45.,90.])*np.pi/180.
In [74]:
colors = ['r','b','k']
labels = [r'$\phi_{rel} = 0^{\circ}$',r'$\phi_{rel} = 45^{\circ}$',r'$\phi_{rel} = 90^{\circ}$']

fig6, axs = plt.subplots(nrows=2,ncols=3,figsize=(24,12))

for i in range(2):
    for j in range(3):
        for k in range(3):
                
            phi_rel = phi_in[j] + phi_rot[k]
            if phi_rel == 2*np.pi:
                phi_rel = 0.
                    
            axs[i,j].plot(delays*1.e9,dB(F_beams[j][i,:,nearest_idx(theta,theta_in[1]),nearest_idx(phi,phi_rel)]),lw=2,c=colors[k],label=labels[k])
        
            if phi_rot[k]==0.:
                axs[i,j].plot(delays*1.e9,dB(F_ccbeams[0][i,:,nearest_idx(theta,theta_in[1]),nearest_idx(phi,phi_rel)]),lw=2,c=colors[k],ls='--',
                              label=r'Old Beam ($\phi_{rel} = 0^{\circ}$)')
        
        axs[i,j].set_ylim(-90.,-20.)
        axs[i,j].set_xlim(-500.,500.)
        axs[i,j].tick_params(axis='both',direction='in',which='both',labelsize=16)
        
        axs[i,j].arrow(arrows[0],-50.,dx[j],dy[j],lw=3,ls='--',color='grey')
        
        if i==1:
            axs[i,j].set_xlabel(r'$\tau\,[$ns$]$',fontsize=16)
        else:
            if j==0:
                handles, labels = axs[i,j].get_legend_handles_labels()
                axs[i,j].legend(handles,labels,fontsize=16,ncol=4,loc='upper left',bbox_to_anchor=(0.8,1.2))
    
        if j==0:
            axs[i,j].set_ylabel('[dB]',fontsize=16)

Delay spectra through several pierce points. Note how the delay response peaks not at zero. This happens because CST simulates the phase offset due to the antenna position.

Animation vs delay

This here bit of code renders an animation of the beam vs delay. Play with the subset selection in the below cell and the colorbar range in lines 24 and 25 of the update function in the block after that.

In [42]:
imdata = dB(F_beams[1])#[:,50:151,:,:])
taus = delays#[50:151]
In [43]:
fig_anim, axs = plt.subplots(nrows=1,ncols=2,figsize=(16,8),subplot_kw=dict(projection='polar'))  
xd = AZ.ravel()
yd = ZA.ravel()

scat1 = axs[0].scatter(xd,yd,c=[],cmap='viridis',vmin=-80.,vmax=0.)
scat2 = axs[1].scatter(xd,yd,c=[],cmap='viridis',vmin=-80.,vmax=0.)

def init():  
    axs[0].set_ylim(0, 0.5*np.pi)
    axs[1].set_ylim(0, 0.5*np.pi)
    axs[0].tick_params(axis='both',direction='in',which='both',labelsize=14)
    axs[1].tick_params(axis='both',direction='in',which='both',labelsize=14)
    
    
def update(i):
    label = r'$\tau =$' + '%.2f ns'%(taus[i]*1.e9)
#    print(label)
    # Update the line and the axes (with a new xlabel). Return a tuple of
    # "artists" that have to be redrawn for this frame.
    scat1.set_color(scat1.to_rgba(imdata[0,i,:,:].ravel()))
    scat2.set_color(scat2.to_rgba(imdata[1,i,:,:].ravel()))
    axs[0].set_title(label,y = 1.1, fontsize=18)
    axs[1].set_title(label,y = 1.1, fontsize=18)
    scat1.set_clim(vmin=-80.,vmax=0.)
    scat2.set_clim(vmin=-80.,vmax=0.)
    
    return scat1, scat2, axs

anim = FuncAnimation(fig_anim, update, frames=len(taus), interval=100,init_func=init)
#anim.save('ccbeam_AB_dual-pol.gif', dpi=80, writer='imagemagick')
anim.save('ccbeam_East_baseline_dual-pol_alldelays.mov', dpi=80, writer='ffmpeg')
print("done!")

done!
In [ ]: