Monday, April 23, 2007

Teaching the Monkey to Dance: Scripting AIPS #2

In the first installment of this thrilling saga, I talked about how to name your AIPS scripts and load them into AIPS. Now, we're going to get into the nitty gritty details of how to write your own AIPS scripts.

In AIPS scripts, you have access to all the usual AIPS commands and tasks like listr, calib, and imagr. You also have access to some basic programming tools. To see exactly what you have in your toolbox, type help popsym (POPS is the name of the AIPS scripting language). I've replicated the contents of the popsym help files below for your enjoyment:

Help on POPSYM in AIPS version 31DEC04
Type: Symbols used in the POPS interpretive language


----Arithmetic expressions

+ A + B Add the expression A to B
- A - B Subtract the expression B from A
* A * B Multiply the expression A with B
/ A / B Divide the expression A by B
** A ** B Calculate A to the power B
( ) (A+B)*C Grouping expressions as desired
= A = B Store the value of B into A
, A = 3,5,4 Separator of elements in an array
~ A(i) ~ 1,2,3 Store values in A(i),A(i+1)...
(change only as many as on RHS)
: TO Equivalent to the verb TO
; Separator between AIPS statements

----Logical expressions

> A > B A greater than B
< A < B A less than B
= A = B A equal B (numeric or string)
>= A >= B A equal to or greater than B
<= A <= B A equal to or less than B
<> A <> B A not equal to B (numeric or string)
! A ! B A or B
& A & B A and B
^ ^ A not A

----String expressions

!! A !! B string = string A followed by string B
SUBSTR SUBSTR(A,i,j) string = chars i through j of string A
LENGTH LENGTH(A) position last non-blank in A
CHAR CHAR(A) convert number A to string
VALUE VALUE(A) convert string A to number

----Looping constructions






----Built-in functions

ATAN Arctangent (one argument)
ATAN2 Arctangent (two arguments)
COS Cosine (degrees)
SIN Sine (degrees)
TAN Tangent (degrees)
EXP Exponential
LN Log base e
LOG Log base 10
SQRT Square-root
MAX Maximum i.e. X = MAX (A, B)
MIN Minimum i.e. X = MIN (A, B)
MODULUS Root-square sum of two arguments
MOD(A,B) A - (A/B) * B i.e. remainder of A/B
CEIL(A) Lowest integer >= A
FLOOR(A) Highest integer <= A

----Procedure building verbs

PROC PV Begin building a procedure
PROCEDUR PV Begin building a procedure
LIST pV List a procedure
EDIT PV Edit a procedure
ENDEDIT PV End editing a procedure
ERASE PV Delete line(s) of a procedure
MODIFY PV Modify a line in a procedure
RETURN V Last statement in a procedure
FINISH PV End procedure building

----Variable declarations

SCALAR pV Declare scalars
ARRAY pV Declare arrays
STRING pV Declare strings

----Input/Output functions

PRINT V Print the following keyword value(s)
TYPE V Print the following keyword value(s)
READ V Read value(s) from terminal after # prompt

----Other information

CORE pV Amount of core left in POPS
COMPRESS PV Compress the core area, recovering lost space
and acquiring any new vocabulary
CLRTEMP V Clear the temp data array
DEBUG pV Debug: turns on compiler debug information
DUMP V Dump K array on terminal screen
SCRATCH PV Remove procedures in POPS
$ PV Makes rest of input line a comment

As you can see, you have most of the basic programming tools like if/then statements, for loops, mathematical functions, and comparison operators. The only non-useful bits are in the "Procedure building verbs" section. Some of the verbs in this section (list, edit, endedit, erase, and modify) are left-over from the bad old days and were used for writing procedures using AIPS built-in line editor. Trust me, you don't want to use this editor. Just write your scripts outside AIPS in your favorite text editor and load them into AIPS.

Now that you know what tools you have to work with, let's move on to how AIPS scripts are actually structured. In general, an AIPS script has two parts. The first part declares all the variables one is going to use in a script, while the second part is the actual script itself. Note that you can have more than one script (declaration plus code) inside the file you created in part 1. The convention is to group scripts that work together in one file. For example, my custom calibration routines contain a declaration section, a calibration script, and a script to nuke the current calibration. Start each section with proc scriptname, substituting the name of your script for scriptname and end each section with return; finish.

