Skip to content

Commit fb43e99

Browse files
timstaleytimstaley
committed
Add simple simulation and imaging notebook.
1 parent 00dbd1a commit fb43e99

3 files changed

Lines changed: 183 additions & 5 deletions

File tree

‎docs/source/notebooks.rst‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -6,3 +6,4 @@ Jupyter Notebook Examples
66

77
notebooks/convolutional_kernel_funcs
88
notebooks/anti_aliasing_properties
9+
notebooks/simulate_and_image

‎docs/source/notebooks/anti_aliasing_properties.ipynb‎

Lines changed: 3 additions & 5 deletions
Original file line numberDiff line numberDiff line change
@@ -99,9 +99,7 @@
9999
"ax = fig.gca()\n",
100100
"# plot_section = np.log10(amp[:oversampling])\n",
101101
"ax.plot(np.arange(len(amp))-support, amp)\n",
102-
"ax.grid()\n",
103-
"\n",
104-
"\n"
102+
"ax.grid()"
105103
]
106104
},
107105
{
@@ -112,7 +110,7 @@
112110
},
113111
"outputs": [],
114112
"source": [
115-
"C_grid = np.fft.fftshift(np.fft.ifft2(np.fft.ifftshift(gs_kernel.array)))\n"
113+
"C_grid = np.fft.fftshift(np.fft.ifft2(np.fft.ifftshift(gs_kernel.array)))"
116114
]
117115
},
118116
{
@@ -140,5 +138,5 @@
140138
],
141139
"metadata": {},
142140
"nbformat": 4,
143-
"nbformat_minor": 1
141+
"nbformat_minor": 0
144142
}

‎docs/source/notebooks/simulate_and_image.ipynb‎

