NB: If you want this walkthrough with images, visit Doug's webpage here.

 

A New Guide to DIFMAP

 

DIFMAP is a powerful and expedient program for imaging multi-element arrays. Originally intended for use with continuum VLBI data, it serves just as well when being used on VLA data sets. DIFMAP uses the increasing amounts of RAM in modern computers to minimize time spent reloading data and regenerating models. This tutorial will walk you through the basic steps of using DIFMAP to clean and self-calibrate your data.

 

1. Starting DIFMAP

 

All of the Brandeis workstations should have DIFMAP installed and ready to use on them. If you need to install DIFMAP, refer to the wiki page at http://brandeisastro.pbwiki.com/Installing+Difmap.

 

With DIFMAP installed, open a terminal and navigate to the directory in which your data resides. Because most of us use AIPS to perform the first round of data flagging/calibration, many people prefer to keep files in the /home/aips/FITS/ directory for convenience as this is where AIPS reads in and produces files.

 

Once you have navigated to the proper directory type “difmap” at the command prompt. This will initiate DIFMAP and we can begin loading data.

 

It should be noted that for this tutorial we will be using C Band C Array data for source J1036+1326. This is from Project AT0205, Seg A, source TEX1033+137, taken on 7/20/97 with 370 seconds TOS.

 

 

 

2. Loading Data

 

Now that we have loaded DIFMAP, you will be confronted with a new prompt:

 

0>

 

At this prompt you will want to give your first command to load data:

 

0> observe YOUR_DATA

 

Note that you can simply type obs followed by Tab to fill in the rest. This technique is helpful in saving some time throughout the DIFMAP experience. Learning when, and where, you can use this feature is worth some investigation on the part of the user. At this point you should see something similar to the following:

 

There are 2 IFs, and a total of 2 channels:

 

IF Channel Frequency Freq offset Number of Overall IF

origin at origin per channel channels bandwidth

------------------------------------------------------------- (Hz)

01 1 8.4149e+09 2.5e+07 1 2.5e+07

02 2 8.4649e+09 2.5e+07 1 2.5e+07

 

Polarization(s): RR LL RL LR

Read 123 lines of history.

Reading 21600 visibilities.

 

This will give you all the starting information about your source. Check to make sure all IFs, channels, frequencies, etc are correct. If you see an error such as “multi-source file” you may need to go back to AIPS and be certain to export the correct data. If at any time you want more information about your source you can type:

 

0> header

 

You will also notice that in the startup directory a new file has been created, difmap.log. This is a log file which will record all of the commands that you enter into DIFMAP. It will also include selected output lines though they will be preceded by a ! in order to comment them out. This will allow you to later use the log file as a script to run through your data. You can use this to recover up to the point you were at by typing @difmap.log to execute the script. It should be noted that the log file does not speak with PGPLOT, so any interactive flagging will not take place. For this reason this is not the best method for recovering your work after you exit DIFMAP. More on this later...

 

The next step is to choose the Stokes parameter you are interested in:

 

0> select i

 

Where i refers to the Stokes I parameter. You can select any Stokes parameter (e.g. select q or u) to work with.

3. Examining and Editing Your Data

 

In the beginning, a useful plot to examine is amplitude versus uv radius, otherwise known as a radial plot.

 

0> radplot

 

Because this is the first PGPLOT window you have activated, you will be required to specify to which device you want it to output data to. The best choice for this will be /xw. You can also choose this at any time by typing:

 

0> device /xw

 

Now lets look at your new plot. The following are useful keyboard commands available to you to manipulate the plot. If you desire more, press “h” while on the plot to display all options.

 

X - Quit radplot, you can also right click at any time.

L - Re-display whole plot, it is good to refresh this to add changes that may need to be updated.

n - Highlight next telescope

p - Highlight previous telescope

T - Specify highlighted telescope from keyboard

s - Show the baseline and time of the nearest point to the cursor

A - (Left-mouse-button) Flag the point closest to the cursor

C - Initiate selection of an area to flag.

Z - Select a new amplitude or phase display range.

U - Select a new UV-radius display range.

 

Some Helpful Hints:

 

