Generate plots comparing basefunctions of OSKAR and EveryBeam

generate_basefunction_plots.py

+import numpy as np
+from matplotlib import pyplot as plt
+from generate_oskar_csv import generate_oskar_csv
+import run_oskar
+from read_oskar_beams import read_oskar_beams
+import subprocess
+import shutil
+from astropy.io import fits
+
+
+l_max = 7
+
+plt.figure(figsize=(10,6))
+
+for em_idx in range(2):
+#for em_idx in [0]:
+
+    for basefunction_idx in range(l_max * (l_max+2)):
+    #for basefunction_idx in [0]:
+        generate_oskar_csv(basefunction_idx,em_idx)
+
+        # run oskar
+
+        run_oskar.main()
+
+        #d1 = fits.getdata('basefunctions_050_MHz_S0000_TIME_SEP_CHAN_SEP_PHASE_XX.fits')[0,0]
+        #d2 = fits.getdata('basefunctions_050_MHz_S0000_TIME_SEP_CHAN_SEP_PHASE_XY.fits')[0,0]
+
+        plt.clf()
+
+        #plt.subplot(1,2,1)
+        #plt.imshow(d1/np.pi*180)
+        #plt.colorbar()
+
+        #plt.subplot(1,2,2)
+        #plt.imshow(d2/np.pi*180)
+        #plt.colorbar()
+
+        #plt.savefig('phases{}'.format(basefunction_idx))
+
+        A = read_oskar_beams()
+        N = A.shape[0]
+
+        print(N)
+
+        v = np.zeros((N,N,4))
+
+        for i in range(N):
+            x = 2.0*i/(N-1) - 1.0
+            for j in range(N):
+                y = 2.0*j/(N-1) - 1.0
+                theta = np.arcsin(np.sqrt(x*x + y*y));
+                phi = np.arctan2(y,x);
+
+                e_theta = np.array([np.cos(phi), np.sin(phi)])
+                e_phi  = np.array([-np.sin(phi), np.cos(phi)])
+
+                T = np.stack((e_theta, e_phi))
+                v[i,j,0] = np.abs(np.dot(A[i,j,:,:], T))[0,0]
+                v[i,j,1] = np.abs(np.dot(A[i,j,:,:], T))[0,1]
+                v[i,j,2] = np.angle(np.dot(A[i,j,:,:],T))[0,0]
+                v[i,j,3] = np.angle(np.dot(A[i,j,:,:],T))[0,1]
+
+        plt.subplot(2,4,1)
+        plt.imshow(v[:,:,0].T, clim=(0,.25), origin='lower')
+        #plt.imshow(v[:,:,0].T, origin='lower')
+        plt.colorbar()
+        plt.title('abs(Etheta)')
+        plt.ylabel('OSKAR')
+
+        plt.subplot(2,4,2)
+        plt.imshow(v[:,:,1].T, clim=(0,.25), origin='lower')
+        plt.colorbar()
+        plt.title('abs(Ephi)')
+
+        plt.subplot(2,4,3)
+        plt.imshow(v[:,:,2].T, clim=(-np.pi, np.pi), cmap='twilight', origin='lower')
+        plt.colorbar()
+        plt.title('angle(Etheta)')
+
+        plt.subplot(2,4,4)
+        plt.imshow(v[:,:,3].T, clim=(-np.pi, np.pi), cmap='twilight', origin='lower')
+        plt.colorbar()
+        plt.title('angle(Ephi)')
+
+        # run EveryBeam
+
+        #generate_oskar_csv(basefunction_idx, em_idx)
+
+        l = int(np.sqrt(basefunction_idx+1))
+        m = basefunction_idx-l*l+1-l
+        s = em_idx
+
+        generate_oskar_csv(l*l-1+l-m, em_idx)
+
+        subprocess.call(["/home/vdtol/src/EveryBeam/oskar/oskar_csv_to_hdf5.py", "telescope.tm", "oskar.h5"])
+
+        shutil.copy("oskar.h5", "/home/vdtol/share/everybeam")
+
+        subprocess.call(['./everybeamtest'])
+
+        A = np.load('response.npy')
+
+        A *= (-1)**(m+1)
+
+        N = A.shape[0]
+
+        print(N)
+
+        v = np.zeros((N,N,4))
+
+        for i in range(N):
+            x = 2.0*i/(N-1) - 1.0
+            for j in range(N):
+                y = 2.0*j/(N-1) - 1.0
+                theta = np.arcsin(np.sqrt(x*x + y*y));
+                phi = np.arctan2(y,x);
+
+                e_theta = np.array([np.cos(phi), np.sin(phi)])
+                e_phi  = np.array([-np.sin(phi), np.cos(phi)])
+
+                T = np.stack((e_theta, e_phi))
+                T = np.eye(2)
+                v[i,j,0] = np.abs(np.dot(A[i,j,:,:], T))[0,0]
+                v[i,j,1] = np.abs(np.dot(A[i,j,:,:], T))[0,1]
+                v[i,j,2] = np.angle(np.dot(A[i,j,:,:],T))[0,0]
+                v[i,j,3] = np.angle(np.dot(A[i,j,:,:],T))[0,1]
+
+        plt.subplot(2,4,5)
+        plt.imshow(v[:,:,0].T, clim=(0,.25), origin='lower')
+        plt.colorbar()
+        plt.title('abs(Etheta)')
+        plt.ylabel('EveryBeam')
+
+        plt.subplot(2,4,6)
+        plt.imshow(v[:,:,1].T, clim=(0,.25), origin='lower')
+        plt.colorbar()
+        plt.title('abs(Ephi)')
+
+        plt.subplot(2,4,7)
+        plt.imshow(v[:,:,2].T, clim=(-np.pi, np.pi), cmap='twilight', origin='lower')
+        plt.colorbar()
+        plt.title('angle(Ex)')
+
+        plt.subplot(2,4,8)
+        plt.imshow(v[:,:,3].T, clim=(-np.pi, np.pi), cmap='twilight', origin='lower')
+        plt.colorbar()
+        plt.title('angle(Ey)')
+
+        #with open('telescope.tm/element_pattern_spherical_wave_x_te_re_0_50.txt') as f:
+            #coeff_line = f.readline()
+            #plt.gcf().suptitle('telescope.tm/element_pattern_spherical_wave_x_te_re_0_50.txt\n{}'.format(coeff_line))
+
+        plt.gcf().suptitle('l = {}, m = {}, s = {}'.format(l,m,s))
+
+        plt.savefig('basefunction{}-{}'.format(basefunction_idx,em_idx))

