fixed airy beam

Author: zacharymartinot

IonRIME/instrument_setup.py

@@ -62,6 +62,7 @@ def __init__(self, parameters, antenna='a'):
             'paper_feko': 'jones_save/PAPER/',
             'hera_legacy': 'jones_save/',
             'analytic_dipole': 'analytic_dipole/',
+            'airy_dipole': 'jones_save/airy_dipole/',
             'sim_group_test0': 'jones_save/sim_group_test0/'
         }
 
@@ -77,6 +78,7 @@ def __init__(self, parameters, antenna='a'):
             'paper_feko': self.paper_feko,
             'hera_legacy': self.hera_legacy,
             'analytic_dipole': self.analytic_dipole,
+            'airy_dipole': self.airy_dipole,
             'sim_group_test0': self.sim_group_test0,
         }
 
@@ -156,5 +158,9 @@ def hera_legacy(self, parameters):
     def analytic_dipole(self, parameters):
         return irf.analytic_dipole_setup(self.nside, self.nfreq,sigma=self.dipole_hpbw, z0_cza=self.z0_cza)
 
+    def airy_dipole(self, parameters):
+        a = 14.6/2. # HERA dish radius ...or was it 14.0?
+        return irf.airy_dipole_setup(self.nside, self.nu_axis, a, z0_cza=self.z0_cza)
+
     def sim_group_test0(self, parameters):
         return sgtest0.make_ijones_spectrum(parameters)

IonRIME/ionRIME_funcs.py

@@ -3,6 +3,8 @@
 import healpy as hp
 from datetime import datetime
 
+import scipy.special as special
+
 import scipy.interpolate as interpolate
 
 import numba_funcs as irnf
@@ -659,6 +661,77 @@ def transform_basis(nside, jones, z0_cza, R_z0):
 
     return jones_out
 
+def AiryBeam(th, lmbda, a):
+    k = 2. * np.pi / lmbda
+    B = np.abs(2. * special.jv(1, k * a * np.sin(th)) / (k * a * np.sin(th)))
+    B[np.where(np.abs(th) < 1e-10)[0]] = 1
+    B[np.where(th > np.pi/2.)[0]] = 0.
+    return B
+
+def airy_dipole_setup(nside, nu_axis, a, z0_cza=None):
+    def transform_basis(nside, jones, z0_cza, R_z0):
+
+        npix = hp.nside2npix(nside)
+        hpxidx = np.arange(npix)
+        cza, ra = hp.pix2ang(nside, hpxidx)
+
+        fR = R_z0
+
+        tb, pb = rotate_sphr_coords(fR, cza, ra)
+
+        cza_v = t_hat_cart(cza, ra)
+        ra_v = p_hat_cart(cza, ra)
+
+        tb_v = t_hat_cart(tb, pb)
+
+        fRcza_v = np.einsum('ab...,b...->a...', fR, cza_v)
+        fRra_v = np.einsum('ab...,b...->a...', fR, ra_v)
+
+        cosX = np.einsum('a...,a...', fRcza_v, tb_v)
+        sinX = np.einsum('a...,a...', fRra_v, tb_v)
+
+
+        basis_rot = np.array([[cosX, sinX],[-sinX, cosX]])
+        basis_rot = np.transpose(basis_rot,(2,0,1))
+
+        return np.einsum('...ab,...bc->...ac', jones, basis_rot)
+
+    if z0_cza is None:
+        z0_cza = np.radians(120.72)
+
+    nfreq = len(nu_axis)
+
+    npix = hp.nside2npix(nside)
+    hpxidx = np.arange(npix)
+    th, phi = hp.pix2ang(nside, hpxidx)
+
+    R_z0 = hp.rotator.Rotator(rot=[0,-np.degrees(z0_cza)])
+
+    th_l, phi_l = R_z0(th, phi)
+    phi_l[phi_l < 0] += 2. * np.pi
+
+    ct,st = np.cos(th_l), np.sin(th_l)
+    cp,sp = np.cos(phi_l), np.sin(phi_l)
+
+    jones_dipole = np.array([
+            [ct * cp, -sp],
+            [ct * sp, cp]
+        ], dtype=np.complex128).transpose(2,0,1)
+
+    jones_c = transform_basis(nside, jones_dipole, z0_cza, np.array(R_z0.mat))
+
+    c = 299792458.
+    # a = 14.6/2.
+
+    B = np.zeros((nfreq, npix,2,2))
+    for nu_i, nu in enumerate(nu_axis):
+        B[nu_i] = np.broadcast_to(AiryBeam(th_l, c/nu, a), (2,2,npix)).T
+
+    jones_c = np.broadcast_to(jones_c, (nfreq, npix, 2, 2))
+    jones_out = B * jones_c
+
+    return jones_out
+
 def PAPER_transform_basis(nside, jones, z0_cza, R_z0):
 
     npix = hp.nside2npix(nside)