1. Some people like to go through AIPS without doing any calibration to use DIFMAP to track down problems with the antennas. If this is your fancy, this is a good place to start looking.

2. Check often! Your radial plot will change as you improve your phase and amplitude calibrations throughout the DIFMAP process. Check often to keep your eye on the changes!

3. Antennas with severe problems will be easily seen and dealt with in DIFMAP. If you have an antenna with bad phases or amplitudes, you can identify it in your radial plot and sometimes fix it in the same plot. Use the “c” command to pick a window to delete if you can get the entire telescope in one pass. If not, continue on for other options for interactive flagging.

4. IT IS MOST IMPORTANT to remember that editing in the radial plot will be easy and tempting, but that you will be making BASELINE based edits, whereas systematic errors are more likely ANTENNA based. Be careful to not flag good data!

 

With our radial plot inspected, the next location for valuable source information is the UV plot. You can access this with:

 

0> uvplot

 

All of the commands for this plot are similar to the radplot commands.

 

Another useful plot is a time interval plot of the data:

 

0> tplot

 

A dot will appear for each time interval, for each telescope for which you have data. Data will be colored differently depending on what has been done to that data:

 

Green = No edits.

Yellow = Flagged.

Blue = Flagged by selfcal or in corplot.

Red = All antenna data flagged.

 

Lastly, we can look at a directional cut of the amplitude/phase along any radial line in the UV plane using the command:

 

0> projplot 20

 

Here the 20 is the angle (measured east from north) along which DIFMAP will display the projected amplitude and phase. The options are again similar to radplot except now you can use:

 

? to select projection angle

< and > to cycle through angles

 

Lastly there is the vplot command:

 

0> vplot 2

 

The number in front designates how many baselines will be displayed per page. You can press “p” to cycle back 2 baselines or “n” to cycle forwards. You can also use “spacebar” to switch from station to baseline editing. This is the best place to do antenna/baseline specific editing. Use the “t” function to quickly locate the antenna you are interested in.

 

 

4. The First Save

 

At this point you will have made all of the initial editing and visual inspections you need to make of your data. This is an excellent opportunity for you to make your first save.

 

0> save YOUR_DATA.1

 

It is helpful to utilize a method of saving that will let you move backwards and forwards in your work, especially when dealing with complex sources that may require multiple efforts to perfect. This should output several files, including .par, .mod, .fits, .uvf, and sometimes .win if you have any active windows. When you are ready to work with your data again, simply type:

 

0>@YOUR_DATA.1.par

 

This will reset you back to where you were at the time of your save. This is very helpful should you make a small but damaging mistake to your data and not want to start from square one.

 

 

5. Last Tidbits of Setup

 

At this point you are almost ready to work with your data. The last items you need to take care of before imaging relate to how you wish to construct your image.

 

0> uvweight #,#

 

This command will set the uv weighting scheme for your data. What you choose here will impact the image you will create as it sets the weighting the data receives when we do the FFT of the uv data into an image. The two most common weighting schemes are

 

0> uvweight 2,0 (or 2,-1) – uniform weighting

 

Uniform weighting is a good choice when trying to see sharp details of a source, though this is not the only reason one might utilize it. With uniform weighting, the data are weighted inversely to the number of visibilities occurring within the bins (unifoirm selects bins with 2 uv-grid pixel sizes.)

 

0> uvweight 0,-2 – natural weighting

 

Natural weighting is when we weigh all data equally. This provides a higher sensitivity than uniform weighting, but at the cost of the resolution of uniform weighting. You will now note that the first number is the bin size and the second number is the weighting of the visibilities by the errors. -1 or -2 is a good choice.

 

With our weighting chosen we now only need to specify the size and scale of our maps. The size of the map is specified as:

 

0> mapsize 1024,20

 

The first number dictates the number of pixels on each axis of the image. This number MUST be a power of 2 for the FFT to function, thus good scales are 1024, 2048, and 4096. At sizes of 8192 and even 16384, most computers will slow down significantly or even seize up. 16384 pixels is not recommended unless you have over 4 GB or RAM, sufficient swap drive space, and a strong parallel CPU which you might want to over clock to assist with speed problems. For most situations maps this size are unnecessary.

 

