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 3 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