Looking at the pops language commands, you only three options for variables. (No advanced data structures for you!) Scalars are just numbers, while strings are just letters. Arrays can either be a collection of numbers or letters. To declare a scalar, you just say scalar anumber. Strings need to have their length declared (just like Fortran): string*8 aname. This creates an eight character string. Now declaring an array is slightly tricky. First, you decide what sort of array you want, then you append a number telling AIPS how big to make the array. To create an array of 15 strings, you say string*8 obj_name(15). Variable names are limited to about nine characters, so choose them wisely!

Now to remember what you were thinking when you wrote a script originally you can put an asterisk (*) at the beginning of any line to indicate that the entire line is a comment. You can also put a dollar sign ($) anywhere in a line to indicate that the rest of the line is a comment.

Note that AIPS does NOT pay attention to the first line of a script file, so always start your files with a comment (put a * on the first line). In addition, AIPS gets cranky if you have any tab characters in your script. (This may just be an issue relating to how my editor of choice, emacs, stores files.) If everything looks fine in your script and AIPS still won't read it in, look for tabs. Usually if you move the cursors over the line of code with your arrow keys you can see where the tabs are because the arrow keys move the cursor more than one space.

Well, that's enough for now. Next time, I'll go through a simple example demonstrating everything I talked about in this entry.

Read more!

Thursday, April 19, 2007

TV Trouble?

Having trouble with your TV? If you're on a funny network setup, the AIPS TV display can run slow or refuse to open. Try starting AIPS with aips tv=local and see if that fixes the problem.

Read more!

Tuesday, April 17, 2007

Self-Cal Wackiness

UPDATE (4/18/07): I think I've tracked the errors down to very bad amplitude solutions for antenna 1 for about an hour of one day. Still working on a solution. I'm tempted to flag antenna 1 for that time frame since it seems so self-contained, but that makes me a bit nervous.

UPDATE (4/23/07): I ended up flagging antenna 1 for the time where the solutions weren't converging and one scan on antenna 28. Now the data look fabulous. These antennas had phase errors in the phase calibrators during the relevant periods of time, so I'm feeling better about flagging them.

This isn't completely relevant to this blog because I'm doing it in Miriad, but I'm getting some funny artifacts when I amplitude+phase self-cal in Miriad. The phase self-cal goes fine, but I get these cheap carpet patterns when I amplitude+phase self-cal. The first image below is no cleaning, phase self-cal only. I can get rid of the swirls with cleaning. The second image is no cleaning, amplitude+phase self-cal. Weird! Any idea of the origin of this pattern? Has anyone else seen something like this?

I'm wondering if I can ignore the amplitude self-calibration and just do the phase self-calibration.....

Read more!

Monday, April 16, 2007

Flux Conservation in Self Calibration?

Any one know what self calibration really does to your flux? I need fluxes of sources that are mostly point-ish, and I like self calibrating because it makes them look more like circles and less like amorphous blobs. But I wonder if I can trust my fluxes after self cal? If I do only phase self cal, should I be ok?

I found in someone's Ph.D. thesis that your fluxes are no longer "absolute" after self calibration. I guess because you're no longer comparing everything with 3C286, you're comparing your target source with itself.

Anyway, when i just do an imean on my images, the total flux density doesn't change a huge amount. Maybe 20% (well, sometimes I get factors of 2, but I avoid those cases). I don't really want to be introducing 20% errors into my flux measurements, though....
Maybe I shouldn't be self caling if I want robust fluxes?

Anyone really understand what cparm(2) does in calib? Is it only important for A&P self cal, or P self cal too?

Read more!

Sunday, April 15, 2007

Self Cal according to Claudia (as inspired by Crystal Brogan)

I'm pasting Claudia's notes below, about how she self calibrated her SMA data. I always like to see how people do their self cals! I'll probably post later some more details about how I personally do it.

Clean in IMAGR
To get a circle (as opposed to a square):
to get another one, keep hitting c until you get a circle (may take a few times)
When cleaning interactively, in order not to have to redraw the boxes every time:

