diff --git a/generate_basefunction_plots.py b/generate_basefunction_plots.py
new file mode 100644
index 0000000000000000000000000000000000000000..bd2fdb93bcbb21d90c2b92e4164a437bbf35e4de
--- /dev/null
+++ b/generate_basefunction_plots.py
@@ -0,0 +1,156 @@
+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))
diff --git a/generate_oskar_csv.py b/generate_oskar_csv.py
index 38ce94ab35dbe61d9ee6f2c4e429400c3d2e72a7..91d2d42195a9efc5245a1232cf49fbed2725a94f 100755
--- a/generate_oskar_csv.py
+++ b/generate_oskar_csv.py
@@ -9,31 +9,36 @@ 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):
-# Parse all freqs
-for freq in freqs:
+ l_max = int(np.sqrt(basefunction_idx+1))
- 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
+ # Parse all freqs
+ for freq in freqs:
- i = 0
- for l in range(l_max):
- n = (l+1) * 2 + 1
- A[l,:n,:] = data[i:i+n,:].reshape(n,2,2)
- i += n
+ 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[basefunction_idx,em_idx] = 1.0
- 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.real(A[:,:,i,j]/s)
+ i = 0
+ for l in range(l_max):
+ n = (l+1) * 2 + 1
+ A[l,:n,:] = data[i:i+n,:].reshape(n,2,2)
+ i += n
- np.savetxt(os.path.join(output_dir,filename), alpha, delimiter=',')
+ 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.real(A[:,:,i,j]/s)
+
+ np.savetxt(os.path.join(output_dir,filename), alpha, delimiter=', ')
+
+if __name__ == "__main__":
+ generate_oskar_csv(0)
diff --git a/main.cpp b/main.cpp
index e89932572e6f17579cfd8b5781f9fdb13542c3ef..a81ca0cc0564042075ad376c2f708fed024c7fac 100644
--- a/main.cpp
+++ b/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]);
}
}
diff --git a/makeplots.py b/makeplots.py
index 08051f0daaca50984daf7d4bbdd23c8e2f30debc..8685370ac51b1dbe09bf1582f48f3db577faba32 100755
--- a/makeplots.py
+++ b/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')
diff --git a/oskar_csv_to_hdf5.py b/oskar_csv_to_hdf5.py
new file mode 100644
index 0000000000000000000000000000000000000000..becae4df60c5d4ab491cc6d7cffadf0ae33bdee4
--- /dev/null
+++ b/oskar_csv_to_hdf5.py
@@ -0,0 +1,93 @@
+#!/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!")
diff --git a/plot_phases.py b/plot_phases.py
new file mode 100644
index 0000000000000000000000000000000000000000..ecad810adf48ce703fc9d913bc1b46106a673580
--- /dev/null
+++ b/plot_phases.py
@@ -0,0 +1,15 @@
+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()
diff --git a/read_oskar_beams.py b/read_oskar_beams.py
index a92935d4ca9af6763e903d5ae8a585cda05615ea..d57b3c52b87b46771319d6ac216cf119ae0d7518 100644
--- a/read_oskar_beams.py
+++ b/read_oskar_beams.py
@@ -4,38 +4,70 @@ 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,:,:].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,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
+ 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()
-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,:,:]
- if A is None:
- N = d.shape[-1]
- A = np.zeros((N,N,2,2), dtype=np.complex128)
- A[:,:,i,j] = d
-
-
-
-print(N)
-
-v = np.zeros((N,N))
-
-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_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])
diff --git a/run_oskar.py b/run_oskar.py
index fcc3e25a73a13dab5317f64dabb0a328fb4d6e4f..3d2aabfa8fd6a3de6ce3f07cba29b237c329f9c3 100644
--- a/run_oskar.py
+++ b/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__':