diff --git a/demo/multi_resolve/beam_to_solid_angle.py b/demo/multi_resolve/beam_to_solid_angle.py new file mode 100644 index 0000000000000000000000000000000000000000..dfea55cf6de78497e7f31ff378aa16e69b9dc47a --- /dev/null +++ b/demo/multi_resolve/beam_to_solid_angle.py @@ -0,0 +1,9 @@ +import numpy as np + + +def pix_to_beam(beam_size): + beam_solid_angle = np.pi * (beam_size / 2) ** 2 + cell_size = beam_size / 5 + pixel_solid_angle = cell_size ** 2 + pix_to_beam = beam_solid_angle / pixel_solid_angle + return pix_to_beam diff --git a/demo/multi_resolve/mosaic_correct.py b/demo/multi_resolve/mosaic_correct.py index 85b51e90498beac2941be7bde8058865a542ee56..c199199a6966530e01c695200f43a344017f0551 100644 --- a/demo/multi_resolve/mosaic_correct.py +++ b/demo/multi_resolve/mosaic_correct.py @@ -60,10 +60,11 @@ sdom = sky.target[3] x_direction = coords(sdom.shape[0], sdom.distances[0]) y_direction = coords(sdom.shape[1], sdom.distances[1]) sky_coordinates = np.array(np.meshgrid( - x_direction, y_direction, indexing='xy')) + # -x_direction, y_direction, indexing='xy' + x_direction, y_direction, indexing='xy' +)) - -output_directory = f"output/combined_corrected_{npix}" +output_directory = f"output/combined_corrected_{npix}_old" sky_center = SkyCoord(center_ra, center_dec, frame=center_frame) @@ -84,12 +85,7 @@ for fldid, oo in enumerate(all_obs): beam = ift.makeField(sdom, beam) # beam_operator = ift.DiagonalOperator(beam) - beam_direction = f'fld{fldid}' - beam_directions[beam_direction] = dict( - dx=dx, - dy=dy, - beam=beam - ) + beam_directions[f'fld{fldid}'] = dict(dx=dx, dy=dy, beam=beam) # Used for the dtypes @@ -131,10 +127,7 @@ def build_response(field_key, dx, dy, obs, sky_dtype): domain_dtype=obs.vis.dtype ) - # FIXME: This response operator makes unnecessary multiple beam pattern - # multiplications. As the SKY_BEAMER calculates the beam pattern for all fields. - # Hence, we throw away all other fields. - return UPCAST_TO_STOKES_I @ RADIO_RESPONSE @ FIELD_EXTRACTOR @ SKY_BEAMER @ REDUCER + return UPCAST_TO_STOKES_I @ RADIO_RESPONSE @ FIELD_EXTRACTOR def build_energy(response, obs): @@ -149,7 +142,8 @@ responses = [] energies = [] for kk, vv, obs in zip(beam_directions.keys(), beam_directions.values(), all_obs): R = build_response(kk, vv['dx'], vv['dy'], obs, tmp_sky.dtype) - responses.append(R) + + responses.append(R @ SKY_BEAMER @ REDUCER) energies.append(build_energy(R, obs)) @@ -159,12 +153,16 @@ for ii in range(1): tot_response_adjoint = np.zeros_like(f.val) for rr, oo in zip(responses, all_obs): - tot_response_adjoint += rr.adjoint(oo.vis).val + dirty = rr.adjoint(oo.vis).val + tot_response_adjoint += dirty + plt.imshow(dirty[0, 0, 0], origin='lower') + plt.show() + plt.imshow(tot_response_adjoint[0, 0, 0], origin='lower') plt.show() lh = reduce(lambda x, y: x+y, energies) -lh = lh @ sky +lh = lh @ SKY_BEAMER @ REDUCER @ sky def callback(samples, i): diff --git a/demo/multi_resolve/ms2rve.py b/demo/multi_resolve/ms2rve.py index 47cd3e00fe67bf08c6aa27601bfd603c1ba4313e..189a67c3a6aca37cec76d2ccb6c78345aa2625c5 100644 --- a/demo/multi_resolve/ms2rve.py +++ b/demo/multi_resolve/ms2rve.py @@ -19,7 +19,11 @@ outdir = args.output_dir makedirs(outdir, exist_ok=True) -name = input.split('/')[-1] + '_fld{field:02d}' -obs = rve.ms2observations(input, "DATA", True, 0) +spw = 1 + +name = input.split('/')[-1] + '_fld{field:02d}_spw{spw}' +obs = rve.ms2observations(input, "DATA", True, spw) for ii, o in enumerate(obs): - o.save(join(outdir, name.format(field=ii)), False) + if o is None: + continue + o.save(join(outdir, name.format(field=ii, spw=spw)), False) diff --git a/demo/multi_resolve/multi_obs.cfg b/demo/multi_resolve/multi_obs.cfg index 9903815c3f19d52c674c0b5feaef19be7587eb07..a45185b4a0a4f10f494548fa4294adba5ca75b00 100644 --- a/demo/multi_resolve/multi_obs.cfg +++ b/demo/multi_resolve/multi_obs.cfg @@ -1,8 +1,8 @@ [sky] freq mode = single polarization=I -space npix x = 256 -space npix y = 256 +space npix x = 128 +space npix y = 128 space fov x = 25as space fov y = 25as # space fov x = 30.72as