A short outline of the algorithm: start with an initial model. For each pointing, multiply the model with the primary beam pattern; convolve with the synthesized beam of that pointing. Subtract the observation and calculate the chi^2 of the residual. Then apply the formulae of the Maximum Entropy Method, in which unobserved data are constrained by a "reference map". To find a step, use Lagrange Multipliers and a Newton-Raphson iterative method. Drive the sum of the chi^2 values of all maps to a value such that on average the residual deviates by one sigma from the prediction and such that the total flux converges toward a wanted value.
The many keywords below come in five groups:
- in, beam, pbpar, maxcorr, rms -- specify interferometer inputs - region#, planes, initial, reference, flux -- set up constraints - out, mode, beamsize, save, saveiter -- specify outputs - maxiters, tol, measure, lm0, meslev -- control parametersObligatory keywords are: in, beam, rms and out; the rest is for fine tuning.
Inputs for interferometer data are: a list of miriad datasets with input observations, a list of the corresponding synthesized beams, an optional list of primary beam descriptions and a list of values for the rms noise in each input map.
For each of the input maps it is possible to select a region outside of which the model is set to zero (see the region and reference keywords). A list of common planes can be given to simplify the selection of channels. If datasets are not on the same frequency grid, use miriad task regrid to make them fit.
The output is a miriad dataset with the model. It is also possible (and advisable) to write output datasets containing the residuals of the input map, a combined residual (without primary beam corrections), the model convolved with a "clean beam", and a final map (the sum of the primary beam corrected residual and the convolved model). The reference coordinate (crval) of the output maps is that of the first in the list of input maps. Note: the mask of the model (outside which it is zero) is only written to the mask item of the 'convolved' output map.
Calculating the secondary outputs (residuals, convolved and final map) is equivalent to using 'restor' after a 'clean' operation. They can also be calculated from a read-in model map (see the mode keyword).
Most intermediate results (like the evolving model and residuals) can be written by using the save= keyword.
The program has three other independent uses: 1) to write a linear mosaic of the primary beam corrected, variance-weighted input maps; 2) to write a noise map, i.e. the harmonic mean of the variances; 3) to write primary-beam-corrected input maps. A choice is made using the mode= keyword.
An example of standard simple usage is:
mosaic in=map* beam=beam* rms=1,... flux=-10 out=out mode=m,f,r,c
TAKE NOTE OF THE FOLLOWING:
Beware: the value given for the rms is of paramount importance in determining convergence. In fact, it is the parameter that has the most influence on the reliability of the result. If it is too low or too high convergence may not occur.
A current limitation is that the crvals of the input maps must be an integer number of pixels offset from each other. The program checks this and tells you what input to use for the invert keyword offset= to get it right. The reason for this limitation is that mosaic does not regrid the maps.
For large fields (larger than about one degree) geometrical projection effects may become important. Maps made with invert will not be on the same tangential plane. Usually this can be ignored. Strictly speaking, however, they should be reprojected (using the miriad task regrid). Unfortunately, miriad does not differentiate between pointing center and projection center, and this causes trouble. So, after using regrid, you need to fix the crval and crpix of the reprojected datasets so that they indicate the pointing center of the observation.
As the model is forced to be non-zero everywhere inside the selected region, it can also be non-zero where there is no actual emission. Thus, it may be even more important than in CLEAN to define a region where the signal is expected. Use the region#= or the reference= keyword (see there for details). One way to create a region is to first run mosaic with mode=m,c, without region selection; then take the resulting convolved model, and fiddle a bit with it; next use it as input to the reference keyword, making use of the clipping option. A choice for the region that is too small will result in missed structure in the model. Or, the bias flux that is spread around in the map may become too high. A choice for the region that is too large can result in "structure" being found in the noise.
It is furthermore important to put in a good value for the total flux (if only as initial estimate). If you don't have this, don't be surprised if the answer looks bad or convergence does not occur. You might try to iterate by changing the initial input value until it does not vary much from the value at the last iteration.
A good reference map improves convergence. Without it, in principle the final result should be the same, but you may have to wait longer if the source structure is complex. The reference map gives the algorithm reasonable data on spatial frequencies that were not observed. There are several possible ways to create a reference map (or an initial estimate):
a) a flat map (default, usually good enough and preferable) b) the sum of the primary beam patterns (for sources whose diameter is comparable to the primary beam) (often unreliable) c) a single-dish map gridded to the interferometer gridspacing d) a map of the same source at a different frequency e) the result of a previous or lower-resolution observationThe latter three cases require you to specify the reference map as input. If any a-priori information is present, use this to create an initial estimate and reference map.
Officially, miriad has no header variable that represents the pointing center. The uv-variables "obsra" and "obsdec" come closest. They are defined as the phase center at the epoch of the observation. To get the phase center at the epoch of the map requires a precessional correction, which is not done (yet?). Thus, "crval1, crval2" are used instead as the coordinates of the pointing center. For miriad datasets this should almost alway be OK, but it is a non-guaranteed use of the variables. For maps imported from aips, check whether crval/crpix indeed give the pointing center. If "crval1, crval2" does not give the pointing center, the primary beam is offset and the algorithm will fail.
Be very cautious using this method. If you don't understand it, it is easy to obtain wrong output maps. Don't blame the program if you fail at the first try, fiddle with the parameters or ask for help (remember the need for an accurate rms).
- "" => tel - none => no primary beam correction - tel => use formula for telescope as described below if telescope not known, make a gaussian with fwhm from header item "pbfwhm" - gauss,fwhm => create a gaussian primary beam (cut off at 0.05) - airy,fwhm => create an airy disk - airyc,fwhm => create an airy disk, but only out to first zero - ... => repeat the last given value for remaining inputs If fwhm > 0 => make primary beam of given type with this fwhm (if restfreq>10 GHz units are arcsec, else arcmin) If fwhm = 0 => read fwhm from header item "pbfwhm"; if not present, get fwhm from "restfreq" and "telescop"If "pbpar=tel", header items "telescop" and "restfreq" are used to make the primary beam: a polynomial for the WSRT, VLA and ATCA, a gaussian for FST and BIMA; these functions are cut off at a level of about 0.05 (see keyword maxcorr). If "restfreq" is not present, but the type of the "freq" axis is "FREQ", then crval, crpix and cdelt are used to get the frequency of plane 1. Example 1: with two input maps, pbpar=gauss,120,gauss,120 makes gaussian primary beams with fwhm 120 arcsec. Example 2: pbpar=airy,120,tel,gauss,0,gauss,118 with four input maps makes an airy disk with fwhm 120 arcsec for the first, uses the telescope formula for the second, makes a gaussian with fwhm found from "restfreq" and "telescop" for the third and makes a gaussian with fwhm 118 arcsec for the fourth.