to create a boxfile, set a filename in oboxfile parameter in imagr (with the logical name ahead of it, same as when you're reading in a fits file)

Then, to use it once you've created in, set boxfile=oboxfile. This will read in the file, allow you to modify it, then spit it out again (with any modifications). Can leave this set in imagr through all the steps that follow

REBOX _not_ TVBOX to add boxes when you're working with a file--otherwise it will delete everything you already have

Stop cleaning when the image looks about uniform. A good way to decide when to stop is when the min and max (see pasted in example below) are about equal.

IMAGR2: I Polarization model processed
IMAGR2: Field 1 min = -15.4 MilliJy,max = 16.1 MilliJy

In a self-cal cycle, can tell that your image is improving if it includes more visibilities than the last iteration.

NB: The best self-cal model (the model is the image you're making of the UV dataset that comes out of a calib run) is the one that contains the most flux--but beware of boxing dubious sources, since anything you put into the model will create a source around it. If a dubious source gets weaker when you box another, more probably real source, that's a good indication
that it's not real.

Herringbone pattern can be due to amplitude errors, poor UV coverage, the beam shape being poorly defined.

i) To get started, make an image of your original UV database (say BASFIT from UVLSF)and clean it as described above--probably only with on box.

ii) Task 'calib'
Do phase-only first: CALIB SOLMODE ='P '/Solution mode
Feed in as your model (get2n) the map made in (i)
Pick a short solution interval for phase, since phase varies quickly as a function of the atmosphere (say 2
CALIB SOLINT = 2.00 /Soln. inter. (min)
CALIB APARM(1) = 4 /Min. no antennas
CALIB APARM(7) = 3.0 /SNR cutoff

Can use a lower S/N cutoff than the default b/c phase is constrained by phase closure.
set OUTCLASS for calib to 'calib' to distinguish the UV files that will come outalso add a _P to the filename to
mark it as phase-cal only


iii) image your _P.CALIB UV file that 'calib' has produced, using imagr

iv) Run calib again. Input is your ORIGINAL UV data. The map (in2name) is the map you made in (iii). see
examples below.
116 CALIB /TIMERANG = beginning to end
**this is very important. This makes the self-calibration insensitive to errors early on--eg boxing sources that it later becomes apparent you shouldn't have boxed. and it retains flexibility. And means that your original UV data file will accumulate SN tables, and you can correctly copy the last version over to the line data and apply it there. DON'T do iterative self-cal. This is most important with low S/N and poor UV coverage data--certainly the case for SMA.

Repeat steps (iii) and (iv) until the phase is as good as it can get.

Once it is:

v) do an amplitude calibration
NB: amplitude calibration is inherently much less constrained than phase calibration. Never do an amplitude calibration until the phases are fixed as best they can be.
This time, use your last, most correct, .CALIB UV file as an input to task CALIB
set CALIB APARM(7) = /SNR cutoff parameter back to =5 (since amplitude is much less constrained)

Unlike phase, which varies quickly as a function of atmosphere, amplitude should vary slowly, and the variation should be primarily due to instrumental effects
Want to solve for amplitude only once per scan (scan=length of time on source)
(for s255n data, 20 minutes)
set SOLINT=scan length

CALIB SOLMODE ='A&P '/Solution mode

vi) After amplitude calibration, do uvplt as a reality check. If you have crazy amplitude errors, specific baselines may jump way up. (Remember that you're actually changing the amplitudes with the amplitude calibration)

vii) SMA weights are based on tsys. Will (particularly with low S/N, poor UV coverage) get a better image if you apply the weights, rather than pretending all antennas have same tsys (which at SMA they emphatically don't).

So, make an image by running imagr, with your last calib phase file as an input, and applying the SN table (which will be the SN table from the amplitude calibration)
To do this:
(for example, if this is your most last phase cal iteration)
149 IMAGR DOCALIB = 2.00000E+00

Read more!

Thursday, April 12, 2007

Bandwidth Smearing

Once your radio data is calibrated you can SPLIT off the source and apply these calibrations. During SPLIT you'll want to average together as many channels as possible. The number of channels you can average together will be limited by bandwidth smearing which causes sources to smear out in the radial direction. Bandwidth smearing becomes significant when the product of (\frac{\delta \nu}{\nu})*(# of beams from the image center) reaches unity where \delta \nu is your channel size and \nu is your total bandwidth. So, if you have a really small beam you won't be able to average together very many channels during SPLIT if you want to avoid smearing sources that are far from the center of your field. For VLA data the correlator doesn't give you very many channels so this doesn't matter as much as with GMRT data which always gives 128 channels. I find that with 610 MHz GMRT data I can only average every 2 channels when I SPLIT off the source if I want to avoid radial smearing out to the edge of my primary beam (HPBW of 0.7 degrees). Your source dataset will most likely have greater than one channel when you're ready to run IMAGR but you only want one image in the end so you'll need to set NCHAV and CHINC to be the total number of channels in your dataset. IMAGR is smart and it knows how to combine these channels during the imaging process so that bandwidth smearing isn't an issue. (So, it's actually not imperative that you average any channels together when you SPLIT off your source data but it does help IMAGR to run faster if you do.)

