diff --git a/meqsilhouette/driver/run_meqsilhouette.py b/meqsilhouette/driver/run_meqsilhouette.py index 6d656e4..620e364 100644 --- a/meqsilhouette/driver/run_meqsilhouette.py +++ b/meqsilhouette/driver/run_meqsilhouette.py @@ -203,7 +203,7 @@ def run_meqsilhouette(config=None): ### TROPOSPHERE COMPONENTS ### combined_phase_errors = 0 #init for trop combo choice - additive_noises = None + if parameters['trop_enabled']: info('Tropospheric module is enabled, applying corruptions...') if parameters['trop_wetonly']: @@ -215,16 +215,11 @@ def run_meqsilhouette(config=None): info('TROPOSPHERE ATTENUATE: attenuating signal using PWV-derived opacity...') sim_coord.trop_opacity_attenuate() - if parameters['trop_noise']: - info('TROPOSPHERE NOISE: adding sky noise from non-zero PWV...') - additive_noises = sim_coord.trop_add_sky_noise() - if parameters['trop_mean_delay']: info('TROPOSPHERE DELAY: computing mean delay (time-variability from elevation changes)...') sim_coord.trop_calc_mean_delays() combined_phase_errors += sim_coord.phasedelay_alltimes - if parameters['trop_turbulence']: info('TROPOSPHERE TURBULENCE: computing Kolmogorov turbulence phase errors...') sim_coord.trop_generate_turbulence_phase_errors() @@ -245,9 +240,6 @@ def run_meqsilhouette(config=None): sim_coord.trop_plots() info('Generated troposphere plots') - ### POPULATE MS WITH SIGMA AND WEIGHT ESTIMATORS ### - sim_coord.add_weights(additive_noises) - ### PARALLACTIC ANGLE AND POLARIZATION LEAKAGE ### if parameters['uvjones_d_on']: info('Introducing parallactic angle rotation and polarization leakage effects') @@ -271,11 +263,9 @@ def run_meqsilhouette(config=None): info('Generating bandpass plots...') sim_coord.make_bandpass_plots() - ### THERMAL NOISE ### - if parameters['add_thermal_noise']: - info('Adding thermal noise...') - sim_coord.add_receiver_noise() - info('Thermal noise added.') + ### Add all noise components and fill in the weight columns + # TODO must do the job of trop_add_sky_noise(), add_weights(additive_noises), add_receiver_noise() here + sim_coord.add_noise(parameters['trop_noise'], parameters['add_thermal_noise']) ### IMAGING, PLOTTING, DATA EXPORT ### if parameters['make_image']: diff --git a/meqsilhouette/framework/SimCoordinator.py b/meqsilhouette/framework/SimCoordinator.py index 443dc90..e545dc5 100644 --- a/meqsilhouette/framework/SimCoordinator.py +++ b/meqsilhouette/framework/SimCoordinator.py @@ -114,28 +114,25 @@ def __init__(self, msname, output_column, input_fitsimage, input_fitspol, input_ self.thermal_noise_enabled = thermal_noise_enabled self.receiver_rms = np.zeros(self.data.shape, dtype='float') - self.compute_weights() - tab.close() # close main MS table - + ### troposphere information self.trop_enabled = trop_enabled - if (self.trop_enabled): - self.trop_wetonly = trop_wetonly - self.average_pwv = pwv - self.average_gpress = gpress - self.average_gtemp = gtemp - self.coherence_time = coherence_time - self.fixdelay_max_picosec = fixdelay_max_picosec - self.elevation_tropshape = np.expand_dims(np.swapaxes(self.elevation, 0, 1), 1) # reshaped for troposphere operations - self.opacity, self.sky_temp = self.trop_return_opacity_sky_temp() - self.transmission = np.exp(-1*self.opacity) - - # Set some optional arrays to None. These will be filled later depending upon the user request. - self.transmission_matrix = None - self.turb_phase_errors = None - self.delay_alltimes = None - self.sky_noise = None + self.trop_wetonly = trop_wetonly + self.average_pwv = pwv + self.average_gpress = gpress + self.average_gtemp = gtemp + self.coherence_time = coherence_time + self.fixdelay_max_picosec = fixdelay_max_picosec + self.elevation_tropshape = np.expand_dims(np.swapaxes(self.elevation, 0, 1), 1) # reshaped for troposphere operations + self.opacity, self.sky_temp = self.trop_return_opacity_sky_temp() + self.transmission = np.exp(-1*self.opacity) + + # Set some optional arrays to None. These will be filled later depending upon the user request. + self.transmission_matrix = None + self.turb_phase_errors = None + self.delay_alltimes = None + self.sky_noise = None ### bandpass information self.bandpass_table = bandpass_table @@ -1091,6 +1088,76 @@ def add_gjones_manual(self): self.save_data() + ################################## + # Add noise components + def add_noise(self, tropnoise, thermalnoise): + """ Add sky and receiver noise components to the visibilities and populate weight columns""" + + # compute receiver rms noise + for a0 in range(self.Nant): + for a1 in range(self.Nant): + if a1 > a0: + self.receiver_rms[self.baseline_dict[(a0,a1)]] = (1/self.corr_eff) * np.sqrt(self.SEFD[a0] * \ + self.SEFD[a1] / float(2 * self.tint * self.chan_width)) + + # compute sky sigma estimator (i.e. sky rms noise) and realise sky noise + if tropnoise: + sefd_matrix = 2 * Boltzmann / self.dish_area * (1e26*self.sky_temp * (1. - np.exp(-1.0 * self.opacity / np.sin(self.elevation_tropshape)))) + self.sky_noise = np.zeros(self.data.shape, dtype='complex') + sky_sigma_estimator = np.zeros(self.data.shape) + + for a0 in range(self.Nant): + for a1 in range(self.Nant): + if a1 > a0: + rms = (1/self.corr_eff) * np.sqrt(sefd_matrix[:, :, a0] * sefd_matrix[:, :, a1] / (float(2 * self.tint * self.chan_width))) + self.temp_rms = rms + rms = np.expand_dims(rms, 2) + rms = rms * np.ones((1, 1, 4)) + self.sky_noise[self.baseline_dict[(a0, a1)]] = self.rng_atm.normal(0.0, rms) + 1j * self.rng_atm.normal(0.0, rms) + sky_sigma_estimator[self.baseline_dict[(a0, a1)]] = rms + np.save(II('$OUTDIR')+'/atm_output/sky_noise_timestamp_%d'%(self.timestamp), self.sky_noise) + + # add sky noise rms to receiver rms if tropnoise is set + try: + for tind in range(self.nchunks): + self.receiver_rms[tind*self.chunksize:(tind+1)*self.chunksize] = np.sqrt(np.power(self.receiver_rms[tind*self.chunksize:(tind+1)*self.chunksize], 2) + \ + np.power(sky_sigma_estimator[tind*self.chunksize:(tind+1)*self.chunksize], 2)) + except MemoryError: + abort("Arrays too large to be held in memory. Aborting execution.") + + # use the receiver rms to realise thermal noise + if thermalnoise: + info('Realise thermal noise from receiver rms...') + self.thermal_noise = np.zeros(self.data.shape, dtype='complex') + size = (self.time_unique.shape[0], self.chan_freq.shape[0], 4) + for a0 in range(self.Nant): + for a1 in range(self.Nant): + if a1 > a0: + rms = self.receiver_rms[self.baseline_dict[(a0,a1)]] + self.thermal_noise[self.baseline_dict[(a0, a1)]] = self.rng_predict.normal(0.0, rms, size=size) + 1j * self.rng_predict.normal(0.0, rms, size=size) + + np.save(II('$OUTDIR')+'/receiver_noise_timestamp_%d'%(self.timestamp), self.thermal_noise) + + # add both the noise components to the data + try: + info('Applying thermal noise to data...') + for tind in range(self.nchunks): + self.data[tind*self.chunksize:(tind+1)*self.chunksize] += self.thermal_noise[tind*self.chunksize:(tind+1)*self.chunksize] + \ + self.sky_noise[tind*self.chunksize:(tind+1)*self.chunksize] + self.save_data() + except MemoryError: + abort("Arrays too large to be held in memory. Aborting execution.") + + # populate MS weight and sigma columns using the receiver rms (with or without sky noise) + tab = pt.table(self.msname, readonly=False,ack=False) + tab.putcol("SIGMA", self.receiver_rms[:,0,:]) + if 'SIGMA_SPECTRUM' in tab.colnames(): + tab.putcol("SIGMA_SPECTRUM", self.receiver_rms) + tab.putcol("WEIGHT", 1/self.receiver_rms[:,0,:]**2) + if 'WEIGHT_SPECTRUM' in tab.colnames(): + tab.putcol("WEIGHT_SPECTRUM", 1/self.receiver_rms**2) + tab.close() + ############################# ##### General MS plots ##### #############################