Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Mauchian #1215

Merged
merged 6 commits into from
Jul 31, 2020
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
11 changes: 11 additions & 0 deletions caracal/schema/line_schema.yml
Original file line number Diff line number Diff line change
Expand Up @@ -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.
Expand Down
23 changes: 16 additions & 7 deletions caracal/workers/line_worker.py
Original file line number Diff line number Diff line change
Expand Up @@ -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))
Expand All @@ -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:
Expand Down Expand Up @@ -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]))
Expand Down
5 changes: 4 additions & 1 deletion caracal/workers/selfcal_worker.py
Original file line number Diff line number Diff line change
Expand Up @@ -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]:
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This isjust cosmetics of course, but in the interest of readability (who will think of the future generations!), I would prefer to phrase that as

local_rms = config[key]['cleanmask_localrms'][min(num-1, len(config[key]['cleanmask_localrms'])-1)]
if local_rms:
    image_opts.update({
        "local-rms": local_rms, 
        "local-rms-window": config[key]['cleanmask_localrms_window'][min(num-1, len(config[key]['cleanmask_localrms'])-1)]
    })

...since that if-else construct is a big eyeful to take in, and recurs a lot.

Even better, define a helper function in there:

def n_1_or_last(seq, num):
    return seq[min(num-1, len(seq)-1)]

...since I presume this treatment of lists happens all over the place?

Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

yeah that would be neater -- we can sneak that into the next PR

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)
Expand Down