-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy path3C277_full.html
2021 lines (1966 loc) · 151 KB
/
3C277_full.html
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
<!DOCTYPE HTML>
<!--
Hyperspace by HTML5 UP
html5up.net | @ajlkn
Free for personal and commercial use under the CCA 3.0 license (html5up.net/license)
-->
<html>
<head>
<title>ERIS2024 - 3C277.1 full tutorial</title>
<meta charset="utf-8" />
<meta name="viewport" content="width=device-width, initial-scale=1, user-scalable=no" />
<link rel="stylesheet" href="assets/css/main.css" />
<link rel="stylesheet" href="assets/css/prism.css">
<noscript><link rel="stylesheet" href="assets/css/noscript.css" /></noscript>
<script type="text/javascript" async src="assets/MathJax/MathJax.js?config=TeX-MML-AM_CHTML">
</script>
<link href='https://fonts.googleapis.com/css?family=Inconsolata' rel='stylesheet'>
<script type="text/x-mathjax-config">MathJax.Hub.Config({tex2jax: {inlineMath: [['$','$'], ['\\(','\\)']]}});
</script>
<script type="text/javascript">
var pcs = document.lastModified.split(" ")[0].split("/");
var date = pcs[1] + '/' + pcs[0] + '/' + pcs[2];
onload = function(){
document.getElementById("lastModified").innerHTML = "Page last modified: " + date;
}
</script>
<script>
// When the user clicks on <div>, open the popup
function myFunction(a) {
var popup = document.getElementById(a);
popup.classList.toggle("show");
}
</script>
</head>
<body class="is-preload">
<!-- Header -->
<header id="header">
<a href="index.html" class="title">ERIS 2024 - Tutorials</a>
<nav>
<ul>
<li><a href="index.html">Home</a></li>
<li><a href="intro_data.html">Intro. to data</a></li>
<li><a href="calibration.html">Calibration</a></li>
<li><a href="imaging.html" >Imaging</a></li>
<li><a href="selfcalibration.html">Self calibration</a></li>
<li><a href="adv_imaging.html">Advanced imaging</a></li>
<li><a href="3C277_full.html" class="active">Full tutorial</a></li>
</ul>
</nav>
</header>
<!-- Wrapper -->
<div id="wrapper">
<!-- Main -->
<section id="main" class="wrapper">
<div class="inner">
<h1 class="major">3C277.1 full tutorial</h1>
<!-- Text -->
<section>
<p>
This guide provides a basic example for calibrating e-MERLIN C-Band (5 GHz) continuum data in CASA. It has been tested in CASA 6.5. It is designed to run using an averaged data set <code class="language-none">all_avg.ms</code> (1.7 G), in which case you will need at least 5.5 G disk space to complete this tutorial. This tutorial is based upon the DARA 3C277.1 tutorial but is using new 2019 observations of the sources and some updates for ERIS 2024. If you find bugs in the code please contact Jack Radcliffe ([email protected]).
</p>
<p>
<i>If you have a very slow laptop, you are advised to skip some of the <code class="language-none">plotms</code> steps (see the plots linked to this page instead) and especially, do not use the <code class="language-none">plotms</code> options to write out a png. This is not used in this web page, but if you run the script as provided, you are advised to comment out the calls to <code class="language-none">plotms</code>or at least the part which writes a png, see tutors for guidance.</i>
</p>
</section>
<hr/>
<section>
<h2>Table of Contents</h2>
<p>
<b>Important:</b> The numbers in brackets respond to the steps in the accompanying python scripts.
</p>
<p>
<h3><a href="#intro">Part 1 - Introduction to interferometric data</a></h3>
<ol>
<li><a href="#guidance">Guidance for calibrating 3C277.1 in CASA</a></li>
<ol class="alt">
<li><a href="#data">The data and supporting material</a></li>
<li><a href="#use_guide">How to use this guide</a></li>
</ol>
<li><a href="#data_inspection">Data inspection and flagging</a></li>
<ol class="alt">
<li><a href="#Check_data">Check data: listobs and plotants (1-3)</a></li>
<li><a href="#late_on_source">Identify 'late on source' bad data (4)</a></li>
<li><a href="#flag_start">Flag the bad data at the start of each scan (5)</a></li>
<li><a href="#flag_end_chan">Flag the bad end channels (6-7)</a></li>
<li><a href="#identify_bad_data">Identify and flag remaining bad data (8-10)</a></li>
</ol>
</ol>
<h3><a href="#calibration">Part 2 - Calibration</a></h3>
<ol start='3'>
<li><a href="#flux_scaling">Setting the flux of the primary flux scale calibrator (11)</a></li>
<li><a href="#bandpass">Bandpass calibration</a></li>
<ol class="alt">
<li><a href="#decorrelation">Decorrelation and solution intervals (12)</a></li>
<li><a href="#pre_bandpass">Pre-bandpass calibration (13)</a></li>
<li><a href="#init_bandpass">Derive bandpass solutions (14)</a></li>
</ol>
<li><a href="#delay">Time-dependent delay calibration (15)</a></li>
<ol class="alt">
<li><a href="#inspect_bandpass">Inspect effects (16)</a></li>
</ol>
<li><a href="#gaincal">Gain calibration (17)</a></li>
<ol class="alt">
<li><a href="#time_dep_phase">Time-dependent phase calibration of all calibration sources</a></li>
<li><a href="#time_dep_amp">Time-dependent amplitude calibration of all calibration sources</a></li>
</ol>
<li><a href="#determine_flux">Determining true flux densities of other calibration sources (18-19)</a></li>
<li><a href="#apply_soln_plot">Apply solutions to target (20-21)</a></li>
</ol>
</ol>
<h3><a href="#imaging">Part 3 - Imaging</a></h3>
<ol start="9">
<li><a href="#imaging101">Imaging 101</a></li>
<li><a href="#determining">Determining imaging parameters (1-2)
<li><a href="#first_image">First image of 3C277.1 (2)</a></li>
<li><a href="#image_properties">Measuring image properties (3)</a></li>
</ol>
<h3><a href="#sc">Part 4 - Self calibration</a></h3>
<ol start="13">
<li><a href="#inspect_self_cal_p">Inspect to choose interval for initial phase self-calibration (4-5)</a></li>
<li><a href="#self_cal_derive_p1">Derive phase self-calibration solutions (step 6)</a></li>
<li><a href="#self_cal_apply_p1">Apply first phase solutions and image again (7)</a></li>
<li><a href="#self_cal_p2">Phase-self-calibrate again (8)</a></li>
<li><a href="#apply_im_p">Apply the solutions and re-image (9)</a></li>
<li><a href="#inspect_self_cal_a">Choose solution interval for amplitude self-calibration (10)</a></li>
<li><a href="#self_cal_a">Amplitude self-cal (11)</a></li>
<li><a href="#apply_im_p">Apply amplitude and phase solutions and image again (12)</a></li>
<li><a href="#end_sc">Summary</a></li>
</ol>
<h3><a href="#adv_im">Part 5 - Advanced imaging</a></h3>
<ol start="22">
<li><a href="#reweigh">Reweigh the data according to sensitivity (13)</a></li>
<li><a href="#multiscale">Try multi-scale clean (14)</a></li>
<li><a href="#taper">Taper the visibilities during imaging (15)</a></li>
<li><a href="#briggs">Different imaging weights: Briggs robust (16)</a></li>
<li><a href="#summary">Summary and conclusions</a></li>
</ol>
</p>
</section>
<hr/>
<section>
<h2 class="alt"><a name="intro"></a>Part 1 - Introduction to interferometric data</h2>
<h3 class="alt"><a name="guidance"></a>1. Guidance for calibrating 3C277.1 in CASA</h3>
<h4><a name="data"></a>1A. The data and supporting material</h4>
<p>
For this workshop, we need the following files. Please double check that you have the following in your current directory:
<ul>
<li><code class="language-bash">all_avg.ms</code> - the data (after conversion to a measurement set).</li>
<li><code class="language-bash">3C286_C.clean.model.tt0</code> - used for fluxscaling of the data set.</li>
<li><code class="language-bash">all_avg_all.flags.txt</code> - all flags</li>
<li><code class="language-bash">3C277.1_cal_outline2024.py</code> - the script that we shall input values into to calibrate these data.</li>
<li><code class="language-bash">3C277.1_cal_all2024.py.tar.gz</code> - cheat script for calibating these data.</li>
<li><code class="language-bash">3C277.1_imaging_outline2024.py</code> - the script that we shall input values into to image and self-calibrate the data. </li>
<li><code class="language-bash">3C277.1_imaging_all2024.py.tar.gz</code> - cheat script for imaging and calibrating these data.</li>
<li><code class="language-bash">1252+5634.ms.tar.gz</code> - measurement set if the calibration goes wrong</li>
</ul>
Note that some of these are possibly tarred (i.e., end in <code class="language-none">.tar.gz</code> or <code class="language-none">.tar</code>) so you need to use <code class="language-bash">tar -xzvf <filename></code> to untar these.
</p>
</section>
<hr/>
<section>
<h4><a name="use_guide"></a>1B. How to use this guide</h4>
<p>
This user guide presents inputs for CASA tasks e.g., <code class="language-bash">gaincal</code> to be entered interactively at a terminal prompt for the calibration of the averaged data and initial target imaging. This is useful to examine your options and see how to check for syntax errors. You can run the task directly by typing <code class="language-bash">gaincal</code> at the prompt; if you want to change parameters and run again, if the task has written a calibration table, delete it (make sure you get the right one) before re-running. You can review these by typing:
</p>
<pre><code class="language-python"># In CASA
!more gaincal.last</code></pre>
<p>
This will output something like:
</p>
<pre><code class="language-none">vis = 'all_avg.ms'
caltable = 'bpcals_precal.ap1'
field = '1407+284'
...
#gaincal( vis='all_avg.ms',caltable='bpcals_precal.ap1',field='1407+284',spw='',intent='',selectdata=True,solint='120s',combine='',
preavg=-1.0,refant='Mk2',refantmode='flex',minblperant=2,minsnr=2,solnorm=True,gaintype='G',smodel=[],calmode='ap',solmode='',rmsthresh=[],
corrdepflags=False,append=True,docallib=False,parang=False,timerange='',uvrange='',antenna='',scan='',observation='',msselect='',normtype='mean',
splinetime=3600.0,npointaver=3,phasewrap=180.0,callib='',gaintable=['all_avg.K', 'bpcals_precal.p1'],gainfield=[''],interp=[],spwmap=[] )</code></pre>
<p>
As you can see in the script, the second format (without the <code class="language-python">#</code>) is the form to use in the script, with a comma after every parameter value. When you are happy with the values which you have set for each parameter, enter them in the script <code class="language-bash">3C277.1_cal_outline2024.py</code> (or whatever you have called it). You can leave out parameters for which you have not set values; the defaults will work for these (see e.g., <code class="language-python">help(gaincal)</code> for what these are), but feel free to experiment if time allows.
</p>
<p>
The script already contains most plotting inputs in order to save time. Make and inspect the plots; change the inputs if you want, zoom in etc. There are links to copies of the plots (or parts of plots) but only click on them once you have tried yourself or if you are stuck.
</p>
<p>
The parameters which should be set explicitly or which need explanation are given below, with <code class="language-python">'**'</code> if they need values. Steps that you need to take action on are in bullet points.
</p>
<p>
The data set that we are using today is a 2019 observation of 3C277.1, a bright radio AGN with two-sided jets. Within the data data set, named <code class="language-bash">all_avg.ms</code> it contains four fields that we are going to use to calibrate these data.:
</p>
<p>
<table>
<thead>
<tr>
<th><i>e</i>-MERLIN name</th>
<th>Common name</th>
<th>Role</th>
</tr>
</thead>
<tbody>
<tr>
<td>1252+5634</td>
<td>3C277.1</td>
<td>target</td>
</tr>
<tr>
<td>1302+5748</td>
<td>J1302+5748</td>
<td>phase reference source</td>
</tr>
<tr>
<td>0319+415</td>
<td>3C84</td>
<td>bandpass, polarisation leakage cal. (usually 10-20 Jy at 5-6 GHz)</td>
</tr>
<tr>
<td>1331+3030</td>
<td>3C286</td>
<td>flux, polarisation angle cal. (flux density known accurately, see Step XX)</td>
</tr>
</tbody>
</table>
These data were taken in full polarization with 0.512 GHz bandwidth centred on 5072 MHz, in 4 adjacent spectral windows (spw's, also known as IFs in AIPS). Each 128-MHz spw contained 512 channels, each 0.25 MHz wide, and 1 s integrations were used (the data are later averaged). We will show in the next section how we obtained this information!
</p>
<p>
These data have already been preprocessing from the original raw fits-idi files from the correlator which included the following steps.
<ol>
<li>Conversion from fitsidi to a CASA compatible measurement set.</li>
<li>Sorted and recalculated the $uvw$ (visibility coordinates), concatenate all data and adjust scan numbers to increment with each source change.</li>
<li>Averaged every 8 channels and 4s.</li>
</ol>
</p>
</section>
<hr/>
<section>
<h3 class="alt">2. <a name="data_inspection">Data inspection and flagging</a></h3>
<h4 class="alt">2A. <a name="Check_data">Check data: listobs and plotants (steps 1-3)</a></h4>
<p class="alt"><a href="#"><font size="-5">← back to top</font></a></p>
<p>
In this part, we shall see how we can get the information from our data set. One of the golden rules of data reduction is to <b>know your data</b>, as different arrays, configurations and frequencies require different calibration strategies.
</p>
<p>
<ul>
<li>Firstly, check that you have <code class="language-bash">all_avg.ms</code> in a directory with enough space and start CASA.</li>
<li>Enter the parameter in step 1 of the script to specify the measurement set. The task <code class="language-bash">listobs</code> will summarise the information given above and write it to a file named <code class="language-none">all_avg.listobs.txt</code>.</li>
</ul>
</p>
<pre class="line-numbers" data-start="83"><code class="language-python"> listobs(vis='**',
overwrite=True,
listfile='all_avg.listobs.txt')</code></pre>
<p>
<b>Important:</b> To make it easier, the line numbers corresponding to the script <code class="language-bash">3C277.1_cal_outline2024.py</code> are included in the excerpts from the script shown in this tutorial.
</p>
<p>
<ul>
<li>When you have entered the parameter, we want to execute the step so that we run the task. In CASA you want to enter the following:
<pre><code class="language-python">runsteps=[1];
execfile('3C277.1_cal_outline2024.py')</code></pre></li>
<li>Check the CASA logger to ensure there are no errors (look for the word SEVERE) and check your current directory for a new file called <code class="language-none">all_avg.listobs.txt</code>.</li>
<li>Open this file using your favourite text editor (e.g., gedit, emacs, vim) and inspect it.</li>
</ul>
<p>
Selected entries from this <code class="language-bash">listobs</code> file, ordered by time, show the various scans during the observation,
</p>
<pre><code class="language-none"> Date Timerange (UTC) Scan FldId FieldName nRows SpwIds Average Interval(s) ScanIntent
01-Aug-2019/23:20:10.5 - 00:00:00.2 1 3 0319+4130 35820 [0,1,2,3] [4, 4, 4, 4]
02-Aug-2019/00:00:04.0 - 00:00:07.0 2 2 1302+5748 60 [0,1,2,3] [3, 3, 3, 3]
00:04:54.5 - 00:06:01.0 3 1 1252+5634 848 [0,1,2,3] [4.14, 4.14, 4.14, 4.14]
00:06:04.0 - 00:07:59.5 4 2 1302+5748 1484 [0,1,2,3] [3.96, 3.96, 3.96, 3.96]
...
09:30:03.0 - 10:00:00.2 183 3 0319+4130 12428 [0,1,2,3] [4, 4, 4, 4]
10:00:03.0 - 10:30:00.2 184 0 1331+3030 16980 [0,1,2,3] [4, 4, 4, 4]
10:30:03.0 - 10:32:00.3 185 2 1302+5748 364 [0,1,2,3] [4, 4, 4, 4]
10:32:03.0 - 10:36:01.0 186 1 1252+5634 3220 [0,1,2,3] [4.02, 4.02, 4.02, 4.02]
10:36:03.0 - 10:38:00.3 187 2 1302+5748 1436 [0,1,2,3] [4, 4, 4, 4]
...</code></pre>
<p>
There are 4 fields:
</p>
<pre><code class="language-none">Fields: 4
ID Code Name RA Decl Epoch nRows
0 ACAL 1331+3030 13:31:08.287300 +30.30.32.95900 J2000 46256
1 1252+5634 12:52:26.285900 +56.34.19.48800 J2000 560340
2 1302+5748 13:02:52.465277 +57.48.37.60932 J2000 242296
3 CAL 0319+4130 03:19:48.160110 +41.30.42.10330 J2000 91856</code></pre>
<p>
Four spectral windows (each with 64 channels) with all four polarisation products:
</p>
<pre><code class="language-none">Spectral Windows: (4 unique spectral windows and 1 unique polarization setups)
SpwID Name #Chans Frame Ch0(MHz) ChanWid(kHz) TotBW(kHz) CtrFreq(MHz) Corrs
0 none 64 GEO 4817.000 2000.000 128000.0 4880.0000 RR RL LR LL
1 none 64 GEO 4945.000 2000.000 128000.0 5008.0000 RR RL LR LL
2 none 64 GEO 5073.000 2000.000 128000.0 5136.0000 RR RL LR LL
3 none 64 GEO 5201.000 2000.000 128000.0 5264.0000 RR RL LR LL</code></pre>
<p>
And 6 antennas participating:
</p>
<pre><code class="language-none">Antennas: 6:
ID Name Station Diam. Long. Lat. Offset from array center (m) ITRF Geocentric coordinates (m)
East North Elevation x y z
0 Mk2 Mk2 24.0 m -002.18.14.1 +53.02.57.3 19618.7284 20856.7583 6908.7107 3822846.760000 -153802.280000 5086285.900000
1 Kn Kn 25.0 m -002.59.49.7 +52.36.17.2 -26823.2185 -28465.4973 7055.9694 3860084.921000 -202105.036000 5056568.917000
2 De De 25.0 m -002.08.40.0 +51.54.49.9 30300.7876 -105129.6730 7263.6278 3923442.718000 -146914.386000 5009755.292000
3 Pi Pi 25.0 m -002.26.43.3 +53.06.14.9 10141.4322 26944.5297 6845.6479 3817549.874000 -163031.121000 5089896.642000
4 Da Da 25.0 m -002.32.08.4 +52.58.17.2 4093.0417 12222.9915 6904.6753 3829087.869000 -169568.939000 5081082.350000
5 Cm Cm 32.0 m +000.02.13.7 +51.58.49.3 176453.7144 -97751.3908 7233.2945 3920356.264000 2541.973000 5014284.480000</code></pre>
<p>
In the next step, we are going to plot the positions of the antennae present in our observation. This is a handy way of getting a feel on the expected resolution and will assist when we try to determine a reference antenna.
</p>
<p>
<ul>
<li>Take a look at step 2 and enter one parameter to plot the antenna positions (optionally, a second to write this to a png).<pre class="line-numbers" data-start="93"><code class="language-python"> plotants(vis='**', figfile='**')</code></pre></li>
<li>Execute step 2 by setting <code class='language-python'>runsteps</code> and using <code class='language-python'>execfile</code> like in step 1.</li>
</ul>
</p>
<p>
You should end up with a plot that looks similar to that below.
</p>
<div class="col-12"><span class="image fit"><img src="plots/antpos.png" alt="" /></span></div>
<p>
Consider what would make a good reference antenna. Although Cambridge has the largest diameter, it has no short baselines. We will use Mk2 for the reference antenna. Ideally we want the reference antenna to look through similar atmospheric conditions so we typically choose the one with the shortest baselines, or we pick the most sensitive antenna.
</p>
<p>
Next we want to plot the $uv$ coverage of the data for the phase-reference source so that we can see what Fourier scales our interferometer is sampling and whether we need to adopt our calibration strategy to adapt for snapshot imaging.
</p>
<p>
<ul>
<li>Check out step 3 in the script. You need to enter several parameters but some have been done for you. Check out the annotated <code class="language-bash">listobs</code> output above to try to identify these but remember to also use the <code class="language-python">help</code> function,
<pre class="line-numbers" data-start="100"><code class="language-python"> plotms(vis='**',
xaxis='**',
yaxis='**',
coloraxis='spw',
avgchannel='8',
field='1302+5748',
plotfile='',
showgui=gui, overwrite=True)</code></pre></li>
<li>Execute step 3 and check the new <code class="language-bash">plotms</code> graphical user interface (GUI) that should have appeared on your screen.</li>
</ul>
</p>
<p>
Note that we averaged in channels a bit so that this doesn't destroy your computer. If this takes a long time, move on and look at the plot below. The $u$ and $v$ coordinates represent (wavelength/projected baseline) as the Earth rotates whilst the source is observed. Thus, as each spw is at a different wavelength, it samples a different part of the $uv$ plane, improving the aperture coverage of the source, allowing Multi-Frequency Synthesis (MFS) of continuum sources.
</p>
<div class="col-12"><span class="image fit"><img src="plots/1302+5748_uv.png" alt="" /></span></div>
<p>
With our initial data inspection under way, we are now going to move onto finding and removing bad data that can adversely affect our calibration later down the line.
</p>
</section>
<hr/>
<section>
<h4 class="alt"><a name="late_on_source"></a>2B. Identify 'late on source' bad data (step 4)</h4>
<p class="alt"><a href="#"><font size="-5">← back to top</font></a></p>
<p>
The first check for bad data that we are going to look for is whether the antennae are on-source for the whole duration of the scan. This is important to check, especially for heterogeneous arrays, where the different slewing rates means that different antennae come onto source at different times. Often this is flagged automatically by the observatory but it is not always the case. The flagging of these data is called 'quacking' (don't ask me why!). Let's first plot the data to see if we need to quack the data!
</p>
<p>
<ul>
<li>Check out step 4 and use <code class="language-bash">plotms</code> to plot the phase-reference source amplitudes as a function of time. To save time we want to select just a few central channels and average them. This needs to be enough to give good S/N but not so many that phase errors across the band cause decorrelation. Some hints on the proper parameter selection are given in the code block below.</li>
</ul>
</p>
<pre class="line-numbers" data-start="119"><code class="language-python"> plotms(vis='**',
field='**', #phase reference source
spw='**', #plot just inner 1/2 channels for each spw (check listobs for the # channels)
avgchannel='**', #average these channels together
xaxis='**',yaxis='**', #plot amp versus time
antenna=antref+'&Pi', #uses Mk2 as a reference antenna for comparison but other w/ short baselines are ok
correlation='**,**', #just select the parallel polarisations as polarised intensity is fainted
coloraxis='corr', #colour by correlation
showgui=gui,
overwrite=True, plotfile='')</code></pre>
<p>
<ul>
<li>Once you have inputted the parameters, execute step 4 and the plotms GUI should have opened</li>
<li>Zoom in onto just a few scans and you should see a plot similar to that shown below.</li>
</ul>
</p>
<div class="col-12"><span class="image fit"><img src="plots/pre-quack.png" alt="" /></span></div>
<p>
In this case, the pre-applied observatory flags which removed data when the antennae was off source was not perfect, resulting in the single integration at the start which needs to be flagged. This is roughly the same amount of time at the start of each scan is bad for all sources and antennas, so the phase reference is a good source to examine. (In a few cases there is extra bad data, ignore that for now).
<ul>
<li>Estimate this time interval using the plot and we shall then enter that value in next part.</li>
</ul>
</p>
</section>
<hr/>
<section>
<h4 class="alt"><a name="flag_start"></a>2C. Flag the bad data at the start of each scan (step 5)</h4>
<p class="alt"><a href="#"><font size="-5">← back to top</font></a></p>
<p>
One of the most important rules for good data reduction is good bookkeeping i.e., knowing exactly what has been done to these data, and how to return back to how the data looked in a previous step. Often you may find that you have entered an incorrect parameter and this has messed up your data, such as accidentally flagging all your data. CASA has some inbuilt methods to help with this bookkeeping. One of these is the task <code class='language-python'>flagmanager</code> which can be used to remember all the flags at a certain point in the calibration steps. This means that you can use this task to revert back to how the flags were before you changed them.
</p>
<p>
<ul>
<li>Take a look at step 5, we are going to input parameters for both tasks before executing the step this time.</li>
<li>The <code class="language-python">flagmanager</code> task will be used to <b>save</b> the current state of the flags before we remove that bad data we identified at the start of the scans. Enter the parameters that will save the flags.
<pre class="line-numbers" data-start="136"><code class="language-python"> flagmanager(vis='**',
mode='**',
versionname='pre_quack')</code></pre></li>
</ul>
</p>
<p>
If the parameters are set correctly, then we will generate a backup named <code class="language-bash">pre_quack</code> in a <code class="language-bash">all_avg.ms.flagversions</code> file. After we execute step 5, you can check that it exists by running <code class="language-bash">flagmanager</code> again using <code class="language-python">mode='list'</code> instead.
</p>
<p>
Next, we need to actually do the quacking or flagging of these channels. We use the CASA task <code class="language-bash">flagdata</code> to do this along with the special mode called <code class="language-bash">quack</code>.
</p>
<p>
<ul>
<li>Fill in the <code class="language-bash">flagdata</code> parameters, inserting the time that the antennas were off source that we found earlier.
<pre class="line-numbers" data-start="141"><code class="language-python"> flagdata(vis='**',
mode='**',
quackinterval='**') #needs to be a integer</code></pre></li>
<li>Execute step 5. The script will also re-run the <code class="language-python">plotms</code> command from step 4 to show that the data is removed (if your GUI is still open we could tick the reload option and replot it.</li>
</ul>
</p>
<p>
This step should produce a plot similar to that below (when you zoom in) which shows the bad data removed.
</p>
<div class="col-12"><span class="image fit"><img src="plots/quacked.png" alt="" /></span></div>
</section>
<hr/>
<section>
<h4 class="alt"><a name="#flag_end_chan"></a>2D. Flag the bad end channels (steps 6-7)</h4>
<p class="alt"><a href="#"><font size="-5">← back to top</font></a></p>
<p>
Due to the receiver engineering and filters on the receivers, the edges of each spw in the frequency domain are of lower sensitivity and are normally be flagged. We often correct for this effect (called bandpass calibration), but the low S/N data on these edge channels can result in errant solutions and noise being injected into these data. This effect is instrumental and so is the same for all sources, polarizations and antennas. However, different spectral windows may be different. We need to look at a very bright source which has enough signal in each channel in order to see this effect. In this case we can use the source 0319+415.
</p>
<p>
To obtain the correct parameters for the flagging part, we want to identify the channels where the amplitudes versus frequencies are low so we want to plot these data.
</p>
<p>
<ul>
<li>Look at step 6. For the <code class="language-bash">plotms</code> command below. Add in the correct parameters to plot the frequency versus time for 0319+415.
<pre class="line-numbers" data-start="164"><code class="language-python"> plotms(vis='**',
field='**',
xaxis='**',yaxis='**', #plot amplitude against channel
antenna=antref+'&Kn', #select one baseline
correlation='**,**', #just parallel hands
avgtime='**', #set to a very large number
avgscan='**', #average all scans
iteraxis='**', coloraxis='**', #view each spectral window in turn, and colour
showgui=gui,
overwrite=True, plotfile='')</code></pre></li>
</ul>
</p>
<p>
The use of <code class="language-bash">iteraxis</code> allows us to scroll between the spectral windows (using the green arrow buttons in the GUI). Below is the plot for the first spectral window.
</p>
<div class="col-12"><span class="image fit"><img src="plots/identify_end_chans_Spw0.png" alt="" /></span></div>
<p>
From looking at this plot we are going to flag the channels dropping off to zero at the edge. These are going to be approximately the first 5 channels (0-4) and the last three channels (61-63).
</p>
<p>
<ul>
<li>Do the same for the other three spectral windows and note down which channels to flag (remember that we are using zero indexing so the first channel is 0).</li>
</ul>
</p>
<p>
With the channels for all spectral windows identified, we want to now actually flag them. However, remember that we should always back-up our flags should things go wrong.
</p>
<p>
<ul>
<li>Look at step 7. Back up the flags and enter the appropriate parameters in <code class="language-bash">flagdata</code>. Note the odd CASA notation for selecting spectral windows and channels. We have provided the "bad"/insensitive channels in spw 0 with the correct syntax, so use that as a guide for the other spectral windows. You can add multiple spectral windows in the same string using a comma e.g., <code class="language-bash">spw='0:0~3,1:0~4'</code>, will flag channels 0 to 3 in spectral window 0 and 0 to 4 in spectral window 1.
<pre class="line-numbers"><code class="language-python"> flagmanager(vis='all_avg.ms',
mode='save',
versionname='pre_endchans')
# End chans
flagdata(vis='all_avg.ms',
mode='manual',
spw='0:0~4;61~63,***', #enter the ranges to be flagged for the other spw
flagbackup=False)</code></pre></li>
</ul>
</p>
<p>
<ul>
<li>Execute step 7 and then re-plot in <code class="language-bash">plotms</code> (the code should already do this replotting for you!). You should get a plot similar to that below</li>
</ul>
</p>
<div class="col-12"><span class="image fit"><img src="plots/0319+4130_flagged_avg_amp-chan.png" alt="" /></span></div>
</section>
<hr/>
<section>
<h4 class="alt"><a name="identify_bad_data"></a>2E. Identify and flag remaining bad data (steps 8-10)</h4>
<p class="alt"><a href="#"><font size="-5">← back to top</font></a></p>
<p>
The final step in the flagging process is to closely inspect your data and look for other bad data. This can take many forms such as correlator issues (which causes low amplitudes on some baselines), antenna problems such as warm loads and Radio Frequency Interference (RFI) causing huge spikes in the data (plus many many more).
</p>
<p>
You usually have to go through all baselines, averaging channels and plotting amplitude against time and frequency. This can be very tedious so there has been considerable efforts to produce automatic flaggers (e.g., AOFlagger, rflag) which can do this automatically. However, you still want to judge whether these algorithms are working properly and not flagging decent data. Therefore having some experience doing it manually can give you experience on what exactly we consider as 'good' and 'bad' data.
</p>
<p>
To save time (and your sanity!), we are not going to do it all here. We are just going to look at one source, 0319+4130, and write down a few commands. This is done in the form of a list file, which has no commas and the only spaces must be a single space between each parameter for a given flagging command line. Some errors affect all the baselines to a given antenna; others (e.g., correlator problems) may only affect one baseline. Usually all spw will be affected so you can just inspect one but check afterwards that all are 'clean'. It can be hard to see what is good or bad for a faint target but if the phase-reference source is bad for longer than a single scan then usually the target will be bad for that intervening time.
</p>
<p>
Before we begin to flag, remember we should start by backing up the flags again!
</p>
<p>
<ul>
<li>Take a look at step 8 and enter the correct parameters for <code class="language-python">flagmanager</code><pre class="line-numbers" data-start='216'><code class="language-python"> flagmanager(vis='**',
mode='**',
versionname='pre_timeflagging')</code></pre></li>
<li>Take a look at the <code class="language-python">plotms</code> commands in the file. These will plot all the sources and save them to a file. The final command will put the amplitude versus time for 0319+4130, the source we want to flag now.</li>
<li>Use the <code class="language-bash">plotms</code> window to inspect 0319+4130, paging (iterating) through the baselines.</li>
<li>Zoom in and identify the time ranges and antennae where we have bad data. Insert these into a text file and save that text file as <code class="language-bash">0319+4130_flags.txt</code>. Some example commands are given below and a plot showing some bad data on 0319+4130 after zooming in is shown</li>
</ul>
</p>
<div class="col-12"><span class="image fit"><img src="plots/example_flagging.png" alt="" /></span></div>
<p>
Here are some of the flagging commands in list format for the data shown below,
</p>
<pre><code class="language-none">mode='manual' field='0319+4130' timerange='2019/08/02/09:30:40~2019/08/02/09:50:36'
mode='manual' field='0319+4130' timerange='2019/08/02/15:00:00~2019/08/02/15:05:40'
mode='manual' field='0319+4130' antenna='Cm' timerange='2019/08/02/15:18:18~2019/08/02/15:18:30'</code></pre>
<p>
By paging through all the baselines, you can figure out which bad data affect which antennas (i.e., the same bad data is on all baselines to that antenna. These commands are applied in <code class="language-bash">flagdata</code> using <code class="language-bash">mode='list'</code>, <code class="language-bash">inpfile='**'</code> where <code class="language-bash">inpfile</code> is the name of your flagging file.
</p>
<p>
<ul>
<li>Take a look at step 9, and enter the correct inputs into <code class="language-bash">flagdata</code> for you to apply the text file you have created. This step will re-plot the source to ensure that the data is now flagged<pre class="line-numbers" data-start='262'><code class="language-python"> flagdata(vis='all_avg.ms',
mode='**',
inpfile='**')</code></pre></li>
<li>Execute step 9 and inspect the plotting window. You should get a plot similar to that shown below.</li>
</ul>
</p>
<div class="col-12"><span class="image fit"><img src="plots/0319+4130_flagged_avg_amp-time.png" alt="" /></span></div>
<p>
<ul>
<li><b>If you have time,</b> expand your investigations to the other sources and enter flagging commands into your flagging file. You can then re-run step 9 to apply the new flags.</li>
<li>If you are running out of time, take a look at <code class="language-bash">all_avg_all.flags.txt</code> as the <code class="language-bash">inpfile</code> command and see how all the previous flagging commands can be applied at once.</li>
<li>Execute step 10 to apply these flags to ensure that your data is ready for the next tutorial: Calibration</li>
</ul>
</p>
<p>
This step should also make some plots named <code class="language-bash"><source_name>_flagged_avg_amp-time.png</code>. Inspect these to check that all the bad data has been removed. A plot of amplitude versus time for the phase calibrator (1302+5748) is given below.
</p>
<div class="col-12"><span class="image fit"><img src="plots/post_flagging_1302+5748.png" alt="" /></span></div>
</section>
<hr/>
<section>
<h2 class="alt"><a name="calibration"></a>Part 2 - Calibration</h2>
<h3 class="alt"><a name="flux_scaling"></a>3. Setting the flux of the primary flux scale calibrator (step 11)</h3>
<p class="alt"><a href="#"><font size="-5">← back to top</font></a></p>
<p>
The first thing that we are going to calibrate is to set the flux density of the flux scale calibrators (which is in this case 3C286 and 3C84). This is important as the correlator outputs the visibilities with a arbitrary scaling so it is not in physical units i.e., Janskys. However, the visibilities amplitudes are correct relative to each other. This means that we can observe a standard source with a physically known flux density (similar to the zero points of a temperature scale) which then can be used to rescale the visibilities and put all fields onto this scale.
</p>
<p>
In this observation, we observed 1331+3030 (= 3C286) is an almost non-variable radio galaxy and its total flux density is very well known as a function of time and frequency (<a href="http://adsabs.harvard.edu/abs/2013ApJS..204...19P" target="_blank">Perley & Butler 2013</a>). However, this flux density standard was derived through VLA observations and the source is somewhat resolved by e-MERLIN.
</p>
<p>
To mitigate this, we use a model made from a ~12 hr observation centred on 5.5 GHz. This is scaled to the current observing frequency range (around 5 GHz) using <code class="language-bash">setjy</code> which derives the appropriate total flux density, spectral index and curvature using parameters from Perley & Butler (2013).
</p>
<p>
However, a few percent of the flux is resolved-out by e-MERLIN; the simplest way (at present) to account for this is later to scale the fluxes derived for the other calibration sources by 0.938 (based on calculations for e-MERLIN by Fenech et al.).
</p>
<p>
The model is a set of delta functions and <code class="language-bash">setjy</code> enters these in the Measurement Set so that their Fourier Transform can be used as a $uv$-plane model to compare later with the actual visibility amplitudes. If you make a mistake (e.g. set a model for the wrong source) use task <code class="language-bash">delmod</code> to remove it.
</p>
<p>
<ul>
<li>Take a look at step 11. Enter the correct parameters for the <code class="language-bash">setjy</code> task and flux scaling of 3C286.
<pre class="line-numbers" data-start="317"><code class="language-python"> setjy(vis='**',
field='**',
standard='Perley-Butler 2013',
model='**') # You should have downloaded and decompressed this model </code></pre></li>
<li>Have a look at the other commands in step 11 which are completed for you. By default, <code class="language-bash">setjy</code> scales the model for each channel, as seen if you plot the model amplitudes against $uv$ distance, coloured by spectral window. The first <code class="language-bash">plotms</code> command (shown below) will plot the model versus $uv$ distance. The model only appears for the baselines present for the short time interval for which 1331+3030 was observed here. <pre class="line-numbers" data-start="323"><code class="language-python"> plotms(vis='all_avg.ms', field='1331+3030', xaxis='uvwave',
yaxis='amp', coloraxis='spw',ydatacolumn='model',
correlation='LL,RR',symbolsize=4,
showgui=gui,overwrite=True,
plotfile='3C286_model.png')</code></pre></li>
</ul>
<div class="col-12"><span class="image fit"><img src="plots/3C286_model.png" alt="" /></span></div>
<p>
<ul>
<li>While this is not normal for e-MERLIN data, we are also setting the flux density of the bandpass calibrator (0319+4130) in this step. Normally, we would derive it's flux density from 3C286 and then derive the bandpass calibration twice (once without taking into account spectral indices and then taking it into account). This source is pretty well known so, to save extra calibration steps, we are going to set its flux density here and the spectral index, <code class="language-python">spix</code>, with the following code (done for you).<pre class="line-numbers" data-start="329"><code class="language-python"> setjy(vis='all_avg.ms',
field='0319+4130',
standard='manual',
fluxdensity=[30.1388, 0.0, 0.0, 0.0],
spix = 0.6, reffreq= '5.07GHz')</code></pre></li>
<li>This source is not resolved by e-MERLIN therefore the model that we are going to use is that of a point source. Hence, the following <code class="language-python">plotms</code> step should have a slope of 0.6 in amplitude versus frequency with a flux density of $30.1388\,\mathrm{Jy}$ at $5.07\,\mathrm{GHz}$<pre class="line-numbers" data-start="336"><code class="language-python"> plotms(vis='all_avg.ms', field='0319+4130', xaxis='freq',
yaxis='amp', coloraxis='spw',ydatacolumn='model',
correlation='LL,RR',symbolsize=4,
showgui=gui,overwrite=True,
plotfile='3C84_model.png')</code></pre></li>
</ul>
</p>
<div class="col-12"><span class="image fit"><img src="plots/3C84_model.png" alt="" /></span></div>
</section>
<hr/>
<section>
<h3 class="alt"><a name="bandpass"></a>4. Bandpass calibration</h3>
<h4 class="alt"><a name="decorrelation"></a>4A. Investigating decorrelation and solution intervals (step 12)</h4>
<p class="alt"><a href="#"><font size="-5">← back to top</font></a></p>
<p>
We have set the primary fluxes of the main flux / bandpass calibrators and now we can move onto calibrating for the instrumental/atmospheric effects. If you remember from your interferometry introduction, we have used a model of the locations of our telescopes which corrects for the geometric delay and ensures that our signal from the source interferes constructively between all antennae. However, the atmosphere, electronics/engineering of the antennae, and clock/timing errors changes the phase and prevents perfect constructive interference. We therefore have to correct for this.
</p>
<p>
To correct for this, our standard calibration assumes that your calibrator sources are point-like (i.e., flat in amplitude and zero in phase). This means that any deviation must be due to an error. Secondly, there is an assumption that the errors are antenna independent i.e., there must be a common factor/error on each baseline to the same antenna. This is sensible as each antenna will have different clocks (in the case of VLBI), different atmospheric line-of-sights, different LNAs etc etc. This also means that we should be able to derive a per antenna solution and prevents us from simply making our visibilities looking like a point source (if we solved per baseline).
</p>
<p>
We are going to calibrate the bandpass calibrator (0319+4130) first with the goal of generating a bandpass calibration table. The bandpass is the response of the antenna (both amplitude and phase) against frequency. This changes quickly with frequency, but is not expected to change against time (at least during the observation) as it is typically due to filters and receiver sensitivities within the antenna. This means we need to get solutions on each channel so we need a very bright source (0319 is 20-30 Jy so very bright) in order to get enough signal on each small bandwidth of the channel. However, before we can derive this correction, we need to first calibrate the other sources of error affecting this source first in order to isolate the bandpass correction.
</p>
<p>
To begin, we are going to correct for the phase errors, $\Delta\Theta$, which can primarily be parameterised by a constant plus higher order terms: $$\Delta\Theta = \Theta(t,\nu) + \frac{\partial \Theta}{\partial \nu} + \frac{\partial \Theta}{\partial t}$$
We call $\Theta(t,\nu)$ the <b>phase</b> error, the slope in frequency space, $\frac{\partial \Theta}{\partial \nu}$, the <b>delay</b> error and the slope against time, $\frac{\partial \Theta}{\partial t}$, the <b>rate</b> error. The rate error and higher order terms are only typically needed for long baseline or low frequencies where more turbulent atmospheric conditions and separate clocks cause errors. Most modern interferometers including e-MERLIN often only need to consider the delay and phase terms.
</p>
<p>
<ul>
<li>Take a look at step 12. It is advised for this step to copy the individual <code class="language-python">plotms</code> commands directly into CASA!</li>
<li>Copy the first command into CASA which will plot the phase versus time for a bright source i.e., our bandpass calibrator (so we can see the visibilities above the noise). Note that we are selecting just one scan of this source.
<pre class="line-numbers" data-start='352'><code class="language-python"> plotms(vis='all_avg.ms', field='0319+4130',
xaxis='frequency',yaxis='phase',
gridrows=5,gridcols=1,iteraxis='baseline',
xselfscale=True,xsharedaxis=True,
ydatacolumn='data',antenna=antref+'&*',
avgtime='2400',correlation='LL,RR',
coloraxis='corr', scan='92',
plotfile='',plotrange=[-1,-1,-180,180],
showgui=gui, # change to True to inspect by interval or antenna
overwrite=True)</code></pre></li>
</ul>
</p>
<p>
The y-axis spans $-180$ to $+180$ deg and some data has a slope of about a full 360 deg across the 512 MHz bandwidth. Looking at the baseline the Cambridge, think about the following questions.
</p>
<p>
<b>QUESTION 1:</b> What is the apparent delay error this corresponds to? What will be the effect on the amplitudes?<div class="popup" onclick="myFunction('Q1')"><b>Click for answer</b><span class="popuptext" id="Q1">A change of $2\pi$ radians in $x$ Hz results from a delay error of $1/x$ sec. So a change of 360 deg in 512 MHz is produced by a delay error of $1/(512\times 10^{6})\,\mathrm{sec} \sim 2\,\mathrm{ns}$. If data are vector-averaged across a phase slope, the amplitudes will be reduced, by a larger amount the greater the phase change in the averaging interval, so if the combined phase change on a baseline is worse in one polarisation than the other, that polarization will have a lower averaged apparent amplitude.</span></div>
</p>
<p>
<b>QUESTION 2:</b> Delay corrections are calculated relative to a reference antenna which has its phase assumed to be zero, so the corrections factorised per other antenna include any correction due to the refant. Roughly what do you expect the magnitude of the Cm corrections to be? Do you expect the two hands of polarisation to have the same sign?
<div class="popup" onclick="myFunction('Q2')"><b>Click for answer</b><span class="popuptext" id="Q2">
The slopes on Mk2 correspond to about 1 turn in 512 MHz, in the same sense for both LL and RR. The slopes on Cm correspond to less than half a turn in the same sense for one polarisation, and about a full turn in the opposite sense for the other. So the combination of errors, which will be solved entirely by correcting Cm, corresponds to magnitudes of $\lesssim 1\,\mathrm{ns}$ and $\sim 4\,\mathrm{ns}$ for the two polarisations, with opposite signs. The slopes, and hence the actual values, vary between spw, see the delay correction plot.</span></div>
</p>
<div class="col-12"><span class="image fit"><img src="plots/pre-delay_mk2_cm.png" alt="" /></span></div>
<p>
<ul>
<li>Have a look at the second and third <code class="language-python">plotms</code> commands which shows what happens if we average these data without correcting for the delay slope. The first command plots the amplitude of one channel and one correlation on the longest baseline.<pre class="line-numbers" data-start="365"><code class="language-python"> plotms(vis='all_avg.ms', field='0319+4130',
spw='0~3:55', ydatacolumn='data', # Just one channel per spw
yaxis='amp', xaxis='time', plotrange=[-1,-1,0,0.025], # Plot amplitude v. time with fixed y scale
avgtime='60',correlation='RR', antenna=antref+'&Cm', # 60-sec average, one baseline/pol
coloraxis='spw',plotfile='',
showgui=gui,
overwrite=True)</code></pre>The next command plots the same but now we are averaging the channels together, <pre class="line-numbers" data-start="375"><code class="language-python"> plotms(vis='all_avg.ms', field='0319+4130',
xaxis='time', yaxis='amp',ydatacolumn='data',
spw='0~3',avgchannel='64',
antenna=antref+'&Cm', avgtime='60',correlation='RR',
plotrange=[-1,-1,0,0.025],coloraxis='spw',
plotfile='',
showgui=gui, # change to T to inspect closely
overwrite=True)</code></pre></li>
<li><b>QUESTION 3:</b> Take a look at the two plots and compare. Which has the higher amplitudes - the channel-averaged data or the single channel? Why?<br/><div class="popup" onclick="myFunction('Q3')"><b>Click for answer</b><span class="popuptext" id="Q3">The single channel has the higher amplitudes. The source 3C84 is very bright so there is plenty of signal in each single channel. But averaging over all channels when there is a large phase slope causes the amplitudes to decorrelate so the channel-averaged data has depressed amplitudes. A similar effect is seen if you have data with phase errors as a function of time and average over too long a time period.</span></div></li>
</ul>
</p>
<p>
We have plotted the same amplitude versus time plots below in python for the four scans of the bandpass calibrator. As you can see, the effects of decorrelation of the amplitudes is more extreme when the delay error is larger (e.g., in scan 4)
</p>
<div class="col-12"><span class="image fit"><img src="plots/delay_average_amp.png" alt="" /></span></div>
<p>
Finally, in preparation for the calibration of the bandpass source, we want to derive the time-dependent phase solution interval. A solution interval is just the amount of time you average together to try and derive a correction. The selection of a solution interval can be tricky but it is all about timescales. On an integration to integration basis (few seconds), your visibilities will be dominated by thermal noise while on longer timescales, the general trends will be dominated by the instrumental / line-of-sight effects that we want to calibrate out. We therefore want to average enough data to track these effects and derive corrections but not over-average that we do not pick up on the detailed trends. We shall show this next
</p>
<p>
<ul>
<li>Take a look at the final <code class="language-python">plotms</code> command in step 12,<pre class="line-numbers" data-start="386"><code class="language-python"> plotms(vis='all_avg.ms',
field='0319+4130',
spw='0~3:55', correlation='RR',
antenna=antref+'&Cm',scan='1',
xaxis='time', yaxis='phase',ydatacolumn='data',
coloraxis='spw',plotfile='',
showgui=gui, overwrite=True)</code></pre></li>
<li>We are just plotting the phase versus time of a single scan, single correlation and single channel of the bandpass calibrator. You should see a plot like that below (colorised by spw)</li>
</ul>
</p>
<div class="col-12"><span class="image fit"><img src="plots/0319+4130_1ch_phase-time.png" alt="" /></span></div>
<p>
Some things to note here:
<ol>
<li>On the 20-30s timescales we have trackable phase fluctutations caused by the atmosphere which causes the phases to drift.</li>
<li>These fluctuations are similar across each spw indicating that the atmosphere affects our frequency range equally (meaning that we could fix the delay and combine spw for higher S/N solutions</li>
<li>The delays are causing the absolute phase offsets between the spectral windows, but, crucially, these seem to not change with time (otherwise the gap between the spectral windows would change with time). Because of this, it is suitable to select a long solution interval (e.g., 600s in this case) for the delay calibration to maximise the S/N of our solutions.</li>
</ol>
<p>
While we have decided on a solution interval for delays, we need to do the same for the phase term. In the plot below, we have overplotted the phase with different averaging intervals of just one spectral window. You can see that we start being unable to trace the general trends causes by the atmosphere when the averaging is too severe. Howevever, we can also see that there's negligible thermal noise on a intergration to integration basis, so we could actually use a very short solution interval such as the 30s shown in the plot. This would track the atmospheric fluctuations well.
</p>
<div class="col-12"><span class="image fit"><img src="plots/phase_solution_interval.png" alt="" /></span></div>
</section>
<hr/>
<section>
<h4 class="alt">4B. <a name="pre_bandpass">Pre-bandpass calibration (step 13)</a></h4>
<p>
With our solution intervals in hand, we are going to now conduct the gain calibration for us to isolate the bandpass corrections. This will allow us to calibrate out the atmosphere, clock errors etc etc. which changes our amplitudes and phases over time.
</p>
<p>
<ul>
<li>Again, for this step, it is <b>advisable</b> to copy and paste your tasks directly into CASA to do each calibration step bit by bit. If you are having syntax errors, remember that the variables <code class="language-python">gui = True</code> and <code class="language-python">antref = 'Mk2'</code> should have been set in CASA.</li>
<li>Have a look at step 13. We shall begin by calibrating the delays by using the CASA task <code class="language-python">gaincal</code>. We saw from the figure in the previous section that the absolute offset between the phases of different spws remained almost the same hence the delay calibration did not change too much (hence can use a long solution interval). Fill in the missing code to generate the delay solutions for the bandpass calibrator. Remember to use <code class="language-python">help(gaincal)</code>!<pre class="line-numbers" data-start="404"><code class="language-python"> gaincal(vis='all_avg.ms',
gaintype='**',
field='**',
caltable='bpcal.K',
spw='0~3',solint='**',
refant=antref,
refantmode='strict',
minblperant=2,
minsnr=2)</code></pre></li>
</ul>
</p>
<p>
The <code class="language-python">gaincal</code> task will generate a calibration table (name <code class="language-bash">bpcal.K</code> in this case) which will include the delay solutions that, when applied to the visiblities, will remove that phase slope versus frequency we saw earlier. However, we <b>always</b> want to plot the solutions to ensure that they are decent. Good solutions are normally smoothly varying without jumps.
</p>
<p>
<ul>
<li>Take a look at the next <code class="language-bash">plotms</code> command and enter it into CASA<pre class="line-numbers" data-start="414"><code class="language-python"> plotms(vis='bpcal.K',
xaxis='time',
yaxis='delay',
gridrows=3,gridcols=2,
iteraxis='antenna', coloraxis='spw',
xselfscale=True,
yselfscale=True,
plotfile='',
showgui=gui, overwrite=True)</code></pre></li>
</ul>
</p>
<p>
This should produce a plot like that shown below.
</p>
<div class="col-12"><span class="image fit"><img src="plots/bpcal.K.png" alt="" /></span></div>
<p>
As you can see, the solutions are smooth and generally follow each other. Note the large delay correction for the Cm telescope in the last scan. We saw this in Section 4A where we plotted the amplitude decorrelation via averaging the channels. The huge delay (shown by the phase wrapping quickly), corresponds to the 35ns solutions for this scan on this telescope.
</p>
<p>
With the delay calibration for 0319+4130 complete, the offsets between the spectral windows in the frequency space would be removed and now we can concentrate on the time-dependent phase calibration which will remove the effect of the atmosphere along our line-of-sight to the source.
</p>
<p>
<ul>
<li>Take a look at the second <code class="language-bash">gaincal</code> command and enter in the parameters to perform a phase only calibration of the bandpass calibrator. Note that we include the delay calibration table in the <code class="language-bash">gaintable</code> parameter. CASA conducts calibration on the fly hence the delay table will be applied to the visibilities before the phase calibration is solved for.<pre class="line-numbers" data-start="427"><code class="language-python"> gaincal(vis='all_avg.ms',
calmode='**',
field='**',
caltable='bpcal_precal.p1',
solint='**',
refant=antref,
refantmode='strict',
minblperant=2,
gaintable=['bpcal.K'],
minsnr=2)</pre></code></li>
<li>Again we want to plot the solutions to ensure that they are not just noise. <pre class="line-numbers" data-start="427"><code class="language-python"> plotms(vis='bpcal_precal.p1',
xaxis='time',
yaxis='phase',
gridrows=3,gridcols=2,
iteraxis='antenna', coloraxis='spw',
xselfscale=True,
yselfscale=True,
plotfile='',
showgui=gui, overwrite=True)</code></pre></li>
</ul>
</p>
<div class="col-12"><span class="image fit"><img src="plots/bpcal_precal.p1.png" alt="" /></span></div>
<p>
These solutions look great (I used a 30s interval`) and seem to follow the phases (and fluctuation timescales) we had seen before. Finally, we have one calibration step missing. We've done delays and phases and it's just amplitudes that are left!
</p>
<p>
<ul>
<li>Look at the final <code class="language-python">gaincal</code> command and enter the correct parameters for amplitude only calibration of the calibrator. You will need to enter the <code class="language-bash">gaintable</code> parameter in this time. Note that we did not talk about the amplitude solution interval but typically these change slower than phases but need a higher S/N for good solutions (see 'closure amplitudes'). In this case, we will use 120s.<pre class="line-numbers" data-start="451"><code class="language-python"> gaincal(vis='all_avg.ms',
calmode='**',
field='**',
caltable='bpcal_precal.a1',
solint='120s',
solnorm=False,
refant=antref,
refantmode='strict',
minblperant=2,
gaintable=['**','**'],
minsnr=2)</code></pre></li>
<li>Finally, we want to plot these solutions to ensure they look ok and not noisy.<pre class="line-numbers" data-start="465"><code class="language-python"> plotms(vis='bpcal_precal.a1',
xaxis='time',
yaxis='amp',
gridrows=3,gridcols=2,
iteraxis='antenna', coloraxis='spw',
xselfscale=True,
yselfscale=True,
plotfile='',
showgui=gui, overwrite=True)</code></pre></li>
</ul>
</p>
<div class="col-12"><span class="image fit"><img src="plots/bpcal_precal.a1.png" alt="" /></span></div>
<p>
These solutions look good. Note that the amplitudes are very low and that is to rescale the visibilities to the flux density that we had set right at the start of this section. We will do the same later for the other sources which we don't know the true flux densities of.
</p>
</section>
<hr/>
<section>
<h4 class="alt"><a name="init_bandpass">4C. Derive bandpass solutions (step 14)</a></h4>
<p>
With all other corruptions now having solutions through the three calibration tables we have derived. The final step is to derive the bandpass calibration table using a special task called <code class="language-python">bandpass</code>.
</p>
<p>
<ul>
<li>Take a look at step 14, enter the correct parameters for the bandpass calibration and then enter this into the CASA prompt.<pre class="line-numbers" data-start="482"><code class="language-python"> bandpass(vis='all_avg.ms',
caltable='bpcal.B1',
field='**',
fillgaps=16,
solint='inf',combine='scan',
solnorm='**',
refant=antref,
bandtype='B',
minblperant=2,
gaintable=['**','**','**'],
minsnr=3)</code></pre></li>
</ul>
</p>
<p>
You may find a few complaints about failed solutions but these should be for the end channels that are flagged anyways.
<ul>
<li>Inspect the solutions using the <code class="language-python">plotms</code> commands at the end of step 14<pre class="line-numbers" data-start="495"><code class="language-python"> plotms(vis='bpcal.B1',
xaxis='freq',
yaxis='amp',
gridrows=3,gridcols=2,
iteraxis='antenna', coloraxis='corr',
xselfscale=True,
yselfscale=True,
plotfile='',
showgui=gui, overwrite=True)</code></pre> for amplitude and, <pre class="line-numbers" data-start="506"><code class="language-python"> plotms(vis='bpcal.B1',
xaxis='freq',
yaxis='phase',
gridrows=3,gridcols=2,
iteraxis='antenna', coloraxis='corr',
xselfscale=True,
yselfscale=True,
plotfile='',
showgui=gui, overwrite=True)</code></pre> for phase.</li>
</ul>
</p>
<p>
These should produce plots like those below which will correct for the receiver response on each antenna (note the close comparison to the plots when we flagged the edge channels). This error is antenna-based and direction-independent and so we can use it for all other sources as it should not change with time. Note that the spectral index of this source is taken into account here as it was specified in step 3A (normally you'd do the bandpass calibration twice with the second time using an estimate of the spectral index).
</p>
<div class="col-12"><span class="image fit"><img src="plots/bpcal.B1_amp.png" alt="" /></span></div>
<div class="col-12"><span class="image fit"><img src="plots/bpcal.B1_phase.png" alt="" /></span></div>
</section>
<hr/>
<section>
<h3 class="alt"><a id="delay"></a>Time-dependent delay calibration (15)</h3>
<p class="alt"><a href="#"><font size="-5">← back to top</font></a></p>
<p>
We've calibrated the bandpass calibrator in order to obtain the bandpass corrections. We now need to derive the time-dependent phase and amplitude corrections for all of the calibrators at once. This will allow us to perform phase referencing (remember this from the lecture) so we can observe the target source. Same as before, we will start with the delays and then will move onto the phase and amplitude corrections.
</p>
<p>
<ul>
<li>Take a look at step 15 and enter the parameters for both the gaincal command and the plotting command so that we can do delay calibration for all calibrator sources. The solution interval we used for the bandpass calibrator should be applicable here too. Remember that the bandpass calibration table should also be applied now!<pre class="line-numbers" data-start="524"><code class="language-python"> gaincal(vis='all_avg.ms',
gaintype='K',
field='**',
caltable='all_avg.K',
spw='0~3',solint='**',
refant=antref,
refantmode='strict',
gaintable=['**'],
minblperant=2,
minsnr=2)
plotms(vis='all_avg.K',
xaxis='**',
yaxis='**',
gridrows=3,gridcols=2,
iteraxis='antenna', coloraxis='spw',
xselfscale=True,
yselfscale=True,
plotfile='',
showgui=gui, overwrite=True)</code></pre></li>
<li>Once you are happy execute step 15 and look at the output of <code class="language-python">plotms</code></li>
</ul>
</p>
<p>
You should end up with a plot that looks like this below. The solutions are smoothly varying with only a few jumps which shows that our solutions are ok.
</p>
<div class="col-12"><span class="image fit"><img src="plots/all_avg.K.png" alt="" /></span></div>
</section>
<hr/>
<section>
<h4 class="alt"><a id="inspect_bandpass"></a>4A. Inspect effects of calibration (16)</h4>
<p>
<b>Important:</b> if you are short on time, skip this step and just read the information and plots shown here
</p>
<p>
In this step, we shall apply the delay and bandpass solutions to the calibration sources as a test, to check that they have the desired effect of producing a flat bandpass for the phase-reference source (not on the fly this time!). There is no need to apply the time-dependent calibration as we have not derived this for the phase-reference source yet anyway.
</p>
<p>
To apply these we shall use a task called <code class="language-bash">applycal</code>. Some notes about <code class="language-bash">applycal</code>:
<ol>
<li>The first time <code class="language-python">applycal</code> is run, it creates a CORRECTED data column in the MS, which is quite slow.</li>
<li>Each time <code class="language-python">applycal</code> is run, it replaces the entries in the CORRECTED column, i.e. the corrected column is not cumulative, but you can give <code class="language-python">applycal</code> a cumulative list of gain tables.</li>
<li><code class="language-python">applymode</code> can be used to flag data with failed solutions but we do not want that here as this is just a test.</li>
<li>The <code class="language-python">interp</code> parameter can take two values for each gaintable to be applied, the first determining interpolation in time and the second in frequency. The times of data to be corrected may coincide with, or be bracketed by the times of calibration solutions, in which case <code class="language-python">linear</code> interpolation can be used, the default, normally used if applying solutions derived from a source to itself, or from phase ref to target in alternating scans. However, if solutions are to be extrapolated in time from one source to another which are not interleaved, then nearest is used. A similar argument covers frequency. There are additional modes not used here (see <code class="language-python">help(applycal)</code>).</li>
</ol>
The delay correction ('K') table contains entries for all sources except the target. so we can use linear interpolation for this one. The bandpass table only has one entry so we cannot interpolate hence we'd want to select the nearest solutions.
</p>
<p>
<ul>
<li>Fill in the <code class="language-python">applycal</code> step in step 16 and look at the following <code class="language-python">plotms</code> command which will plot the 'before and after' phase reference phase and amplitude against frequency.<pre class="line-numbers" data-start="555"><code class="language-python"> applycal(vis='all_avg.ms',
field='1302+5748,0319+4130,1331+3030',
calwt=False,
applymode='calonly',
gaintable=['**','**'],
interp=['**','**,**'])</code></pre></li>
<li>Execute step 16 and take a look at the plots saved to the current directory.</li>
</ul>
<p>
The plot below shows the uncorrected and corrected amplitudes and phases for the Mk2-Cm baseline of the phase reference source (1302+5748) that is coloured by scan. It has been replotted in python for clarity!
</p>
<div class="col-12"><span class="image fit"><img src="plots/applycal_effects.png" alt="" /></span></div>
<p>
In the figure, each scan is coloured separately and shows a random offset in phase and some discrepancies in amplitude compared to other scans, but for each individual time interval it is quite flat. So we can average in frequency across each spw and plot the data to show remaining, time-dependent errors.
</p>
</section>
<hr/>
<section>
<h3 class="alt"><a name="gaincal"></a>6. Gain calibration (step 17)</h3>
<h4 class="alt"><a name="time_dep_phase"></a>6A. Time dependent phase calibration</h4>
<p class="alt"><a href="#"><font size="-5">← back to top</font></a></p>
<p>
As with the bandpass calibrator, we now move onto correcting for the phase fluctuations caused by tubulence in the atmosphere. This will allow us to phase reference and see our target source as the phase calibrator is approximately along the same line-of-sight.
</p>
<p>
<ul>
<li>Take a look at step 17 and insert the parameters that will allow us to conduct time-dependent phase calibration of all calibrator sources<pre class="line-numbers" data-start="623"><code class="language-python"> gaincal(vis='all_avg.ms',
calmode='p',
caltable='calsources.p1',
field='**,**,**',
solint='**',
refant=antref,
refantmode='strict',
minblperant=2,minsnr=1,
gaintable=['all_avg.K','bpcal.B1'],
interp=['linear','nearest,nearest'])</code></pre></li>
<li>Copy this into the CASA propmt and run to generate the calibration table and have a look at the logger</li>
</ul>
</p>
<p>
The logger should show something like (if you used 16s solution intervals like me):<pre><code class="language-none">2022-08-12 10:08:20 INFO Writing solutions to table: calsources.p1
2022-08-12 10:08:21 INFO calibrater::solve Finished solving.
2022-08-12 10:08:21 INFO gaincal::::casa Calibration solve statistics per spw: (expected/attempted/succeeded):
2022-08-12 10:08:22 INFO gaincal::::casa Spw 0: 1988/1572/1571
2022-08-12 10:08:22 INFO gaincal::::casa Spw 1: 1988/1572/1571
2022-08-12 10:08:22 INFO gaincal::::casa Spw 2: 1988/1572/1571
2022-08-12 10:08:22 INFO gaincal::::casa Spw 3: 1988/1572/1571
</code></pre>
</p>
<p>
This means that although there are 1988 16s intervals one the calibration sources data, about a fifth (1988−1572=416) of them have been partly or totally flagged. But almost all (99.93%) of the intervals which do have enough unflagged ('attempted') data gave good solutions ('succeeded'). There are no further messages about completely flagged data but in the terminal you see a few messages when there are fewer than the requested 2 baselines per antenna unflagged:<pre><code class="language-none">Found no unflagged data at: (time=2019/08/02/21:00:04.5 field=0 spw=3 chan=0)
Insufficient unflagged antennas to proceed with this solve.
(time=2019/08/02/21:17:19.2 field=0 spw=3 chan=0)</code></pre>
</p>
<p>
Here, there are very few predicable failures. If there were many or unexpected failures, investigate - perhaps there are unflagged bad data, or pre-calibration has not been applied, or an inappropriate solution interval was used.
</p>
<p>
<ul>
<li>Plot the solutions and you should get something which tracks the phases well (as seen below)</li>
</ul>
</p>
<div class="col-12"><span class="image fit"><img src="plots/calsources.p1.png" alt="" /></span></div>
</section>
<hr/>
<section>
<h4 class="alt"><a name="time_dep_phase"></a>6B. Time dependent amplitude calibration</h4>
<p class="alt"><a href="#"><font size="-5">← back to top</font></a></p>
<p>
Next we are going to derive time-dependent amplitude solutions for all the calibration sources, applying the delay, bandpass tables and the time-dependent phase solution table.
</p>
<p>
<ul>
<li>Take a look at the second <code class="language-python">gaincal</code> entry in step 17 and Iisert a suitable solution interval, the tables to apply and interpolation modes. Note that if you set a solint longer than a single scan on a source, the data will still only be averaged within each scan unless an additional 'combine' parameter is set (not needed here).<pre class="line-numbers"
data-start="648"><code class="language-python"> gaincal(vis='all_avg.ms',
calmode='a',
caltable='calsources.a1',
field='**,**,**',
solint='**',
solnorm=False,
refant=antref,
minblperant=2,minsnr=2,
gaintable=['**','**','**'],
interp=['linear','nearest,nearest','**'])</code></pre></li>
<li>The logger shows a similar proportion of expected:attempted solution intervals as for phase, but those with failed phase solutions are not attempted so all the attempted solutions work. Plot these solutions to ensure they look ok</li>
</ul>
</p>
<div class="col-12"><span class="image fit"><img src="plots/calsources.a1.png" alt="" /></span></div>
<p>
The plot above shows the time dependent amplitude solutions on just the L polarisation and coloured by spw. The solutions for each separate source look similar but there are big differences between sources, since they have all (apart from 1331+3030 and 0319+4130) been compared with a 1 Jy model but they have different flux densities.
</p>
</section>
<hr/>