Read more!

Wednesday, April 11, 2007

Teaching the Monkey to Dance: Scripting AIPS #1

Is your wrist tired from typing and re-typing all those AIPS commands and parameters? Then this tutorial is for you! My goal is to relay my hard-won AIPS scripting knowledge to the rest of the world. Enjoy! For the first post, I'm show you how to name your scripts and get them into AIPS.

First, we're going to learn how to get AIPS to figure out where your scripts are. Fire up AIPS in the directory where you're going to put your test scripts. To tell AIPS where to look for your scripts set version="PWD". Now AIPS will look at the directory you started it up to for your scripts. Now we need something to call your script. Unfortunately, you can't just call the script mywonderfulsupergreataipsscript. First, the name needs to be fairly short (try for less than 9 characters). (Note that it also needs to start with a letter NOT a number. Thanks for the tip, Emily!) Second, you need to take the obvious step (or at least obvious if you're AIPS) of appending your AIPS user number in hex to the end of the filename. Fortunately, the task EHEX in AIPS can help you with that. Put in your user number and it will return a four digit number that corresponds to your user number in hex. Append the LAST THREE digits of the number to your file name. Third, you need to have the filename all in uppercase letters because AIPS only understands uppercase.

For example, my user number is 333. EHEX tells me that this is the equivalent of the hex value 0099. Therefore, I call my first script HELLO.099.

Now comes the fun part, what to put in your new script? Below I give some sample code that will tell AIPS to say hi. Type this into the HELLO.099 file.

proc sayhello

print 'Hello'

return; finish

Now in AIPS, type run hello. You should see the contents of HELLO.099 printed out on the console. Then type sayhello and AIPS should say 'HELLO'. If so, congratulations! If not, go back and check on your script. If you're still having trouble, leave a comment.

Next time: We're going to get into the nitty-gritty of actually writing a script.

Read more!

HTML Tips for Code

The pre HTML tag might be helpful. You can't use it in the comments, but you can use it in posts. The following HTML code should generate a cross made out of the word stuff:

stuff stuff

Read more!

Speeding Up Self Cal (a tiny bit)

When doing self calibration, you have to do lots of iterations of calib and then imagr. IMAGR has been taking forever for me!

Just realized that imaging goes a lot faster if you use the calibrated UV files which are outputted from CALIB, instead of using the orignal uv file and applying the SN table in IMAGR. Every little bit of speed helps!

Read more!

Imaging with Multiple Fields

There are two reasons you might want to image with multiple facets.
1) To correct for 3D effects, and image degradation that goes along with that.
2) Because there's a really bright source far away from the center of the field that's messing up your image, but you don't want to image all the way out to it. It would be so much nicer if you could just put a small field on it!

Here's some of the tricks of imaging with multiple fields:

You use SETFC to figure out where to put the facets. I think SETFC is really a pretty great program.
The first time you run it, if you set cellsize and imsize to 0, SETFC will actually recommend a pixel size and image size to you! Like this:

1 5 23-MAR-2007 14:08:03 SETFC Task SETFC (release of 31DEC05) begins
1 3 23-MAR-2007 14:08:03 SETFC Found NGC3184 SPLIT Seq 1 Disk: 1 in slot 6
1 4 23-MAR-2007 14:08:20 SETFC SETCEL: recommends IMSIZE 932 CELLSIZE 1.26446
1 4 23-MAR-2007 14:08:20 SETFC SETCEL: returns IMSIZE 1024 CELLSIZE 1.17417
1 2 23-MAR-2007 14:08:20 SETFC ZTXOP2: using translated file name =
1 2 23-MAR-2007 14:08:20 SETFC ZTXOP2: TEST
1 4 23-MAR-2007 14:08:20 SETFC FLYEYE added 7 fields to BOXFILE to .283 deg
1 3 23-MAR-2007 14:08:20 SETFC Searching catalog between .28 and 1.50 degrees radius
1 2 23-MAR-2007 14:08:20 SETFC ZTXOP2: using translated file name =