The second number dictates the size of each pixel, in this case the 20 goes with the default scale of milli-arcseconds or mas. You can change this scaling with the command:

 

0> mapunits arcsec

 

This command would then change the map scale to arcseconds. Be careful as any changes to your mapsize will require you to give the new scale in arcseconds. The three scales used are mas, arcsec, and arcmin. As time passes you will determine what scales and sizes are appropriate for which maps.

 

 

Some Helpful Hints:

 

1. Know your sky! There are two websites available to help you understand the sky around your source:

 

NVSS - http://www.cv.nrao.edu/nvss/postage.shtml

The NRAO VLA Sky Survey (NVSS) is a 1.4 GHz continuum survey covering the entire sky north of -40 deg declination. A detailed description appears in the 1998 May issue of The Astronomical Journal (Condon, J. J., Cotton, W. D., Greisen, E. W., Yin, Q. F., Perley, R. A., Taylor, G. B., & Broderick, J. J. 1998, AJ, 115, 1693). This survey was done primarily in the D Array of the VLA, thus resolution is not terrific, but with 1.8 million sources in the survey, most prominent sources in the sky are included.

 

FIRST - http://third.ucllnl.org/cgi-bin/firstcutout

Faint Images of the Radio Sky at Twenty-cm (FIRST) is a project designed to produce the radio equivalent of the Palomar Observatory Sky Survey over 9,000 square degrees of the North Galactic Cap. The survey uses the VLA in its B-configuration, and acquires 3-minute snapshots covering a hexagonal grid using 2×7 3-MHz frequency channels centered at 1365 and 1435 MHz. The data are edited, self-calibrated, mapped, and CLEANed using an automated pipeline. These maps have 1.8" pixels, a typical rms of 0.15 mJy, and a resolution of 5". At the 1 mJy source detection threshold, there are ~90 sources per square degree, ~35% of which have resolved structure on scales from 2-30".

 

It is important to know what sources are in the sky around your source as they will add confusion to your image and make it more difficult to selfcal all of your phases and amplitudes. You will often be able to notice a “rogue” source in your image as a sharp line extending from off screen into your image. You must expand your map or shift it to include it into your image. More on shifting will be added below.

 

If for some reason FIRST or NVSS images are insufficient/unavailable to locate sources around your primary source, an excellent and quick strategy is to raid the NRAO archive at https://archive.nrao.edu/archive/advquery.jsp and dig up an old L Band C/D array snapshot on your source or one nearby. A quick AIPS calibration and subsequent clean/selfcal in DIFMAP should reveal all nearby sources without having to make ludicrously large maps. Remember, the speed at which you can image is exponentially related to your mapsize. C and D array images can have extremely large pixel sizes, allowing you to fill the 30' primary beam size of an L Band image with a 2048 pixel mapsize.

 

2. The following is a list of information about different arrays/bands and the VLA

 

Band Wavelength (cm) Frequency (GHz) Primary Beam (arcmin) Maximum Resolution (arcsec)

 

P 90 0.30 – 0.34 150 6.0

L 20 1.34 – 1.73 30 1.4

C 6 4.5 – 5.0 9 0.4

X 3.6 8.0 – 8.8 5.4 0.24

U 2 14.4 – 15.4 3 0.14

K 1.3 22 – 24 2 0.08

Q 0.7 40 – 50 1 0.05

 

 

As a general rule of thumb, any time you are looking at an L band image, you will have a relatively large FOV, namely a half degree on the sky. Add to this that radio sources tend to be brighter at low frequencies and you have a very high probability to have sources within your FOV. You should refer back to the first helpful hint as to how to deal with these. When you move to C Band the probability decreases significantly, and for X Band and beyond you will almost always not have to worry about “rogue” sources in your FOV, though it does sometimes happen.

 

For the first time DIFMAP user the following initial mapsizes are recommended:

 

Band Array Mapsize

 

L A 1024,100

C A 1024,40

X A 1024,25

U A 1024,15

K A 1024,10

 

Some sources can be up to, and larger than, 2'. In this case none of these maps will give you a sufficient FOV to work with the whole structure. Again, these starting values are for the initial inspections of the map. They may, and often will, increase as the calibration continues. Though it is rarely the case, be certain not to expand the map far beyond the primary beam of the VLA as this is CPU and RAM being wasted.

 

3. When choosing a cell size, DIFMAP will complain if the size is too coarse (less than 2.5 pixels across the synthesized beam). It is recommended when making your map to have at the least 3 pixels across the beam, but is far better with 5-6 pixels per beam. You can determine this by using the mapl command:

 

0>mapl

Inverting map and beam

Estimated beam: bmin=198.5 mas, bmaj=248.2 mas, bpa=63.01 degrees

Estimated noise=0.669814 mJy/beam.

 

For this case you should look at the bmin value. The largest pixel size I would use would be 50 mas/pixel. A better situation would be 20-30 mas/pixel. If you move your pixels to a far too coarse size, DIFMAP will give you the following error:

 

0>mapl

Inverting map and beam

Your choice of large map pixels excluded 1% of the data.

The y-axis pixel size should ideally be below 53.32 milli-arcsec

Estimated beam: bmin=205.5 mas, bmaj=257.8 mas, bpa=52.17 degrees

Estimated noise=0.677087 mJy/beam.

 

This tells you that 53.32 mas should be the largest pixel you make, but be wary as using coarse pixel sizes may hinder the ability of the program to accurately place model points.

 

6. Selfcalibration and Imaging

 

With all of the above completed, it is now time to take the first look at our source. At the prompt type:

 

0> mapl

 

This command generates an image with an FFT on the uv-data. You should now see one of two things.

 

1. The source is not well centered in your map! At this point you should choose a larger mapsize, even one with a coarse grain to find your source. As an example:

 

0> mapsize 2048,1400

0> mapl

 

 

 

 

Though we have lost some data for having pixels far too coarse, this is temporary. The new view shows us that the source is located far to the north east.

 

 

 

 

 

 

 

At this point we want to shift the source to the center of the map. First, we must get the location of the source from the center of the map. Press “z” and make a small box around your source. This will zoom you in to the size of the box you specified. Your image should look something like the image on the right. Move your cursor over the strongest portion (here the brightest white portion) and hit “v”. This will give you the information about the pixel your cursor is on.

 

Pixel value at x=2e+05 y=1.99e+05 (mas) is 0.0415 Jy/beam

RA = 10 36 26.914, Dec = +13 26 50.800 (2000.0)

 

Using that the pixel value is 200000,199000 from the center, type:

 

0> shift -200000,-199000

 

This will move the source to the center of your map.

 

 

 

 

 

 

2. The other possibility is the same possibility you now face, namely the source is already at the center of your map and looks like the following:

 

 

 

 

 

 

At this point you will want to hit “c” to add color to the map and then hit “l” to refresh it. You should now see a lot of color. The rest of this tutorial will be done in black and white, and both color and grey scale work fine. At this point you should hit “f” to fiddle with the brightness and the contours on the map. Each time you fiddle you need to refresh with “l”. Different brightnesses and contrasts will be achieved by having your cursor on different parts of the map. Continue this process until you get something that looks similar to the following:

 

Now we have just the core of our source visible in the dirty map. If your source is very pointy, using a model to fit the source might be best. However, for this tutorial, we will use clean boxes.

 

To make your first clean box, zoom in tightly around the core and make a box around the brightest region by left clicking your mouse at one corner of the box and then move your mouse until you define a sufficient box, finishing the operation with another left click.

 

Defining the correct sized clean box is a bit of an art form that will come over time. Practice using different clean boxes and monitoring the effects.

 

With our first clean box made we can begin to clean the dirty map.

 

6.1 Clean Boxes and Selfcalibration

 

Clean boxes define the regions where clean components can be put by the program when we ask it to clean the image. The collection of clean components that we find will become our model and be used to self-calibrate the data. Because we will use the clean components for calibration, we don't want clean components to go just anywhere in the image. We only want clean components in regions we trust as being real parts of the source.

 

As we clean the image, the brightest parts will be cleaned away (including their nasty beam effects) leaving the less bright parts of the source. As these parts become more prominent, we will put clean boxes around them as well. Eventually, we will have a complete model of the source.

 

