diff --git a/demo/comparison-oskar/run_oskar_simulation.py b/demo/comparison-oskar/run_oskar_simulation.py new file mode 100755 index 0000000000000000000000000000000000000000..d05aff0b3f826e512e19dd7c60f7b3e648f91ccd --- /dev/null +++ b/demo/comparison-oskar/run_oskar_simulation.py @@ -0,0 +1,122 @@ +#!/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 main(npixels): + """Main function.""" + # Name of the application to run, and a settings file for it. + app1 = 'oskar_sim_interferometer' + app2 = 'oskar_sim_beam_pattern' + settings_path = '_temp_settings.ini' + + # Define some basic observation parameters. + ra0_deg = 0.0 + dec0_deg = -27.0 + length_sec = 3600.0*6 + + # Define base settings dictionary. + common_settings = { + '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, + 'num_time_steps': 60 + + }, + 'telescope': { + 'normalise_beams_at_phase_centre': False, + 'aperture_array/element_pattern/normalise': False + } + } + + # Define frequencies of interest (in MHz). + freqs = [50] + #freqs = [70] + + # Define telescope models to use, and associated overrides for them. + telescopes = { + #'stationresponse': { + 'basefunctions': { + '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': True, + }, + } + + # 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(app1, 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['observation/start_frequency_hz'] = 1e6 * freq + settings['interferometer/ms_filename'] = 'blah.ms' + + # Run the app with the settings file. + subprocess.call([app1, settings_path]) + + # Create the settings file. + open(settings_path, 'w').close() + settings = oskar.SettingsTree(app2, 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 + settings['beam_pattern/coordinate_frame'] = 'Horizon' + settings['beam_pattern/beam_image/size'] = npixels + settings['beam_pattern/station_outputs/fits_image/amp'] = True + settings['beam_pattern/station_outputs/fits_image/phase'] = True + settings['beam_pattern/station_outputs/fits_image/auto_power_real'] = False + + # Run the app with the settings file. + subprocess.call([app2, settings_path]) + + +if __name__ == '__main__': + main(32) diff --git a/demo/comparison-oskar/stationresponse.cpp b/demo/comparison-oskar/stationresponse.cpp index cdb944e0008d247abd4f94b9d1fa38d41ce10c1c..9898f3a5a8244fb1f5beb2857a873642e2085c40 100644 --- a/demo/comparison-oskar/stationresponse.cpp +++ b/demo/comparison-oskar/stationresponse.cpp @@ -20,7 +20,8 @@ int main(int argc, char** argv){ Options options; options.element_response_model = response_model; - casacore::MeasurementSet ms("/home/vdtol/skalowmini/skalowmini-coef1.MS"); +// casacore::MeasurementSet ms("/home/vdtol/skalowmini/skalowmini-coef1.MS"); + casacore::MeasurementSet ms("/home/vdtol//src/EveryBeam/build-release/demo/comparison-oskar/blah.ms"); casacore::ScalarMeasColumn<casacore::MDirection> referenceDirColumn( ms.field(), @@ -88,11 +89,16 @@ int main(int argc, char** argv){ }; real_t freq0 = 50e6; - auto result = station->Response(time, freq, direction, freq0, station0, tile0); - result_arr[i][j][0][0] = result[0][0]; - result_arr[i][j][0][1] = result[0][1]; - result_arr[i][j][1][0] = result[1][0]; - result_arr[i][j][1][1] = result[1][1]; +// auto result = station->Response(time, freq, direction, freq0, station0, tile0); +// result_arr[i][j][0][0] = result[0][0]; +// result_arr[i][j][0][1] = result[0][1]; +// result_arr[i][j][1][0] = result[1][0]; +// result_arr[i][j][1][1] = result[1][1]; + auto result = station->ArrayFactor(time, freq, direction, freq0, station0, tile0); + result_arr[i][j][0][0] = result[0]; + result_arr[i][j][0][1] = 0.0; + result_arr[i][j][1][0] = 0.0; + result_arr[i][j][1][1] = result[1]; } }