generate_oskar_csv.py

@@ -9,20 +9,22 @@ import h5py
 from tqdm import tqdm
 from glob import glob
 
-output_dir = 'test.tm'
+output_dir = 'telescope.tm'
 element_id = 0
 
 
 freqs = [50]
 
-l_max = 1
+def generate_oskar_csv(basefunction_idx, em_idx=0):
+
+    l_max = int(np.sqrt(basefunction_idx+1))
 
     # Parse all freqs
     for freq in freqs:
 
         A = np.zeros( (l_max,2*l_max+1, 2,2), dtype=np.complex128)
         data = np.zeros((l_max * (l_max+2), 4), dtype=np.complex128)
-    data[2,0] = 1.0
+        data[basefunction_idx,em_idx] = 1.0
 
         i = 0
         for l in range(l_max):
@@ -37,3 +39,6 @@ for freq in freqs:
                     alpha = np.real(A[:,:,i,j]/s)
 
                     np.savetxt(os.path.join(output_dir,filename), alpha, delimiter=', ')
+
+if __name__ == "__main__":
+    generate_oskar_csv(0)

main.cpp

@@ -13,8 +13,8 @@ int main() {
     std::string name = "station0_LBA";
 
 
-//     auto model = everybeam::ElementResponseModel::OSKARSphericalWave;
-    auto model = everybeam::ElementResponseModel::Hamaker;
+    auto model = everybeam::ElementResponseModel::OSKARSphericalWave;
+//     auto model = everybeam::ElementResponseModel::Hamaker;
 
 
     // Create station.
@@ -29,7 +29,7 @@ int main() {
     double phi = 0.0;
     std::complex<double> response[2][2];
 
-    constexpr int N=100;
+    constexpr int N=256;
 
     std::vector<std::complex<double>> result(N*N*2*2);
 
@@ -46,11 +46,6 @@ int main() {
             double theta = asin(sqrt(x*x + y*y));
             double phi = atan2(y,x);
 
-            double az = M_PI - phi;
-            double el = M_PI_2 - theta;
-
-//             element_response->response(0, freq, theta, phi, result_arr[i][j]);
-
             element_response->response(0, freq, theta, phi, result_arr[i][j]);
         }
     }

makeplots.py

@@ -5,6 +5,32 @@ from matplotlib import pyplot as plt
 
 
 A = np.load('response.npy')
-plt.imshow(abs(A[:,:,0,0])**2 + abs(A[:,:,0,1])**2)
 
-plt.savefig('response.png')
+plt.subplot(2,2,1)
+plt.imshow(abs(A[:,:,0,0])**2, clim=(0.0,np.nanmax(abs(A[:,:,0,0])**2)))
+print (0.0,np.amax(abs(A[:,:,0,0])**2))
+plt.colorbar()
+plt.subplot(2,2,2)
+plt.imshow(abs(A[:,:,0,1])**2, clim=(0.0,np.nanmax(abs(A[:,:,0,1])**2)))
+plt.colorbar()
+plt.subplot(2,2,3)
+plt.imshow(abs(A[:,:,1,0])**2, clim=(0.0,np.nanmax(abs(A[:,:,1,0])**2)))
+plt.colorbar()
+plt.subplot(2,2,4)
+plt.imshow(abs(A[:,:,1,1])**2, clim=(0.0,np.nanmax(abs(A[:,:,1,1])**2)))
+plt.colorbar()
+
+plt.savefig('response_amp.png')
+
+plt.figure()
+
+plt.subplot(2,2,1)
+plt.imshow(np.angle(A[:,:,0,0])**2)
+plt.subplot(2,2,2)
+plt.imshow(np.angle(A[:,:,0,1])**2)
+plt.subplot(2,2,3)
+plt.imshow(np.angle(A[:,:,1,0])**2)
+plt.subplot(2,2,4)
+plt.imshow(np.angle(A[:,:,1,1])**2)
+
+plt.savefig('response_phase.png')

oskar_csv_to_hdf5.py

+#!/usr/bin/env python3
+
+import argparse
+import os
+import scipy.io as sio
+import os.path
+import numpy as np
+import h5py
+from tqdm import tqdm
+from glob import glob
+import csv
+
+# Parse command-line arguments
+parser = argparse.ArgumentParser()
+parser.add_argument("input_dir", type=str, help="directory containing coefficients")
+parser.add_argument("output_file", type=str, help="name of the HDF5 output file")
+parser.add_argument("-debug", dest="debug", action="store_true")
+parser.set_defaults(input_dir=".")
+parser.set_defaults(output_file="oskar.h5")
+parser.set_defaults(debug=False)
+args = parser.parse_args()
+
+# Read the input directory
+# This directory should contain subdirectories,
+# with names corresponding to frequency.
+input_dir = args.input_dir
+print("Reading coefficients from: " , input_dir)
+
+# Only read element_id 0, because EveryBeam's implementation of the OSKAR model
+# currently supports only a single model for all elements
+element_id = 0
+nr_elements = 1
+
+files = glob(os.path.join(input_dir, 'element_pattern_spherical_wave_?_t?_??_{}_*.txt').format(element_id))
+
+files = [os.path.basename(f) for f in files]
+
+freqs = sorted(set([int(f[:-4].split('_')[-1]) for f in files]))
+
+
+print("Parsing frequencies %d-%d MHz" % (int(freqs[0]), int(freqs[-1])))
+
+# Create HDF5 file
+output_file = args.output_file
+print("Creating HDF5 file: ", output_file)
+hf = h5py.File(output_file, 'w')
+
+# Debugging
+debug = args.debug
+if (debug):
+    freqs = [freqs[0]]
+    tqdm = lambda x: x
+
+# Parse all freqs
+for freq in tqdm(freqs):
+
+    A = None
+
+    for i, pol in enumerate(('x', 'y')):
+        for j, em in enumerate(('e', 'm')):
+            for reim, s in zip(('re', 'im'), (1.0, 1.0j)) :
+                filename = 'element_pattern_spherical_wave_{}_t{}_{}_{}_{}.txt'.format(pol, em, reim, element_id, freq)
+                alpha = np.loadtxt(os.path.join(input_dir,filename), delimiter=',', ndmin=2)
+                l_max = alpha.shape[0]
+                assert(alpha.shape[1] == 2*l_max + 1)
+                if A is None:
+                    A = np.zeros( (l_max,2*l_max+1, 2,2), dtype=np.complex128)
+                else:
+                    assert(A.shape == (l_max,2*l_max+1, 2,2))
+                A[:,:,i,j] += alpha*s
+
+    # Create the output matrix
+    data = np.zeros((nr_elements, l_max * (l_max+2), 4), dtype=np.complex128)
+
+    # Fill the output matrix
+    # The inner dimension are set as follows:
+    # (x_te_re, x_te_im), (x_tm_re, x_tm_im),
+    # (y_te_re, y_te_im), (y_tm_re, y_tm_im)
+
+    i = 0
+    for l in range(l_max):
+        n = (l+1) * 2 + 1
+        data[0,i:i+n,:] = A[l,:n,:].reshape(n,4)
+        i += n
+
+    # Add group for current frequency
+    hf.create_dataset(str(freq), data.shape, data=data)
+
+# Close HDF5 file
+hf.close()
+
+# Done
+print("Done!")

plot_phases.py

+from astropy.io import fits
+from matplotlib import pyplot as plt
+
+d1 = fits.getdata('basefunctions_050_MHz_S0000_TIME_SEP_CHAN_SEP_PHASE_XX.fits')[0,0]
+d2 = fits.getdata('basefunctions_050_MHz_S0000_TIME_SEP_CHAN_SEP_PHASE_XY.fits')[0,0]
+
+plt.subplot(1,2,1)
+plt.imshow(d1/np.pi*180)
+colorbar()
+
+plt.subplot(1,2,2)
+plt.imshow(d2/np.pi*180)
+colorbar()
+
+plt.show()

read_oskar_beams.py

@@ -4,24 +4,31 @@ from astropy.io import fits
 
 import numpy as np
 
-
+def read_oskar_beams():
     A = None
-
     for i, pol1 in enumerate(['X', 'Y']):
         for j, pol2 in enumerate(['X', 'Y']):
             filename_amp = 'basefunctions_050_MHz_S0000_TIME_SEP_CHAN_SEP_AMP_{}{}.fits'.format(pol1, pol2)
             filename_phase = 'basefunctions_050_MHz_S0000_TIME_SEP_CHAN_SEP_PHASE_{}{}.fits'.format(pol1, pol2)
-        d = (fits.getdata(filename_amp) * np.exp(1j*fits.getdata(filename_phase)))[0,0,:,:]
+            d = (fits.getdata(filename_amp) * np.exp(1j*fits.getdata(filename_phase)))[0,0,:,:].T
             if A is None:
                 N = d.shape[-1]
                 A = np.zeros((N,N,2,2), dtype=np.complex128)
             A[:,:,i,j] = d
 
+    return A
+
+
 
+if __name__ == "__main__":
+
+
+    A = read_oskar_beams()
+    N = A.shape[0]
 
     print(N)
 
-v = np.zeros((N,N))
+    v = np.zeros((N,N,4))
 
     for i in range(N):
         x = 2.0*i/(N-1) - 1.0
@@ -29,13 +36,38 @@ for i in range(N):
             y = 2.0*j/(N-1) - 1.0
             theta = np.arcsin(np.sqrt(x*x + y*y));
             phi = np.arctan2(y,x);
-        e_r = np.array([x,y, np.cos(theta)])
-        e_phi = np.array([np.cos(theta)*np.cos(phi), np.cos(theta)*np.sin(phi), -np.sin(theta)])
-        e_theta = np.array([-y, x, 0.0])
-        T = np.stack((e_phi, e_theta)).T[:2,:]
-        v[i,j] = abs(np.dot(A[i,j,:,:], T))[0,1]
-        U,S,V = np.linalg.svd(A[i,j,:,:])
-        v[i,j] = np.angle(A[i,j,0,1])
+
+            e_theta = np.array([np.cos(phi), np.sin(phi)])
+            e_phi  = np.array([-np.sin(phi), np.cos(phi)])
+
+            T = np.stack((e_theta, e_phi)).T
+            v[i,j,0] = np.abs(np.dot(A[i,j,:,:], T))[0,0]
+            v[i,j,1] = np.abs(np.dot(A[i,j,:,:], T))[0,1]
+            v[i,j,2] = np.mod(np.angle(A[i,j,0,0]),2*np.pi)
+            v[i,j,3] = np.angle(A[i,j,0,1])
+
+    plt.subplot(2,2,1)
+    plt.imshow(v[:,:,0], clim=(0,.25))
+    plt.colorbar()
+    plt.title('abs(Etheta)')
+
+    plt.subplot(2,2,2)
+    plt.imshow(v[:,:,1], clim=(0,.25))
+    plt.colorbar()
+    plt.title('abs(Ephi)')
+
+    plt.subplot(2,2,3)
+    plt.imshow(v[:,:,2])
+    plt.colorbar()
+    plt.title('angle(Ex)')
+
+    plt.subplot(2,2,4)
+    plt.imshow(v[:,:,3])
+    plt.colorbar()
+    plt.title('angle(Ey)')
+
+    plt.show()
+
 
 
 

run_oskar.py

@@ -99,12 +99,13 @@ def main():
         'observation': {
             'phase_centre_ra_deg': ra0_deg,
             'phase_centre_dec_deg': dec0_deg,
+            #'pointing_file': 'station_pointing.txt',
             'start_time_utc': get_start_time(ra0_deg, length_sec),
             'length': length_sec
         },
         'telescope': {
             'normalise_beams_at_phase_centre': False,
-            'aperture_array/element_pattern/normalise': True
+            'aperture_array/element_pattern/normalise': False
         },
         'beam_pattern': {
             'coordinate_frame': 'Horizon',
@@ -122,7 +123,7 @@ def main():
     # Define telescope models to use, and associated overrides for them.
     telescopes = {
         'basefunctions': {
-            'telescope/input_directory': 'test.tm',
+            'telescope/input_directory': 'telescope.tm',
             'telescope/aperture_array/element_pattern/enable_numerical': True,
             'telescope/aperture_array/element_pattern/swap_xy': True,
             'telescope/aperture_array/array_pattern/enable': False
@@ -156,13 +157,13 @@ def main():
             subprocess.call([app, settings_path])
 
             # Make plots.
-            title = tel + ' @ ' + str(freq) + ' MHz'
-            make_plots(title=title,
-                       glob_pattern=root_path+'*_AMP*',
-                       out_basename=root_path+'_amp')
-            make_plots(title=title,
-                       glob_pattern=root_path+'*_AUTO_POWER*',
-                       out_basename=root_path+'_stokes')
+            #title = tel + ' @ ' + str(freq) + ' MHz'
+            #make_plots(title=title,
+                       #glob_pattern=root_path+'*_AMP*',
+                       #out_basename=root_path+'_amp')
+            #make_plots(title=title,
+                       #glob_pattern=root_path+'*_AUTO_POWER*',
+                       #out_basename=root_path+'_stokes')
 
 
 if __name__ == '__main__':