When you begin your clean you should start slowly as this will assure that the components will be subtracted gently from the image. Lets start with:

 

0> clean 250,0.03

 

This command will clean 250 components with a loop gain of 0.03, or 3%. The clean components are chosen at whatever is the brightest point in the image. The loop gain determines how strong the subtracted clean component will be. At a loop gain of 0.03, the subtracted clean component will be 1/33^rd the strength of the brightest point in the image for that iteration.

 

HINT:

 

Often times it is efficient to use a command like:

 

0> clean 250,0.03;mapl

 

This command, with a colon in the middle tells the program to run one task immediately after the other. This will display your map as soon as the clean is done. Also, once you have run clean once, any subsequent command, such as:

 

0> clean

 

 

Will result in the same clean as the previous one. You do not need to specify the number of iterations or gain.

 

 

 

 

 

 

 

Now that we have completed our first clean, lets look at the image. By pressing “m” while in the image I can display the model components, represented in the window by small green crosses.

 

Our tight box on the core seems to have cleaned most of the core flux and that flux is nicely centered on the core.

 

 

 

 

 

 

 

 

If we had not centered and kept a tight window, our image might look like the following:

 

In this image the clean box was too large and thus cleaned some flux we might not trust. The bottom right most model component is red, meaning it is a negative component. This and the other outliers should be removed to improve the model. Do not sweat if this happens. Simply move your cursor over the components you do not trust and press “r” to remove them. You can then press “d” to remove the nearest clean box, and redraw the box around the model you like.

 

 

 

 

Once you have achieved the model you trust for your core, it is time to make the first phase self-calibration. With this partial model of the source, DIFMAP will allow us to do a phase calibration against the remaining data. In order to do a phase selfcal, use the command:

 

0> selfcal

 

After you run this command you will want to look at the last portion of the output which:

 

Fit before self-cal, rms=0.029762Jy sigma=1.150860

Fit after self-cal, rms=0.021026Jy sigma=0.806191

 

Keep your eye on your rms! If it is going down it is a good indicator that you are on the right course.

Type “mapl” to see the new selfcaled image. You will see the the structure to the east and the south have disappeared. Do not worry, if they are real, they will come back. For now, we can see that there is more flux at the core to be cleaned, and it would be a good idea to clean the northern component as well. Place a new clean box around the northern component, and run another clean.

 

Over time, this image will also show a southern component you should believe as well. Once they are both cleaned, and you should run another selfcal. This is where you enter the “rinse and repeat” stage of DIFMAP. As you clean more structure, you can perform phase selfcals to bring more antenna phases in line. Over time you will produce a good image.

 

 

 

With this image, after about 4 rounds of cleaning and selfcal, you should see the image on the right. You can make out the two boxes I now have for the northern component, and one in the south. The core box is still there, but invisible in this image. The image now looks fairly noisy, and I am suspicious that there may be a “rogue” object in my primary beam. Lets check FIRST.

 

 

 

 

 

 

 

The FIRST image shows our source in the center. We can see the northern and southern structure we have detected, which is a good start. However, we can see that there are a few sources in the field we should watch out for! Time to make a bigger map and see if they are lingering in our image.

 

REMEMBER: The FIRST survey was done in L Band, so you should not always expect to see all of the same sources as you move to other bands! Many sources are steep spectrum and will be on par with the noise in your image.

 

Lets make a nice big map.

 

0> mapsize 2048,700

 

