Browse Source

LBA at 2.5 cm (12 GHz) addition

master
Benito Marcote 1 month ago
parent
commit
08cddf29ea
6 changed files with 70 additions and 25 deletions
  1. +8
    -14
      app.py
  2. BIN
      assets/waves-2_5cm.png
  3. +1
    -1
      data/network_catalog.inp
  4. +14
    -2
      data/stations_catalog.inp
  5. +2
    -1
      vlbiplanobs/freqsetups.py
  6. +45
    -7
      vlbiplanobs/observation.py

+ 8
- 14
app.py View File

@@ -1417,32 +1417,26 @@ def get_fig_uvplane(obs):


def get_fig_dirty_map(obs):
dirty_map_nat, laxis = obs.get_dirtymap(pixsize=1024, robust='natural')
dirty_map_uni, laxis = obs.get_dirtymap(pixsize=1024, robust='uniform')
dirty_map_nat, laxis = obs.get_dirtymap(pixsize=1024, robust='natural', oversampling=4)
dirty_map_uni, laxis = obs.get_dirtymap(pixsize=1024, robust='uniform', oversampling=4)
fig1 = px.imshow(img=dirty_map_nat, x=laxis, y=laxis[::-1], labels={'x': 'RA (mas)', 'y': 'Dec (mas)'}, \
aspect='equal')
fig2 = px.imshow(img=dirty_map_uni, x=laxis, y=laxis[::-1], labels={'x': 'RA (mas)', 'y': 'Dec (mas)'}, \
aspect='equal')
# for i,f in enumerate([fig1, fig2]):
# f.layout.xaxis.autorange = "reversed"
# f.layout.yaxis.autorange = True
# f['layout']['coloraxis']['showscale'] = False

# fig1.layout.xaxis.autorange = "reversed"
# fig2.layout.xaxis.autorange = "reversed"
fig = make_subplots(rows=1, cols=2, subplot_titles=('Natural weighting', 'Uniform weighting'),
shared_xaxes=True, shared_yaxes=True)
fig.add_trace(fig1.data[0], row=1, col=1)
fig.add_trace(fig2.data[0], row=1, col=2)
fig.update_layout(coloraxis={'showscale': False, 'colorscale': 'Inferno'}, showlegend=False, xaxis={'autorange': "reversed"},
yaxis={'autorange': True}, xaxis2={'autorange': "reversed"}, autosize=False)
mapsize = 30*obs.synthesized_beam()['bmaj'].to(u.mas).value
fig.update_layout(coloraxis={'showscale': False, 'colorscale': 'Inferno'}, showlegend=False,
xaxis={'autorange': False, 'range': [mapsize, -mapsize]},
# This xaxis2 represents the xaxis for fig2.
xaxis2={'autorange': False, 'range': [mapsize, -mapsize]},
yaxis={'autorange': False, 'range': [-mapsize, mapsize]}, autosize=False)
fig.update_xaxes(title_text="RA (mas)", constrain="domain")
fig.update_yaxes(title_text="Dec (mas)", row=1, col=1, scaleanchor="x", scaleratio=1)

return fig
return fig1
# layout : xaxis autorange reversed






BIN
assets/waves-2_5cm.png View File

Before After
Width: 851  |  Height: 315  |  Size: 81 KiB

+ 1
- 1
data/network_catalog.inp View File

@@ -43,7 +43,7 @@ observing_bands = 92cm, 49cm, 30cm, 21cm, 18cm, 13cm, 6cm, 5cm, 3.6cm, 1.3cm, 0.
name = Australian Long Baseline Array
default_antennas = ATCA, Pa, Mo, Ho, Cd, Td, Ww
max_datarate = 1024
observing_bands = 21cm, 18cm, 13cm, 6cm, 5cm, 3.6cm, 1.3cm, 0.7cm
observing_bands = 21cm, 18cm, 13cm, 6cm, 5cm, 3.6cm, 2.5cm, 1.3cm, 0.7cm


[KVN]


+ 14
- 2
data/stations_catalog.inp View File

@@ -85,6 +85,7 @@ SEFD_13 = 106.0
SEFD_6 = 70.0
SEFD_5 = 70.0
SEFD_3.6 = 86.0
SEFD_2.5 = 212
SEFD_1.3 = 106
SEFD_0.7 = 180.0

@@ -233,6 +234,7 @@ SEFD_13 = 400.0
SEFD_6 = 450.0
SEFD_5 = 550.0
SEFD_3.6 = 600.0
SEFD_2.5 = 750
SEFD_1.3 = 750
SEFD_0.7 = 2500.0

