• MyCalib3C193

3C196

A Priori Information

1. 3C196: a point source at 08:13:36.05, +48:13:02.26
2. Flux at 120 MHz ~ 70 Jy
3. The whole purpose is to test the astrometric accuracy of the calibration.

First Steps

1. Get the data from the archive, using getms: spacings 36m and 96m.
   • "./getms.csh" and follow instructions
2. Split subbands, combine in time using WSRT glish scripts.

   for f in 10705866_S0_T[0-9].MS; do
    glish -l MSIVCsplit.g $f;
    rm -rf $f;
   done
   for f in 10705866_S0_T[1-9][0-9].MS; do
    glish -l MSIVCsplit.g $f;
    rm -rf $f;
   done
   let "i=0";
   while [ $i -le 7 ] ; do
    str="10705866_S0_T?_IVC"$i".MS";
    str1="10705866_S0_T1?_IVC"$i".MS";
    str2="10705866_S0_T2?_IVC"$i".MS";
    str3="10705866_S0_T3?_IVC"$i".MS";
    str4="10705866_S0_T4?_IVC"$i".MS";
    str5="10705866_S0_T5?_IVC"$i".MS";
    outfile="out_IVC"$i".ms"
    echo $str;
    echo $str1;
    echo $str2;
    echo $str3;
    echo $str4;
    echo $str5;
    glish -l MSTimeCombine.g  $str $str1 $str2 $str3 $str4 $str5 $outfile;
    rm -rf $str;
    rm -rf $str1;
    rm -rf $str2;
    rm -rf $str3;
    rm -rf $str4;
    rm -rf $str5;
    let "i = $i +1";
   # glish -l applyTaper.g $outfile hanning
   done

3. Flagging: routine flagging with clipping, median flagging and flagging bad antenna no 6.

   include 'autoflag.g'
   infile:=0;
   # clipping threshold
   limm:=1e5;
   # uv clip limit
   limuv:=1;
   ### parse args
   for( a in argv )
   {
     print 'arg: ',a;
     if( a =~ s/ms=// )
       infile:= a;
     else if( a =~ s/minuv=// )
       limuv:=as_integer(a);
     else if( a =~ s/minclip=// )
       limm:=as_float(a);
   }
   print spaste("Preprocessing:::",infile);
   af:=autoflag(infile)
   af.setdata()
   cliprec:=[=]
   cliprec[1] := [expr='ABS XX',min=1e-6,max=limm]
   cliprec[2] := [expr='ABS YY',min=1e-6,max=limm]
   cliprec[3] := [expr='ABS XY',min=1e-6,max=limm]
   cliprec[4] := [expr='ABS YX',min=1e-6,max=limm]
   cliprec[5] := [expr='UVD',min=limuv]
   af.setselect(clip=cliprec)
   af.run(reset=T)
   af.done()
   ## median flag
   af:=autoflag(infile)
   af.setdata()
   af.settimemed(thr=6,hw=5, expr='ABS XX')
   af.settimemed(thr=6,hw=5, expr='ABS XY')
   af.settimemed(thr=6,hw=5, expr='ABS YX')
   af.settimemed(thr=6,hw=5, expr='ABS YY')
   ### frequency median
   #af.setfreqmed(thr=6,hw=5, expr='ABS XX')
   #af.setfreqmed(thr=6,hw=5, expr='ABS YY')
   af.run()
   af.done()
   ## special flags, depenging on the measurement set
   af:=autoflag(infile)
   af.setdata()
   #af.setselect(ant=14)
   af.setselect(ant=6)
   af.run()
   af.done()
   ###################### compress
   include 'ms.g'
   include 'os.g'
   ## append new name to middle
   newms:=infile;
   newms=~s/.MS//g;
   newms=~s/.ms//g;
   newms:=sprintf("%s_M.MS",newms);
   ## check if file already exists
   if (dos.fileexists(newms)) {
     print sprintf("Error: File  %s already present, quitting\n",newms);
     exit
   } else {
   mm:=ms(infile,readonly=F)
   #mm.split(newms,nchan=[4],start=[32],step=[48]);
   mm.split(newms,nchan=[8],start=[200],step=[6]);
   mm.done()
   }
   #### re-run median flagger on new MS
   af:=autoflag(newms)
   af.setdata()
   af.settimemed(thr=6,hw=5, expr='ABS XX')
   af.settimemed(thr=6,hw=5, expr='ABS XY')
   af.settimemed(thr=6,hw=5, expr='ABS YX')
   af.settimemed(thr=6,hw=5, expr='ABS YY')
   af.run()
   af.done()
   ## open with imager
   include 'imager.g'
   ii:=imager(newms);
   ii.done()
   print spaste("Done Preprocessing:::",infile);
   exit

4. We only work with the first subband (IVC0). So only keep it and delete the rest.
5. The original MS has 512 channels. Reduce this to something smaller by using ms.split: We end up with 10 smaller MSes with 8 channels each. Each channel in these MS are average of 6 original channels. (see above script)

Calibration

1. We start off with only 3C193 in our sky model. Calibration is nothing fancy here, just G Jones (diagonal).
2. We calibrate the data, make images and use Duchamp for source extraction. Here's the parameter file for Duchamp:

   imageFile /home/sarod/center_mean.fits
   logFile   logfile.txt
   outFile   results.txt
   spectraFile spectra.ps
   minPix    10
   flagATrous  0
   snrRecon  10.
   snrCut    5.
   threshold 0.300
   minChannels 3
   flagBaseline    0
   flagKarma 1
   karmaFile duchamp.ann
   flagnegative 0

This shows the extracted sources in red.

1. We use the source model built by Duchamp as the updated sky model for the calibration.
2. I repeated the above steps 2,3 about 4 times to get a good sky model.

The above figure shows the area centered around 3c196, with w-projection (left) and without w-projection (right). They are the same right? wrong!!! It turns out, that enabling w-projection gives an error in source locations, especially away from the centre. The figure below is an animation illustrating the difference of using and not using w-projection.

Here is the WRONG way to make an image:

....
....
myimager.setoptions(ftmachine='wproject', wprojplanes=256, padding=1.2 , cache=500000000)
myimager.clean(algorithm="wfclark" , model=my_model , image=my_img , residual=my_res, threshold=[value=0.0, unit="Jy" ], niter=800, interactive=F, async=F, displayprogress=F)
....
...

Here is the way that works for WSRT:

....
....
myimager.setoptions(cache=500000000)
myimager.clean(algorithm="clark" , model=my_model , image=my_img , residual=my_res, threshold=[value=0.0, unit="Jy" ], niter=800, interactive=F, async=F, displayprogress=F)
....
....

1. I have reduced band 0 for spacing 36 and 96 m. The image below shows the mean image and the mean residual image for both spacings. Note, only some sources have been subtracted well. I dont assume non zero spectral indices to any of the sources and no fancy peeling/uvsubtraction has been done here.

1. Anyway, our goal is to image the A-Team. so here they are, imaged from the residual.

Comparison with NVSS

The image below is the NVSS source positions overlaid on the original image. To first approximation, they match very well!

But if you look really closely youll see a small rotation (counterclockwise) of the observed position due to precession.

J2000 or JAPP

The MS has UVW coordinates in JAPP. But it is assumed to be J2000. That is the reason for this rotation. The fix is simple:

 j2convert msin=input.MS type=J2000 force=True

Tiling in Frequency

Here's what happen if you try to solve for too many channels at one go:

12hrx6spacings

Here's the result (comparison) of traditional clean and uvsubtraction+restoration: