Skip to content
Snippets Groups Projects
Commit c1d36eda authored by Bas van der Tol's avatar Bas van der Tol
Browse files

Add scripts and data to run OSKAR for comparison

parent 578c7f5c
No related branches found
No related tags found
No related merge requests found
Showing
with 18400 additions and 0 deletions
#!/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
output_dir = 'test.tm'
element_id = 0
freqs = [50]
l_max = 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
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
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=',')
from matplotlib import pyplot as plt
from astropy.io import fits
import numpy as np
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])
#!/usr/bin/env python3
"""
Generates all-sky zenith-centred beam patterns for SKALA-4 and EDA-2 antennas.
"""
import copy
import glob
import os.path
import re
import subprocess
from astropy.io import fits
from astropy.time import Time, TimeDelta
import matplotlib
matplotlib.use('Agg')
# pylint: disable=wrong-import-position
import matplotlib.pyplot as plt
import numpy
import oskar
def get_start_time(ra0_deg, length_sec):
"""Returns optimal start time for field RA and observation length."""
t = Time('2000-01-01 00:00:00', scale='utc', location=('116.764d', '0d'))
dt_hours = (24.0 - t.sidereal_time('apparent').hour) / 1.0027379
dt_hours += (ra0_deg / 15.0)
start = t + TimeDelta(dt_hours * 3600.0 - length_sec / 2.0, format='sec')
return start.value
def plot_panel(ax, image, title, cmap):
"""Plots a single panel."""
im = ax.imshow(numpy.squeeze(image), cmap=cmap)
plt.colorbar(im, format='%.2e')
plt.tick_params(labelcolor='none', top='off', bottom='off', left='off', right='off')
ax.set_xticks([])
ax.set_yticks([])
ax.invert_yaxis()
ax.set_title(title)
ax.axis('equal')
def make_plots(title, glob_pattern, out_basename):
"""Generates a plot with four panels."""
# Load FITS images matching the glob pattern.
files = glob.glob(glob_pattern)
files.sort() # Must be sorted!
images = []
for file in files:
images.append(fits.getdata(file))
# Set titles and colour maps to use.
cmap = ''
titles = []
if '_AUTO_POWER' in glob_pattern:
# cmap = 'plasma'
cmap = 'Blues_r'
titles = ['Stokes I', 'Stokes Q', 'Stokes U', 'Stokes V']
elif '_AMP' in glob_pattern:
# cmap = 'jet'
cmap = 'CMRmap'
titles = ['XX', 'XY', 'YX', 'YY']
# Sub-plots, one for each Stokes parameter.
fig = plt.figure(figsize=(8, 6))
ax = fig.add_subplot(221, frameon=False)
plot_panel(ax, images[0], titles[0], cmap)
ax = fig.add_subplot(222, frameon=False)
plot_panel(ax, images[1], titles[1], cmap)
ax = fig.add_subplot(223, frameon=False)
plot_panel(ax, images[2], titles[2], cmap)
ax = fig.add_subplot(224, frameon=False)
plot_panel(ax, images[3], titles[3], cmap)
# Add main title.
title = title.replace('_', ' ') # Replace underscores with spaces.
fig.suptitle(title)
fig.tight_layout()
fig.subplots_adjust(top=0.88)
# Save and close.
plt.savefig('%s.png' % (out_basename))
plt.close('all')
def main():
"""Main function."""
# Name of the application to run, and a settings file for it.
app = 'oskar_sim_beam_pattern'
settings_path = '_temp_settings.ini'
# Define some basic observation parameters.
ra0_deg = 0.0
dec0_deg = -27.0
length_sec = 1.0
# Define base settings dictionary.
common_settings = {
'observation': {
'phase_centre_ra_deg': ra0_deg,
'phase_centre_dec_deg': dec0_deg,
'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
},
'beam_pattern': {
'coordinate_frame': 'Horizon',
'beam_image/size': 256,
'station_outputs/fits_image/amp': True,
'station_outputs/fits_image/phase': True,
'station_outputs/fits_image/auto_power_real': False
}
}
# Define frequencies of interest (in MHz).
freqs = [50]
#freqs = [70]
# Define telescope models to use, and associated overrides for them.
telescopes = {
'basefunctions': {
'telescope/input_directory': 'test.tm',
'telescope/aperture_array/element_pattern/enable_numerical': True,
'telescope/aperture_array/element_pattern/swap_xy': True,
'telescope/aperture_array/array_pattern/enable': False
},
}
# Loop over telescope models.
for tel, tel_params in telescopes.items():
# Copy the base settings dictionary.
current_settings = copy.deepcopy(common_settings)
# Update current settings with telescope model parameters.
current_settings.update(tel_params)
# Loop over frequencies.
for freq in freqs:
# Create the settings file.
open(settings_path, 'w').close()
settings = oskar.SettingsTree(app, settings_path)
settings.from_dict(current_settings)
# Update output root path and frequency.
tel_root = re.sub(r'[^\w]', '', tel) # Strip symbols from tel.
root_path = tel_root + ('_%03d_MHz' % freq)
settings['beam_pattern/root_path'] = root_path
settings['observation/start_frequency_hz'] = 1e6 * freq
# Run the app with the settings file.
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')
if __name__ == '__main__':
main()
0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00
0.000000000000000000e+00,0.000000000000000000e+00,1.000000000000000000e+00
0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00
0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00
0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00
0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00
0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00
0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00
0,0
-605.1565,146.9829
-64.79,-577.8771
-878.9322,-123.8241
-237.3169,-340.6209
407.1441,198.459
139.1693,-888.4022
345.5998,890.5183
-203.9332,11.27668
-404.5717,373.44
136.5508,745.0029
-375.0552,-139.6862
254.7167,-507.6426
100.512,-355.1444
-360.9658,161.2952
698.7782,-471.8213
210.4574,-100.3645
-969.0271,190.9689
893.7347,117.048
-844.7219,-321.3391
558.4537,469.6326
-40.34307,530.6472
-148.6963,281.2029
-543.7159,731.1521
256.6313,445.902
125.0,125.0
100.0,100.0
50.0,50.0
25.0,25.0
10.0,10.0
0 -50
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment