casa6
Public
Source
1
+
import numpy as np2
+
import shutil3
+
import os4
+
from io import StringIO5
+
import math, re, sys, logging6
+
from functools import reduce7
+
import time as tm8
+
import datetime as dt9
+
import casatools10
+
msmd = casatools.msmetadata()11
+
tb = casatools.table()12
+
13
+
root = "/export/home/sazed/nschweig/import_casa_testing_py3/CASA6/ANTAB/casa-6.6.4-17-py3.10.el7/bin/"14
+
15
+
#idifile = root + "VLBA.idifits"16
+
idi2ms = root + "idi1.ms"17
+
antabfile = root + "n14c3.antab"18
+
idicopy = "idicopy.ms"19
+
20
+
desc = {'ANTENNA_ID': {'comment': 'ID of antenna in this array',21
+
'dataManagerGroup': 'StandardStMan',22
+
'dataManagerType': 'StandardStMan',23
+
'keywords': {},24
+
'maxlen': 0,25
+
'option': 0,26
+
'valueType': 'int'},27
+
'FEED_ID': {'comment': 'Feed id',28
+
'dataManagerGroup': 'StandardStMan',29
+
'dataManagerType': 'StandardStMan',30
+
'keywords': {},31
+
'maxlen': 0,32
+
'option': 0,33
+
'valueType': 'int'},34
+
'INTERVAL': {'comment': 'Interval for which this set of parameters is accurate',35
+
'dataManagerGroup': 'StandardStMan',36
+
'dataManagerType': 'StandardStMan',37
+
'keywords': {'QuantumUnits': np.array(['s'], dtype='<U1')},38
+
'maxlen': 0,39
+
'option': 0,40
+
'valueType': 'double'},41
+
'SPECTRAL_WINDOW_ID': {'comment': 'ID for this spectral window setup',42
+
'dataManagerGroup': 'StandardStMan',43
+
'dataManagerType': 'StandardStMan',44
+
'keywords': {},45
+
'maxlen': 0,46
+
'option': 0,47
+
'valueType': 'int'},48
+
'TIME': {'comment': 'Midpoint of time for which this set of parameters is accurate',49
+
'dataManagerGroup': 'StandardStMan',50
+
'dataManagerType': 'StandardStMan',51
+
'keywords': {'MEASINFO': {'Ref': 'UTC', 'type': 'epoch'},52
+
'QuantumUnits': np.array(['s'], dtype='<U1')},53
+
'maxlen': 0,54
+
'option': 0,55
+
'valueType': 'double'},56
+
'TSYS': {'comment': 'System temp. for each of the two receptors',57
+
'dataManagerGroup': 'StandardStMan',58
+
'dataManagerType': 'StandardStMan',59
+
'keywords': {'QuantumUnits': np.array(['K'], dtype='<U1')},60
+
'maxlen': 0,61
+
'ndim': -1,62
+
'option': 0,63
+
'valueType': 'float'},64
+
'_define_hypercolumn_': {},65
+
'_keywords_': {},66
+
'_private_keywords_': {}}67
+
dminfo = {'*1': {'COLUMNS': np.array(['ANTENNA_ID', 'FEED_ID', 'INTERVAL', 'SPECTRAL_WINDOW_ID', 'TIME',68
+
'TSYS'], dtype='<U18'),69
+
'NAME': 'StandardStMan',70
+
'SEQNR': 0,71
+
'SPEC': {'BUCKETSIZE': 1152,72
+
'IndexLength': 7486,73
+
'MaxCacheSize': 2,74
+
'PERSCACHESIZE': 2},75
+
'TYPE': 'StandardStMan'}}76
+
desc_gc = {'ANTENNA_ID': {'comment': 'Antenna identifier',77
+
'dataManagerGroup': 'StandardStMan',78
+
'dataManagerType': 'StandardStMan',79
+
'keywords': {},80
+
'maxlen': 0,81
+
'option': 0,82
+
'valueType': 'int'},83
+
'FEED_ID': {'comment': 'Feed identifier',84
+
'dataManagerGroup': 'StandardStMan',85
+
'dataManagerType': 'StandardStMan',86
+
'keywords': {},87
+
'maxlen': 0,88
+
'option': 0,89
+
'valueType': 'int'},90
+
'GAIN': {'comment': 'Gain polynomial',91
+
'dataManagerGroup': 'StandardStMan',92
+
'dataManagerType': 'StandardStMan',93
+
'keywords': {},94
+
'maxlen': 0,95
+
'ndim': -1,96
+
'option': 0,97
+
'valueType': 'float'},98
+
'INTERVAL': {'comment': 'Interval for which this set of parameters is accurate',99
+
'dataManagerGroup': 'StandardStMan',100
+
'dataManagerType': 'StandardStMan',101
+
'keywords': {'QuantumUnits': np.array(['s'], dtype='<U1')},102
+
'maxlen': 0,103
+
'option': 0,104
+
'valueType': 'double'},105
+
'NUM_POLY': {'comment': 'Number of terms in polynomial',106
+
'dataManagerGroup': 'StandardStMan',107
+
'dataManagerType': 'StandardStMan',108
+
'keywords': {},109
+
'maxlen': 0,110
+
'option': 0,111
+
'valueType': 'int'},112
+
'SENSITIVITY': {'comment': 'Antenna sensitivity',113
+
'dataManagerGroup': 'StandardStMan',114
+
'dataManagerType': 'StandardStMan',115
+
'keywords': {'QuantumUnits': np.array(['K/Jy'], dtype='<U4')},116
+
'maxlen': 0,117
+
'ndim': -1,118
+
'option': 0,119
+
'valueType': 'float'},120
+
'SPECTRAL_WINDOW_ID': {'comment': 'Spectral window identifier',121
+
'dataManagerGroup': 'StandardStMan',122
+
'dataManagerType': 'StandardStMan',123
+
'keywords': {},124
+
'maxlen': 0,125
+
'option': 0,126
+
'valueType': 'int'},127
+
'TIME': {'comment': 'Midpoint of time for which this set of parameters is accurate',128
+
'dataManagerGroup': 'StandardStMan',129
+
'dataManagerType': 'StandardStMan',130
+
'keywords': {'MEASINFO': {'Ref': 'UTC', 'type': 'epoch'},131
+
'QuantumUnits': np.array(['s'], dtype='<U1')},132
+
'maxlen': 0,133
+
'option': 0,134
+
'valueType': 'double'},135
+
'TYPE': {'comment': 'Gain curve type',136
+
'dataManagerGroup': 'StandardStMan',137
+
'dataManagerType': 'StandardStMan',138
+
'keywords': {},139
+
'maxlen': 0,140
+
'option': 0,141
+
'valueType': 'string'},142
+
'_define_hypercolumn_': {},143
+
'_keywords_': {},144
+
'_private_keywords_': {}}145
+
dminfo_gc = {'*1': {'COLUMNS': np.array(['ANTENNA_ID', 'FEED_ID', 'GAIN', 'INTERVAL', 'NUM_POLY',146
+
'SENSITIVITY', 'SPECTRAL_WINDOW_ID', 'TIME', 'TYPE'], dtype='<U18'),147
+
'NAME': 'StandardStMan',148
+
'SEQNR': 0,149
+
'SPEC': {'BUCKETSIZE': 1920,150
+
'IndexLength': 142,151
+
'MaxCacheSize': 2,152
+
'PERSCACHESIZE': 2},153
+
'TYPE': 'StandardStMan'}}154
+
155
+
#################################################################################################156
+
# Main task157
+
def appendantab(vis=None, outvis=None, antab=None):158
+
159
+
# get info from vis160
+
msmd.open(vis)161
+
ant_names = msmd.antennastations()162
+
n_band = len(msmd.bandwidths())163
+
spws = msmd.spwfordatadesc()164
+
165
+
# is gc_begin even needed?166
+
gc_begin = msmd.timerangeforobs(0)['begin']['m0']['value'] * 86400167
+
gc_end = msmd.timerangeforobs(0)['end']['m0']['value'] * 86400168
+
gc_interval = gc_end - gc_begin169
+
gc_time = (gc_begin + gc_end) / 2170
+
msmd.close()171
+
172
+
# get start time and end time from ms173
+
tb.open(vis)174
+
times = tb.getcol('TIME')175
+
first_time = times[0]176
+
last_time = times[-1]177
+
tb.close()178
+
179
+
# remove existing data copy180
+
if os.path.exists(outvis):181
+
print("REMOVING EXISTING COPY")182
+
shutil.rmtree(outvis)183
+
184
+
# run the antab parsing and table filling185
+
shutil.copytree(vis, outvis)186
+
antab_interp(vis, outvis, antab, ant_names, n_band, spws,187
+
gc_interval, gc_time, first_time, last_time)188
+
189
+
190
+
191
+
#################################################################################################192
+
# scanner regex from casavlbitools193
+
def s_err(scanner, error):194
+
raise RuntimeError("line: %d" % scanner.line_no)195
+
196
+
def s_keyword(scanner, token):197
+
if len(token) > 9 or '.' in token:198
+
res = ('value', token)199
+
else:200
+
res = ('key', token.upper())201
+
return res202
+
203
+
def s_newline(scanner, token):204
+
scanner.line_no += 1205
+
return None206
+
207
+
def s_quote(scanner, token):208
+
return ('quote', token[1:-1])209
+
210
+
def s_number(scanner, token):211
+
return ('number', float(token))212
+
213
+
def s_angle(scanner, token): # also time214
+
l = token.split(":")215
+
## It's not python to use reduce.216
+
## But I neither remember nor care what is.217
+
val = reduce(lambda acc, x: 60*acc+math.copysign(acc, float(x)), l, 0.0) 218
+
return ('number', val)219
+
220
+
def s_value(scanner, token):221
+
return ('value', token)222
+
223
+
def s_string(scanner, token):224
+
return ('value', token)225
+
226
+
def s_misc(scanner, token):227
+
logging.debug("Misc token %s at line %d" % (str(token), scanner.line_no))228
+
return ('misc', token)229
+
230
+
def s_comment(scanner, token): # was only used for debugging.231
+
scanner.line_no += 1232
+
return ('comment', token)233
+
234
+
scanner = re.Scanner([235
+
("\!.*\n", s_newline),236
+
("[ \t]+", None),237
+
("\n", s_newline),238
+
("\r\n", s_newline),239
+
("=", lambda s, t: ('equal',)),240
+
("'[^'\n]*'", s_quote), 241
+
("\"[^'\n]*\"", s_quote), # Sigh. Double quotes used in freq.dat242
+
("/", lambda s, t: ('end_chunk',)),243
+
(",", lambda s, t: ('comma',)),244
+
("[+-]?[0-9]+:[0-9]+:[0-9]+(.[0-9]*)?", s_angle),245
+
("[+-]?[0-9]+:[0-9]+(.[0-9]*)?", s_angle),246
+
("[+-]?[0-9]*\.[0-9]+([Ee][+-][0-9]{1,3})?(?![A-Za-z_0-9()])", s_number),247
+
("[+-]?[0-9]+\.?(?![A-Za-z_0-9()])", s_number),248
+
## apparently parens and unquoted minus signs are allowed in keywords?249
+
("[A-Za-z.0-9]([()A-Za-z_0-9._+-]+)?", s_keyword), 250
+
(".*", s_misc)251
+
])252
+
253
+
def p_all():254
+
global tok255
+
chunks = []256
+
try:257
+
while True:258
+
if tok=='EOF':break259
+
chunks.append(p_chunk())260
+
except StopIteration:261
+
pass262
+
return chunks263
+
264
+
def p_chunk():265
+
global tok266
+
entries = []267
+
while tok[0] != 'end_chunk':268
+
entries.append(p_item())269
+
logging.debug("p_chunk %s", str(tok))270
+
tok = next(tokIt)271
+
return entries272
+
273
+
def p_item():274
+
global tok275
+
lhs = p_key()276
+
#if tok[0] == 'key':277
+
#print(tok)278
+
if tok==('equal',):279
+
logging.debug("p_item %s", str(tok))280
+
tok = next(tokIt)281
+
rhs = p_rhs()282
+
elif tok[0] in ['value', 'quote', 'number']:283
+
rhs = p_rhs()284
+
else:285
+
rhs = True # for unitary expressions.286
+
return (lhs, rhs)287
+
288
+
def p_key():289
+
global tok290
+
logging.debug("p_key: %s", str(tok))291
+
if tok[0] == 'key':292
+
res = tok293
+
tok = next(tokIt)294
+
else:295
+
raise RuntimeError("Expected key token, got %s" % str(tok))296
+
return res[1]297
+
298
+
def p_rhs():299
+
global tok300
+
val = p_value()301
+
rhs = [val]302
+
while tok == ('comma',):303
+
logging.debug("p_rhs: %s", str(tok))304
+
tok = next(tokIt)305
+
rhs.append(p_value()) # p_value advances tok beyond the value.306
+
if len(rhs)==1:307
+
rhs = rhs[0]308
+
return rhs309
+
310
+
def p_value():311
+
global tok312
+
if tok[0] not in ['value', 'quote', 'number', 'key']:313
+
raise RuntimeError("Unexpected RHS token %s" % str(tok))314
+
val = tok315
+
logging.debug("p_value: %s", str(val))316
+
tok = next(tokIt)317
+
return val[1]318
+
319
+
scanner.line_no = 0320
+
#################################################################################################321
+
# data table Classes322
+
323
+
class TSysTable:324
+
def __init__(self):325
+
self.time = []326
+
self.time_interval = []327
+
self.source_id = []328
+
self.antenna_no = []329
+
self.array = []330
+
self.freqid = []331
+
self.tsys_1 = []332
+
self.tsys_2 = []333
+
self.tant = []334
+
return335
+
336
+
class GainCurveTable:337
+
def __init__(self):338
+
self.antenna_no = []339
+
self.array = []340
+
self.freqid = []341
+
self.nterm = []342
+
self.y_typ = []343
+
self.gain = []344
+
self.sens_1 = []345
+
self.sens_2 = []346
+
return347
+
348
+
##################################################################349
+
350
+
351
+
def read_keyfile(f):352
+
global tok, tokIt353
+
scanner.line_no = 0354
+
355
+
try:356
+
res = scanner.scan(f.read())357
+
#print(res[0])358
+
if res[1]!='':359
+
raise RuntimeError("Unparsed text: %s." % (res[1][:20]))360
+
tokIt = iter(res[0]+['EOF'])361
+
try: 362
+
tok = next(tokIt)363
+
res = p_all()364
+
except StopIteration: # empty file365
+
res = ''366
+
except RuntimeError as txt:367
+
#print("line %d: %s" % (scanner.line_no, txt), file=sys.stderr)368
+
raise RuntimeError369
+
return res370
+
371
+
def update_map(pols, spws, spwmap, index):372
+
idx = 0373
+
if not isinstance(index, (list, tuple)):374
+
index = [index]375
+
pass376
+
for labels in index:377
+
for label in labels.split('|'):378
+
pol = label[0]379
+
rng = label[1:].split(':')380
+
if pol != 'X':381
+
if not pol in pols:382
+
pols.append(pol)383
+
pass384
+
if len(rng) == 1:385
+
rng.append(rng[0])386
+
pass387
+
rng = [int(x) - 1 for x in rng]388
+
for spw in range(rng[0], rng[1] + 1):389
+
if not spw in spws:390
+
spws.append(spw)391
+
pass392
+
spwmap[(pol, spw)] = idx393
+
continue394
+
pass395
+
continue396
+
idx += 1397
+
continue398
+
spws = sorted(spws)399
+
return400
+
401
+
def find_antenna(keys, ignore):402
+
for key in keys[1:]:403
+
if not type(key[1]) is bool:404
+
continue405
+
if key[0] in ignore:406
+
continue407
+
return key[0]408
+
return None409
+
410
+
def skip_values(infp):411
+
for line in infp:412
+
if line.startswith('!'):413
+
continue414
+
if line.strip().endswith('/'):415
+
break416
+
continue417
+
return418
+
419
+
def get_antenna_index(antenna_name, ant_names):420
+
return ant_names.index(antenna_name)421
+
422
+
def mjd_to_date(mjd):423
+
# Method used in FITSDateUtil424
+
mjd_epoch = dt.datetime(1858, 11 , 17, 0, 0, 0)425
+
426
+
mjd_int = int(mjd / 86400)427
+
delta = dt.timedelta(days=mjd_int)428
+
return mjd_epoch + delta429
+
430
+
def mjd_seconds(yy, mm, dd, d=0):431
+
if (mm < 3):432
+
yy -= 1433
+
mm += 12434
+
435
+
dd += d436
+
b = 0437
+
438
+
if (yy>1582 or (yy==1582 and (mm>10 or (mm==10 and dd >= 15)))): 439
+
b = np.floor(yy/100.)440
+
b = 2 - b + int(b/4)441
+
442
+
val = np.floor(365.25*yy) + np.floor(30.6001*(mm+1)) + dd - 679006.0 + b443
+
return val444
+
445
+
446
+
def get_timetuple(ts):447
+
# ts as string with these possible formats:448
+
# hh.hh 449
+
# hh:mm.mm 450
+
# hh:mm:ss.ss451
+
# NOTE: Regexs below will match any number of decimals on the last quantity (e.g. 19.8222222 and 19.8 both work)452
+
if re.match(r"[0-9]{2}\.[0-9]+", ts):453
+
# hh.hh 454
+
tm_hour = int(ts.split('.')[0])455
+
tm_min = math.modf(60*float(ts.split('.')[1]))456
+
tm_sec = int(60 * tm_min[0])457
+
tm_min = int(tm_min[1])458
+
elif re.match(r"[0-9]{2}:[0-9]{2}\.[0-9]+", ts):459
+
# hh:mm.mm 460
+
tm_hour = int(ts.split(':')[0])461
+
tm_min = math.modf(float(ts.split(':')[1]))462
+
tm_sec = int(60 * tm_min[0])463
+
tm_min = int(tm_min[1])464
+
elif re.match(r"[0-9]{2}:[0-9]{2}:[0-9]{2}$", ts):465
+
# hh:mm:ss466
+
tm_hour = int(ts.split(':')[0])467
+
tm_min = int(ts.split(':')[1])468
+
tm_sec = int(ts.split(':')[2])469
+
elif re.match(r"[0-9]{2}:[0-9]{2}:[0-9]{2}\.[0-9]+", ts):470
+
# hh:mm:ss.ss471
+
tm_hour = int(ts.split(':')[0])472
+
tm_min = int(ts.split(':')[1])473
+
tm_sec = float(ts.split(':')[2])474
+
return tm_hour, tm_min, tm_sec475
+
476
+
def process(infp, keys, pols, data, n_band, spws, first_time, last_time, ant_names):477
+
# Get the antenna name478
+
antenna_name = find_antenna(keys[0], ['SRC/SYS'])479
+
# Skip if no antenna name found480
+
if not antenna_name:481
+
print('ANTENNA missing from TSYS group')482
+
skip_values(infp)483
+
return484
+
try:485
+
antenna = get_antenna_index(antenna_name, ant_names)486
+
except:487
+
print('Antenna %s not present in FITS-IDI file' % antenna_name)488
+
skip_values(infp)489
+
return490
+
491
+
keys = dict(keys[0])492
+
spws = []493
+
spwmap = {}494
+
495
+
update_map(pols, spws, spwmap, keys['INDEX'])496
+
if 'INDEX2' in keys:497
+
update_map(pols, spws, spwmap, keys['INDEX2'])498
+
pass499
+
if len(spws) != n_band:500
+
print('INDEX for antenna %s does not match FITS-IDI file'501
+
% antenna_name, file=sys.stderr)502
+
pass503
+
504
+
spws = range(n_band)505
+
timeoff = 0506
+
if 'TIMEOFF' in keys:507
+
timeoff = float(keys['TIMEOFF'])508
+
509
+
for line in infp:510
+
if line.startswith('!'):511
+
continue512
+
513
+
fields = line.split()514
+
515
+
if len(fields) > 1:516
+
tm_yday = int(fields[0])517
+
# Get timestamp data depending on data format518
+
obs_date = mjd_to_date(first_time)519
+
tm_year = obs_date.year520
+
idi_time = tm.strptime(f"{obs_date.year}-{obs_date.month}-{obs_date.day}","%Y-%m-%d")521
+
idi_ref = tm.mktime(idi_time)522
+
523
+
# Get the start time in sec for MJD using the method from the table filler524
+
mjd_sec = mjd_seconds(int(obs_date.year), int(obs_date.month), int(obs_date.day)) * 86400525
+
526
+
tm_hour, tm_min, tm_sec = get_timetuple(fields[1])527
+
days = (tm_hour*3600 + tm_min*60 + tm_sec) / 86400528
+
t = "%dy%03dd%02dh%02dm%02ds" % \529
+
(tm_year, tm_yday, tm_hour, tm_min, tm_sec)530
+
t = tm.mktime(tm.strptime(t, "%Yy%jd%Hh%Mm%Ss"))531
+
days = (t + timeoff - idi_ref) / 86400532
+
533
+
curr_time = mjd_sec + (days*86400)534
+
535
+
source = 0536
+
if (days) > ((last_time-mjd_sec)/86400) or (days) < ((first_time-mjd_sec)/86400):537
+
source = -1538
+
539
+
#print(deg - mjd_sec)540
+
values = fields[2:]541
+
tsys = {'R': [], 'L': []}542
+
for spw in spws:543
+
for pol in ['R', 'L']:544
+
try:545
+
value = float(values[spwmap[(pol, spw)]])546
+
if value == 999.9:547
+
value = float(-999.9)548
+
pass549
+
except:550
+
value = float(-999.9)551
+
pass552
+
tsys[pol].append(value)553
+
continue554
+
continue555
+
if source != -1:556
+
time_val = mjd_sec + (days * 86400)557
+
data.time.append(time_val)558
+
data.time_interval.append(0.0)559
+
data.antenna_no.append(antenna)560
+
data.tsys_1.append(tsys['R'])561
+
data.tsys_2.append(tsys['L'])562
+
563
+
pass564
+
pass565
+
566
+
if line.strip().endswith('/'):567
+
break568
+
continue569
+
return570
+
571
+
gain_keys = [ 'EQUAT', 'ALTAZ', 'ELEV', 'GCNRAO', 'TABLE', 'RCP', 'LCP' ]572
+
573
+
def process_gc(infp, keys, pols, data, ant_names):574
+
antenna_name = find_antenna(keys[0], gain_keys)575
+
if not antenna_name:576
+
print('Antenna missing from GAIN group')577
+
skip_values(infp)578
+
return579
+
try:580
+
antenna = get_antenna_index(antenna_name, ant_names)581
+
except:582
+
print('Antenna %s not present in FITS-IDI file' % antenna_name)583
+
skip_values(infp)584
+
return585
+
keys = dict(keys[0])586
+
587
+
dpfu = {}588
+
try:589
+
dpfu['R'] = keys['DPFU'][0]590
+
dpfu['L'] = keys['DPFU'][1]591
+
pols.append('R')592
+
pols.append('L')593
+
except:594
+
dpfu['R'] = dpfu['L'] = keys['DPFU']595
+
pols.append('X')596
+
pass597
+
try:598
+
value = keys['POLY'][0]599
+
except:600
+
keys['POLY'] = [keys['POLY']]601
+
pass602
+
603
+
y_typ = 0604
+
if 'ELEV' in keys:605
+
y_typ = 1606
+
elif 'EQUAT' in keys:607
+
y_typ = 1608
+
elif 'ALTAZ' in keys:609
+
y_typ = 2610
+
else:611
+
print('Unknown gain curve type for antenna %s' % antenna_name)612
+
return613
+
614
+
poly = keys['POLY']615
+
data.antenna_no.append(antenna)616
+
data.array.append(1)617
+
data.freqid.append(1)618
+
data.y_typ.append(y_typ)619
+
data.nterm.append(len(poly))620
+
data.gain.append(poly)621
+
data.sens_1.append(dpfu['R'])622
+
data.sens_2.append(dpfu['L'])623
+
return624
+
625
+
def antab_interp(vis, outvis, antab, ant_names, n_band, spws,626
+
gc_interval, gc_time, first_time, last_time):627
+
pols = []628
+
data = TSysTable()629
+
data_gc = GainCurveTable()630
+
keys = StringIO()631
+
fp = open(antab, 'r')632
+
633
+
for line in fp:634
+
if line.startswith('!'):635
+
continue636
+
keys.write(line)637
+
if line.strip().endswith('/'):638
+
keys.seek(0)639
+
try:640
+
tsys = read_keyfile(keys)641
+
except RuntimeError:642
+
pass643
+
if tsys and tsys[0] and tsys[0][0][0] == 'TSYS':644
+
process(fp, tsys, pols, data, n_band, spws, first_time, last_time, ant_names)645
+
pass646
+
elif tsys and tsys[0] and tsys[0][0][0] == 'GAIN':647
+
process_gc(fp, tsys, pols, data_gc, ant_names)648
+
keys = StringIO()649
+
continue650
+
continue651
+
652
+
# Create the subtable for syscal653
+
create_syscal_subtable(data, outvis, spws)654
+
create_gaincurve_subtable(data_gc, outvis, spws, gc_interval, gc_time)655
+
656
+
def create_gaincurve_subtable(data, outvis, spws, gc_interval, gc_time):657
+
# Create the subtable for gain curve658
+
tb.create(f'{outvis}/GAIN_CURVE', desc_gc, dminfo=dminfo_gc)659
+
tb.putkeyword('GAIN_CURVE', f'Table: {outvis}/GAIN_CURVE')660
+
661
+
# assemble columns662
+
ANTENNA_ID = []663
+
FEED_ID = []664
+
SPECTRAL_WINDOW_ID = []665
+
TIME = []666
+
INTERVAL = []667
+
TYPE = []668
+
NUM_POLY = []669
+
GAIN = []670
+
SENSITIVITY = [[],[]]671
+
672
+
# Fill the arrays673
+
for i in range(len(data.antenna_no)):674
+
for j in spws:675
+
TIME.append(gc_time)676
+
ANTENNA_ID.append(data.antenna_no[i])677
+
FEED_ID.append(-1)678
+
INTERVAL.append(gc_interval)679
+
TYPE.append(str(data.y_typ[i]))680
+
NUM_POLY.append(data.nterm[i])681
+
# ASK GEORGE ABOUT THIS. SHAPE OF EACH ELEMENT VARIES682
+
# Might be an existing bug...683
+
GAIN.append(data.gain[i])684
+
SPECTRAL_WINDOW_ID.append(int(j))685
+
SENSITIVITY[0].append(data.sens_1[i])686
+
SENSITIVITY[1].append(data.sens_2[i])687
+
688
+
# Add the columns to the table689
+
tb.open(f'{outvis}/GAIN_CURVE', nomodify=False)690
+
tb.addrows(len(ANTENNA_ID))691
+
tb.putcol('TIME', TIME)692
+
tb.putcol('ANTENNA_ID', ANTENNA_ID)693
+
tb.putcol('FEED_ID', FEED_ID)694
+
tb.putcol('INTERVAL', INTERVAL)695
+
tb.putcol('TYPE', TYPE)696
+
tb.putcol('NUM_POLY', NUM_POLY)697
+
#Fill gain row by row698
+
for i in range(len(GAIN)):699
+
# If the gain col shape doesn't match what the table filler does700
+
if len(np.shape(GAIN[i])):701
+
tb.putcell('GAIN', i, [GAIN[i],GAIN[i]])702
+
else:703
+
tb.putcell('GAIN', i, GAIN[i])704
+
tb.putcol('SPECTRAL_WINDOW_ID', SPECTRAL_WINDOW_ID)705
+
tb.putcol('SENSITIVITY', SENSITIVITY)706
+
#tb.flush()707
+
tb.close()708
+
709
+
710
+
def create_syscal_subtable(data, outvis, spws):711
+
# Create the subtable for syscal712
+
tb.create(f'{outvis}/SYSCAL', desc, dminfo=dminfo)713
+
tb.putkeyword('SYSCAL', f'Table: {outvis}/SYSCAL')714
+
715
+
# Assemble columns716
+
ANTENNA_ID = []717
+
FEED_ID = []718
+
INTERVAL = []719
+
SPW_ID = []720
+
TIME = []721
+
TSYS = [[],[]]722
+
723
+
# Fill the arrays724
+
for i in range(len(data.time)):725
+
for j in spws:726
+
ANTENNA_ID.append(data.antenna_no[i])727
+
FEED_ID.append(0)728
+
INTERVAL.append(data.time_interval[i])729
+
SPW_ID.append(int(j))730
+
TIME.append(data.time[i])731
+
#TIME.append(ref_time + data.time[i]*86400)732
+
# Still need to figure out tsys...733
+
idx = int(i/len(spws))734
+
TSYS[0].append(data.tsys_1[idx][j])735
+
TSYS[1].append(data.tsys_2[idx][j])736
+
737
+
# Add the columns to the table738
+
tb.open(f'{outvis}/SYSCAL', nomodify=False)739
+
tb.addrows(len(TIME))740
+
tb.putcol('TIME', TIME)741
+
tb.putcol('ANTENNA_ID', ANTENNA_ID)742
+
tb.putcol('INTERVAL', INTERVAL)743
+
tb.putcol('SPECTRAL_WINDOW_ID', SPW_ID)744
+
tb.putcol('TSYS', TSYS)745
+
#tb.flush()746
+
tb.close()