diff --git a/caracal/schema/line_schema.yml b/caracal/schema/line_schema.yml index 137021c18..3d21c4b60 100644 --- a/caracal/schema/line_schema.yml +++ b/caracal/schema/line_schema.yml @@ -411,6 +411,17 @@ mapping: type: bool required: false example: 'False' + pb_type: + desc: Choose between a Gaussian (gauss) primary beam model or the MeerKAT Mauch et al. (2020) model (mauch). + type: str + enum: ['gauss', 'mauch'] + required: false + example: 'gauss' + dish_size: + desc: Dish diameter in meters. Only used in the Gaussian primary beam model + type: float + required: false + example: '13.5' freq_to_vel: desc: Convert the spectral axis' header keys of the line cube from frequency to velocity in the radio definition, v=c(1-obsfreq/restfreq). No change of spectra reference frame is performed. diff --git a/caracal/workers/line_worker.py b/caracal/workers/line_worker.py index b0c7175dd..193206882 100644 --- a/caracal/workers/line_worker.py +++ b/caracal/workers/line_worker.py @@ -127,7 +127,7 @@ def fix_specsys(filename, specframe): del headcube['specsys'] headcube['specsys3'] = specsys3 -def make_pb_cube(filename, apply_corr): +def make_pb_cube(filename, apply_corr, typ, dish_size): if not os.path.exists(filename): caracal.log.warn( 'Skipping primary beam cube for {0:s}. File does not exist.'.format(filename)) @@ -143,11 +143,17 @@ def make_pb_cube(filename, apply_corr): datacube = np.repeat(datacube, headcube['naxis3'], axis=0) * np.abs(headcube['cdelt1']) - sigma_pb = 17.52 / (headcube['crval3'] + headcube['cdelt3'] * ( - np.arange(headcube['naxis3']) - headcube['crpix3'] + 1)) * 1e+9 / 13.5 / 2.355 - # sigma_pb=headcube['crval3']+headcube['cdelt3']*(np.arange(headcube['naxis3'])-headcube['crpix3']+1) - sigma_pb.resize((sigma_pb.shape[0], 1, 1)) - datacube = np.exp(-datacube**2 / 2 / sigma_pb**2) + freq = (headcube['crval3'] + headcube['cdelt3'] * ( + np.arange(headcube['naxis3']) - headcube['crpix3'] + 1)) + if typ == 'gauss': + sigma_pb = 17.52 / (freq / 1e+9) / dish_size / 2.355 + sigma_pb.resize((sigma_pb.shape[0], 1, 1)) + datacube = np.exp(-datacube**2 / 2 / sigma_pb**2) + elif typ == 'mauch': + FWHM_pb = (57.5/60) * (freq / 1.5e9)**-1 + FWHM_pb.resize((FWHM_pb.shape[0], 1, 1)) + datacube = (np.cos(1.189 * np.pi * (datacube / FWHM_pb)) / ( + 1 - 4 * (1.189 * datacube / FWHM_pb)**2))**2 fits.writeto(filename.replace('image.fits','pb.fits'), datacube, header=headcube, overwrite=True) if apply_corr: @@ -1044,7 +1050,10 @@ def worker(pipeline, recipe, config): recipe.add(make_pb_cube, 'make pb_cube-{0:d}'.format(uu), {'filename': image_cube_list[uu], - 'apply_corr': config['pb_cube']['apply_pb'],}, + 'apply_corr': config['pb_cube']['apply_pb'], + 'typ': config['pb_cube']['pb_type'], + 'dish_size': config['pb_cube']['dish_size'], + }, input=pipeline.input, output=pipeline.output, label='Make primary beam cube for {0:s}'.format(image_cube_list[uu])) diff --git a/caracal/workers/selfcal_worker.py b/caracal/workers/selfcal_worker.py index 84d4d3a08..47c89cb2c 100644 --- a/caracal/workers/selfcal_worker.py +++ b/caracal/workers/selfcal_worker.py @@ -561,8 +561,11 @@ def image(trg, num, img_dir, mslist, field): image_opts.update({ "auto-mask": config[key]['cleanmask_thr'][num-1 if len(config[key]['cleanmask_thr']) >= num else -1], "local-rms": config[key]['cleanmask_localrms'][num-1 if len(config[key]['cleanmask_localrms']) >= num else -1], - "local-rms-window": config[key]['cleanmask_localrms_window'][num-1 if len(config[key]['cleanmask_localrms_window']) >= num else -1], }) + if config[key]['cleanmask_localrms'][num-1 if len(config[key]['cleanmask_localrms']) >= num else -1]: + image_opts.update({ + "local-rms-window": config[key]['cleanmask_localrms_window'][num-1 if len(config[key]['cleanmask_localrms_window']) >= num else -1], + }) elif mask_key == 'sofia': fits_mask = 'masking/{0:s}_{1:s}_{2:d}_clean_mask.fits'.format( prefix,field, num)