I often rerun SETFC with a rounded off version of cellsize (like 1.2 in this case) so I can remember it more easily. CHANGE AS OF 5/09: I"ve really decided that the CELLSIZE suggested by SETFC is far too small. It often just barely Nyquist samples your source-- often giving only 2 or 3 pixels across it. I now try to use CELLSIZEs which are 2 or 3 times smaller than that suggested by SETFC. Also, I don't like imaging over fields bigger than 2048 because than IMAGR averages many pixels together to display the field, and it's hard to tell what's going on. So, I often have imsize of 2048 and a cellsize that gives ~5--6 pixels across a point source.

Now, the other cool thing about SETFC is that you can set an inner tiling of facets, and then you can have some external fields on bright sources (SETFC queries the NVSS database to know where there are bright sources around your source). If you want both internal and outlying fields, you're input will look like this:

AIPS 1: SETFC: Task to make a BOXFILE for input to IMAGR
AIPS 1: Adverbs Values Comments
AIPS 1: ----------------------------------------------------------------
AIPS 1: INNAME 'NGC' UV dataset name (name)
AIPS 1: INCLASS 'SPLIT' UV dataset name (class)
AIPS 1: INSEQ 1 UV dataset name (seq. #)
AIPS 1: INDISK 1 Disk drive #
AIPS 1: SOURCES *all ' ' Source selected
AIPS 1: BCOUNT 1 First field number to use
AIPS 1: disk file to write to (the
AIPS 1: input BOXFILE for IMAGR)
AIPS 1: CELLSIZE 1.2 1.2 (X,Y) size of grid in asec
AIPS 1: IMSIZE 1024 1024 field size
AIPS 1: SHIFT 0 0 Position shift (RA,Dec) asec
AIPS 1: for all fields
AIPS 1: FLUX 0 Minimum component flux =
AIPS 1: (source * beam)
AIPS 1: BPARM .3 10 (1) Inner region radius (deg)
AIPS 1: 0 1.2 (2) Field overlap (pixels)
AIPS 1: 5.000E-04 256 (3) Factor to scale NVSS
AIPS 1: *rest 0 fluxes, 0 -> 1
AIPS 1: (4) Radius NVSS search (deg)
AIPS 1: (5) Flux limit in NVSS (Jy)
AIPS 1: (6) IMSIZE for NVSS fields
AIPS 1: (7) IMSIZE for Sun fields
AIPS 1: (8) Write Clean boxes for
AIPS 1: NVSS fields
AIPS 1: (9) Maximum allowed phase
AIPS 1: error in imaging
AIPS 1: (10) Points per beaam
AIPS 1: PBPARM *all 0 Beam parameters:
AIPS 1: (1) Cutoff; (2) Use (3)-(7)
AIPS 1: (3)-(7) Beam shape parms
AIPS 1: NVSS input file name
AIPS 1: ' ' => AIPS provided.

This tiles 7 "big" (1024x1024) fields in the center of my image (within 0.3 degrees).
And then it puts smaller (256x256) fields on any NVSS sources with flux greater than 5e-4 Jy which are located between 0.3 and 1.2 degrees of the center of my field. In this case, that gave 31 external fields.

As Joe Lazio suggests, you might want to do a quick imaging run of all these external fields to see which ones actually contain sources, and delete from the BOX file any fields that don't really contain much of anything in your images. If you do delete some fields, you'll have to renumber the renaming ones or I think IMAGR will balk at you.

When you're ready to image with IMAGR, here are a few important things to keep in mind:
IMSIZE-- is the minimum field size! So in the above case, I'd set it to 256. IMAGR is smart enough to know to make the bigger internal fields 1024x1024 (because the BOX file tells it to).
NFIELD-- set this to the total number of fields in your BOX file.
BOXFILE-- set this to the file that came out of SETFC
OBOXFILE-- IMAGR will be smart enough to print out the field locations/sizes, AND any clean boxes to this file. Hurrah!

When you do start imaging, if you have DOTV = 1, IMAGR will show you one field at a time, and I think it tries to show you which ever field has the most flux in it (or the most flux remaining in it, after cleaning.)

Read more!

Tuesday, April 10, 2007

tvflg outputs

Quicky: How do I keep TVFLG from making those stupid .TVFLG. files whenever I use it?

Read more!

Monday, April 2, 2007

Number of Channels = One Less Than I Expect

4/2/07: When I write my observe file, I ask for 8 channels. And when I get the data and FILLM, there are always 7 channels there.
Now, I have a new data set in which I asked for 4 channesl, and when I FILLM I am again missing a channel, there are only 3! This is a pretty significant loss of signal; I'd like that channel back.

I think I'm using pretty standard, default parameters in FILLM. Has anyone ever found a way around this problem? (I'm sure it's not a big deal if you're doing bona fide spectral line data).

