Skip to content
Gijs Molenaar edited this page Feb 13, 2014 · 5 revisions

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

1. 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
  1. We only work with the first subband (IVC0). So only keep it and delete the rest.
  2. 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
[[duchamp_exmap.png]] This shows the extracted sources in red. 

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

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. 

[[with_no_wproj.gif]] 

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. 
[[mean_and_res.gif]] 

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

[[3c196_CygA.png]] 

[[3c196_TauA.png]] 

[[3c196_VirA.png]] 


### Comparison with NVSS

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

[[3c196_NVSS.png]] 

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: 

[[3c196_tiling.png]] 


### 12hrx6spacings

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

[[3c196_all.gif]] 

[[3c196_allbw.gif]] 
Clone this wiki locally