This map will cover the entire beam (9-10' for C band). As you can see, there is a source to the south west which has been marked by an over-sized clean box. You should make a nice tight clean box around this source and begin the clean/selfcal process again.

 

 

 

 

 

 

 

 

 

 

 

 

 

 

In this image you can see that there is also found a source to the north west to clean, but that should not be a surprise as it is evident in the FIRST image. At this point we can look at the bar below the image and see that there is very little more positive flux than negative flux in the image. After checking the image for negative/bad components and performing one last selfcal, we can choose to move to natural weighting to find more flux. If you desire maximum resolution, you may wish to stick with uniform weighting throughout. Play around with this to find out what works best in each situation.

To go natural weighting, remember the command:

 

0> uvweight 0,-2

 

Again it is time for rinse and repeat. Over time you will find another source a little north and east of your core as seen in the image to the right. Although there is positive flux to possibly clean, it is at this point in the mapping difficult to separate noise from structure, and so it is time to stop cleaning for a moment.

 

With a map we feel has been cleaned and phase selfcaled sufficiently, it is time to perform an amplitude selfcal.

 

WARNING: Amplitude selfcal is “risky”. You should only perform an amplitude selfcal if you have sufficient signal to noise in your data and a good model. Although phase calibration can also be risky, the process is smoother and thus minimizes the risk. Amplitude self-calibration has the power to “make” the data fit the model, and can thereby “freeze” false components into the data.

 

There are very few “rules of thumb” with respect to amplitude calibration. You need to play with your data and see what works best. To perform your first amplitude calibration, use the command:

 

0> selfcal true,true,1e6

 

This will perform an amplitude selfcal over 10⁶ minutes, otherwise known as your entire observation. At this point you should check your model and be certain you like the changes. If you do not like the effects of your amplitude selfcal, you can use the command:

 

0> uncal false,true

 

This will undo all amplitude corrections in your model.

 

For our map a million minute amplitude selfcal is just fine. Now we can look at our clean map to see the fruits of our labor.

 

NOTE: Checking the clean map as your go through this entire process is recommended, especially with more complicated sources to keep track of your progress.

 

In order to see our clean map we must first set the contour levels of the map with the loglevs command. We at Brandeis enjoy using levels which are multiples of root two. This means that every two contours you have doubled the flux.

 

0> loglevs 0.8,100,sqrt(2)

 

This command will place logarithmic levels starting at 0.8% and work up to 100% in intervals of root two. If you do not specify any loglevs, DIFMAP will defaultly go to:

 

0> loglevs 1,100,2

 

With the loglevs set, we need one more command to see the clean map:

 

0> mapl cln

 

We can now look at the cleaned map. We can see all three sources surrounding our source as well as the nice structure of the source itself.

 

With this done we need to clear our model and start again to try to bring our phases and amplitudes in line even more. Do this with:

 

0> clrmod true,true,true

 

This will clear the entire model and set you back to the beginning except you will keep your phase and amplitude calibrations. You may wish to reset your weighting back to uniform with:

 

0> uvweight 2,0

 

And begin your work again...

 

The idea is to continue doing this until you cannot improve the image in any way. Once you have performed your final pass, it is time to make the final image.

6.2 Clean Boxes and Selfcalibration

 

In order to make our final image, clear your model and make a map which is significantly larger than your source. For example, if your source and all other sources are visible within a mapsize of 1024,30, you may wish to make a map at the scale of 2048,30 in order to pick up model components that may be further away.

 

Open your new map and make one large clean box which covers 75-80% of your image. You can save your windows with the command:

 

0> wwin YOUR_WINDOWS_FILE_NAME

 

and read them back in with:

 

0> rwin YOUR_WINDOWS_FILE_NAME

 

With your new large box in place, we will perform a large clean. This will allow DIFMAP to put clean components essentially anywhere. With all of the phases and amplitudes lined up as best possible, this will not be as haphazard as it would have been if we did it as our first process. Try:

 

0> clean 10000,0.05

 

If at the end of this clean cycle there is still positive flux being cleaned, clean ever more components until positive flux is no longer being cleaned. At this time, switch over to natural weighting if you like, and clean even more.

 

Alternatively, you can use one of Dan Homan's scripts to do the final clean process. You can find them at http://personal.denison.edu/~homand/final_clean.

 

Take a look at your clean map with a mapl clean. If you would like to save your image as a postscript file, zoom to the scale you like with your image and press “k” to keep this framing. At the command prompt kill the color with:

 

0> mapcolor none

 

then prepare the file:

 

0> device file.ps /ps

 

if you are happy, run:

 

0> mapl cln

 

This will now read out and be saved as “file.px” in the directory you started DIFMAP in, to be examined or printed later on.

 

You can get back to color with:

 

0> mapcolor color

 

or

 

0> mapcolor grey

 

and then:

 

0> device /xw

 

lastly, save your work:

 

0> save file

0> exit