UPDATE as of 4/14/07: After talking to my advisor Eric, it appears that Channel 0 actually counts as one of your channels!!! So when you ask for a correlator setup of 4 channels, you actually get 3 channels of original data, and then one of the channels is taken up by the average of these, CH 0. This is really one of the stupider things I've heard yet about radio interferometry. And let this be a warning to you, don't ever ever use the 2-channel correlator setup at the VLA!!!! I guess it's the same as doing continuum observations, but with less bandwidth!

Read more!

Using Calibrator Models

Wisdom from Stefanie Muehle:

BTW, I now got a nice calibration of my spectal line data using the models as recommended in the AIPS cookbook. The thing you have to get your head around is that when you apply a model, you don't _expect_ the phases to be zero, because you are dropping the default assumption that the calibrator is a point source. The clean components and the visibility plots are not much help in checking whether your calibration is ok, either. So, here's the tip from Michael Rupen, one of the gurus at NRAO:
The structure of the phases will depend in detail on the model.  To check whether you're getting something reasonable, try the following:

* If you restrict the uv-range in UVPLT to match the original (no model)  recommendations, you should still get flat (zero) phases.  Do you?

* You can check the phases predicted by the model by using UVSUB with opcode 'MODL' (this may require splitting off 3C xxx to make single source  files).  Try that and see whether the predicted phases match what you observe.

If the post-calibration phases do NOT match the model, there is indeed something wrong, and we'll work with you to track it down.  I hope and expect this is not the case though... :}

So, you basically first split off the calibrators into single-source files (SPLIT), then you substitute in each file the observed phases with those derived from the model images provided by AIPS (UVSUB with optype = 'modl'). When you then plot the phases of these files, they should look like the UVPLT of your calibrated data, only without the noise. At U-band (2cm) the model phases of 0137+331 and 0542+498 are so close to zero that the difference between the model and a point source is lost in the noise. But for 1331+305, the effect was quite spectacular. The model predicts phases up to |\phi| = 60\degr! That doesn't look like one should "force" the phases down to around zero by assuming a point source.

The catch is that I don't know whether the calibration with a model affects the polarization calibration. PCAL should be ok, since we are using the phase calibrator, which is assumed to be a point source. But RLDIF requires the polarization calibrators 0521+166 and 1331+305, for which the phases (RCP and LCP) will be substantially different from zero.

Read more!

Getting TELL to Work

TELL sounds like a really great idea; For example, I might accidentally turn off the TV in the middle of a 6 hour long IMAGR session, and I could turn it back on by telling AIPS dotv =1!

BUT--does anyone know how to get TELL to work? Unfortunately, it appears that you need a prompt to be able to use TELL, and the long routines that one runs usually withold the prompt.

The way I see it, I think there are two ways to get a prompt, and to be able to communicate with AIPS and TELL it something.

1) dowait = false
haven't actually tried this, but i think that it will give you the prompt back after you start IMAGR. However, from quick glances at the help pages, this may limit how IMAGR interacts with the user.


2) open up a new AIPS window and log in with the same user ID. But when I did this, it didn't recognize that I was running IMAGR. Maybe there is something fancy to make AIPS recognize other windows of the same user?

Read more!