@@ -474,6 +476,7 @@ SEFD_13 = 410.0
SEFD_6 = 650.0
SEFD_5 = 700.0
SEFD_3.6 = 630.0
SEFD_2.5 = 630.0
SEFD_1.3 = 1800


@@ -496,6 +499,7 @@ SEFD_13 = 650.0
SEFD_6 = 640.0
SEFD_5 = 640.0
SEFD_3.6 = 560.0
SEFD_2.5 = 1200
SEFD_1.3 = 1200
SEFD_0.7 = 1800.0

@@ -593,7 +597,7 @@ SEFD_1.3 = 910
[Katherine]

station = Katherine
code = Ka
code = Ke
network = LBA
possible_networks = LBA
country = Australia
@@ -605,6 +609,8 @@ min_elevation = 10.0
real_time = no
SEFD_21 = 5000.0
SEFD_5 = 4000.0
SEFD_3.6 = 4800.0
SEFD_2.5 = 4800.0


[Kitt Peak]
@@ -913,6 +919,7 @@ SEFD_13 = 530.0
SEFD_6 = 350.0
SEFD_5 = 350.0
SEFD_3.6 = 430.0
SEFD_2.5 = 1300
SEFD_1.3 = 1300
SEFD_0.7 = 675.0

@@ -1072,6 +1079,7 @@ SEFD_13 = 30.0
SEFD_6 = 110.0
SEFD_5 = 110.0
SEFD_3.6 = 43.0
SEFD_2.5 = 370
SEFD_1.3 = 370
SEFD_0.7 = 810.0

@@ -1530,6 +1538,8 @@ min_elevation = 10.0
real_time = no
SEFD_21 = 3400.0
SEFD_5 = 3800.0
SEFD_3.6 = 3800.0
SEFD_2.5 = 3800.0


[Westerbork]
@@ -1599,7 +1609,7 @@ SEFD_3.6 = 750.0
[Yarragadee]

station = Yarragadee
code = Ya
code = Yg
network = LBA
possible_networks = LBA
country = Australia
@@ -1611,6 +1621,8 @@ min_elevation = 10.0
real_time = no
SEFD_21 = 4000.0
SEFD_5 = 3500.0
SEFD_3.6 = 3800.0
SEFD_2.5 = 3800.0


[Yebes]


+ 2
- 1
vlbiplanobs/freqsetups.py View File

@@ -15,7 +15,8 @@ bands = {'92cm': 'P band (92 cm or 0.33 GHz)', '49cm': 'P band (49 cm or 0.6 GHz
'30cm': 'UHF band (30 cm or 1 GHz)', '21cm': 'L band (21 cm or 1.4 GHz)',
'18cm': 'L band (18 cm or 1.7 GHz)', '13cm': 'S band (13 cm or 2.3 GHz)',
'6cm': 'C band (6cm or 5 GHz)', '5cm': 'M band (5 cm or 6 GHz)',
'3.6cm': 'X band (3.6 cm or 8.3 GHz)', '2cm': 'U band (2 cm or 15 GHz)',
'3.6cm': 'X band (3.6 cm or 8.3 GHz)', '2.5cm': 'U band (2.5 cm or 12 GHz)',
'2cm': 'U band (2 cm or 15 GHz)',
'1.3cm': 'K band (1.3 cm or 23 GHz)', #'0.9cm': 'Ka band (0.9 cm or 33 GHz)',
'0.7cm': 'Q band (0.7 cm or 43 GHz)', '0.3cm': 'W band (0.3 cm or 100 GHz)',
'0.1cm': 'Band 3 (0.1 cm or 300 GHz)'}


+ 45
- 7
vlbiplanobs/observation.py View File

@@ -1,4 +1,5 @@
import numpy as np
import scipy.ndimage
import datetime as dt
from astropy import units as u
from astropy import coordinates as coord
@@ -854,7 +855,7 @@ class Observation(object):
return self._synth_beam


def get_dirtymap(self, pixsize:int = 1024, robust: str = "natural"):
def get_dirtymap(self, pixsize:int = 1024, robust: str = "natural", oversampling: int = 4):
"""Returns the dirty beam produced for the given observation.

Input:
@@ -863,6 +864,10 @@ class Observation(object):
- robust : str
The weighting Briggs robust used to compute the dirty map.
It must be either 'natural' or 'uniform'.
- oversampling : int
Oversampling factor when plotting the dirty map.
Recommended values are between 1 and 10.
CURRENTLY THIS PARAMETER IS IGNORED. WILL BE IMPLEMENTED IN A LATER VERSION

Returns:
- dirty_image : (pixsize x pixsize) np.array
@@ -873,21 +878,54 @@ class Observation(object):
"""
assert robust in ('natural', 'uniform')
assert isinstance(pixsize, int)
# assert isinstance(oversampling, int) and 20 > oversampling >= 1
# Oversampling the dirty beam
uvimg = np.zeros((pixsize, pixsize))

# Gridding uv. Inspired from VNSIM code (https://github.com/ZhenZHAO/VNSIM)
# uvdata = self.get_uv_array()
# uvscaling = pixsize/(2.0*1.025*np.max(np.max(np.abs(uvdata), axis=0)))
# uvimg[pixsize//2 + np.trunc(uvdata[:,0]*uvscaling).astype(int), \
# pixsize//2 + np.trunc(uvdata[:,1]*uvscaling).astype(int)] += 1
#
# # Recovering the requested size for the image
# # if oversampling > 1:
# # uvimg = scipy.ndimage.zoom(uvimg, oversampling, order=1)
# # uvimg = uvimg[int(pixsize*0.5):int(pixsize*1.5), int(pixsize*0.5):int(pixsize*1.5)]
# if robust == 'uniform':
# if np.max(uvimg) > 1:
# uvimg[np.where(uvimg > 0)] = 1
# else:
# print('\n\nWARNING: Uniform dirty map same as natural one.\n\n')
#
# dirty_beam = np.real(np.fft.ifftshift(np.fft.ifft2(np.fft.fftshift(uvimg/np.max(uvimg)))))
# if oversampling > 1:
# dirty_beam = scipy.ndimage.zoom(dirty_beam, oversampling, order=1)
# imgsize = (uvscaling*u.rad/(2*oversampling)).to(u.mas) # angular equivalent size of the resulting image
# return dirty_beam[int(pixsize*(oversampling-1)/2):int(pixsize*(oversampling+1)/2), \
# int(pixsize*(oversampling-1)/2):int(pixsize*(oversampling+1)/2)].T/np.max(dirty_beam), \
# np.linspace(-imgsize, imgsize, pixsize)
# # return dirty_beam.T/np.max(dirty_beam), \
#

oversampling = 1
# Gridding uv. Inspired from VNSIM code (https://github.com/ZhenZHAO/VNSIM)
uvdata = self.get_uv_array()
uvscaling = pixsize/(2*np.max(np.max(uvdata, axis=0)))
uvscaling = (oversampling*pixsize)/(2*1.05*np.max(np.max(np.abs(uvdata), axis=0)))
if robust == 'natural':
for uv in uvdata:
uvimg[int(pixsize//2 + round(uv[0]*uvscaling)), int(pixsize//2 + round(uv[1]*uvscaling))] += 1
uvimg[oversampling*pixsize//2 + np.trunc(uvdata[:,0]*uvscaling).astype(int), \
oversampling*pixsize//2 + np.trunc(uvdata[:,1]*uvscaling).astype(int)] += 1
else:
for uv in uvdata:
uvimg[int(pixsize//2 + round(uv[0]*uvscaling)), int(pixsize//2 + round(uv[1]*uvscaling))] = 1
uvimg[oversampling*pixsize//2 + np.trunc(uvdata[:,0]*uvscaling).astype(int), \
oversampling*pixsize//2 + np.trunc(uvdata[:,1]*uvscaling).astype(int)] = 1

# Recovering the requested size for the image
# uvimg = uvimg[int(pixsize*0.5):int(pixsize*1.5), int(pixsize*0.5):int(pixsize*1.5)]
dirty_beam = np.real(np.fft.ifftshift(np.fft.ifft2(np.fft.fftshift(uvimg/np.max(uvimg)))))
imgsize = (uvscaling*u.rad/2).to(u.mas) # angular equivalent size of the resulting image
return dirty_beam.T/np.max(dirty_beam), np.linspace(-imgsize, imgsize, pixsize)
return dirty_beam[int(pixsize*(oversampling-1)/2):int(pixsize*(oversampling+1)/2), \
int(pixsize*(oversampling-1)/2):int(pixsize*(oversampling+1)/2)].T/np.max(dirty_beam), \
np.linspace(-imgsize, imgsize, pixsize)





Loading…
Cancel
Save