Lines changed: 179 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,179 @@
1+
{
2+
"cells": [
3+
{
4+
"cell_type": "markdown",
5+
"metadata": {},
6+
"source": [
7+
"# Example visibility-simulation and imaging"
8+
]
9+
},
10+
{
11+
"cell_type": "code",
12+
"execution_count": null,
13+
"metadata": {
14+
"collapsed": true
15+
},
16+
"outputs": [],
17+
"source": [
18+
"import numpy as np\n",
19+
"import logging\n",
20+
"import os\n",
21+
"\n",
22+
"import astropy.units as u\n",
23+
"import astropy.constants as const\n",
24+
"import fastimgproto.imager as imager\n",
25+
"import fastimgproto.visibility as visibility\n",
26+
"\n",
27+
"from astropy.coordinates import Angle, SkyCoord\n",
28+
"from astropy.time import Time\n",
29+
"from fastimgproto.gridder.conv_funcs import GaussianSinc\n",
30+
"from fastimgproto.skymodel.helpers import SkyRegion, SkySource\n",
31+
"from fastimgproto.sourcefind.image import SourceFindImage\n",
32+
"from fastimgproto.telescope.readymade import Meerkat\n"
33+
]
34+
},
35+
{
36+
"cell_type": "code",
37+
"execution_count": null,
38+
"metadata": {
39+
"collapsed": true
40+
},
41+
"outputs": [],
42+
"source": [
43+
"import matplotlib.pyplot as plt\n",
44+
"import matplotlib.image as mpimg\n",
45+
"%matplotlib inline\n",
46+
"# Plot image pixels in cartesian ordering (i.e. y-positive == upwards):\n",
47+
"plt.rcParams['image.origin'] = 'lower'\n",
48+
"# Make plots bigger\n",
49+
"plt.rcParams['figure.figsize'] = 10, 10"
50+
]
51+
},
52+
{
53+
"cell_type": "code",
54+
"execution_count": null,
55+
"metadata": {
56+
"collapsed": false
57+
},
58+
"outputs": [],
59+
"source": [
60+
"telescope = Meerkat()\n",
61+
"print(\"Telescope with {} antennae\".format(len(telescope.ant_local_xyz)))\n",
62+
"print(\"Centre: {!r}, {!r}\".format(telescope.lon, telescope.lat))"
63+
]
64+
},
65+
{
66+
"cell_type": "code",
67+
"execution_count": null,
68+
"metadata": {
69+
"collapsed": false
70+
},
71+
"outputs": [],
72+
"source": [
73+
"pointing_centre = SkyCoord(0 * u.deg, -30 * u.deg)\n",
74+
"obs_central_frequency = 3. * u.GHz\n",
75+
"wavelength = const.c / obs_central_frequency\n",
76+
"transit_time = telescope.next_transit(pointing_centre.ra,\n",
77+
" start_time=Time('2017-01-01'))"
78+
]
79+
},
80+
{
81+
"cell_type": "code",
82+
"execution_count": null,
83+
"metadata": {
84+
"collapsed": false
85+
},
86+
"outputs": [],
87+
"source": [
88+
"nstep=1\n",
89+
"obs_times = transit_time + np.linspace(-1, 1, nstep) * u.hr\n",
90+
"print(\"Generating UVW-baselines for {} timesteps\".format(nstep))\n",
91+
"uvw_m = telescope.uvw_tracking_skycoord(pointing_centre, obs_times)\n",
92+
"# From here on we use UVW as multiples of wavelength, lambda:\n",
93+
"uvw_lambda = (uvw_m / wavelength).to(u.dimensionless_unscaled).value\n",
94+
"\n",
95+
"\n",
96+
"\n",
97+
"# Additional source to North-East of pointing centre\n",
98+
"extra_src_position = SkyCoord(ra=pointing_centre.ra + 0.01 * u.deg,\n",
99+
" dec=pointing_centre.dec + 0.01 * u.deg, )\n",
100+
"\n",
101+
"steady_sources = [\n",
102+
" SkySource(pointing_centre, flux=1 * u.Jy),\n",
103+
" SkySource(extra_src_position, flux=0.4 * u.Jy),\n",
104+
"]\n",
105+
"\n",
106+
"# Simulate incoming data; includes transient sources, noise:\n",
107+
"print(\"Simulating visibilities\")\n",
108+
"data_vis = visibility.visibilities_for_source_list(\n",
109+
" pointing_centre,\n",
110+
" source_list = steady_sources, \n",
111+
" uvw = uvw_lambda)\n",
112+
"\n",
113+
"vis_noise_level = 0.1 * u.Jy\n",
114+
"data_vis = visibility.add_gaussian_noise(vis_noise_level, data_vis)\n"
115+
]
116+
},
117+
{
118+
"cell_type": "code",
119+
"execution_count": null,
120+
"metadata": {
121+
"collapsed": false
122+
},
123+
"outputs": [],
124+
"source": [
125+
"print(\"Simulated {} visibilities\".format(len(data_vis)))"
126+
]
127+
},
128+
{
129+
"cell_type": "code",
130+
"execution_count": null,
131+
"metadata": {
132+
"collapsed": true
133+
},
134+
"outputs": [],
135+
"source": [
136+
"image_size=1024 * u.pixel\n",
137+
"cell_size=1 * u.arcsecond"
138+
]
139+
},
140+
{
141+
"cell_type": "code",
142+
"execution_count": null,
143+
"metadata": {
144+
"collapsed": false
145+
},
146+
"outputs": [],
147+
"source": [
148+
"kernel_support = 3\n",
149+
"kernel_func = GaussianSinc(trunc=kernel_support)\n",
150+
"image, beam = imager.image_visibilities(data_vis, uvw_lambda,\n",
151+
" image_size=image_size,\n",
152+
" cell_size=cell_size,\n",
153+
" kernel_func=kernel_func,\n",
154+
" kernel_support=kernel_support,\n",
155+
" kernel_exact=True,\n",
156+
" kernel_oversampling=None\n",
157+
" )"
158+
]
159+
},
160+
{
161+
"cell_type": "code",
162+
"execution_count": null,
163+
"metadata": {
164+
"collapsed": false
165+
},
166+
"outputs": [],
167+
"source": [
168+
"fig, ax1 = plt.subplots(1,1)\n",
169+
"ax1.imshow(np.real(image))\n",
170+
"axlims = 400, 600\n",
171+
"ax1.set_xlim(axlims)\n",
172+
"ax1.set_ylim(axlims)"
173+
]
174+
}
175+
],
176+
"metadata": {},
177+
"nbformat": 4,
178+
"nbformat_minor": 2
179+
}

0 commit comments

Comments
 (0)