This procedure creates a smoothed map of some function on the sky,
sampled at discrete points. It relies on adxy and gcirc from the
IDL astronomy user's library to handle all of the hard
spherical geometry stuff.
map: A blank map and fits header. Use map_init.pro to create
x: The x location (sky coordinates in degrees) of the data
y: The y location (sky coordinates in degrees) of the data
val: The value at each (x,y) point
dval: The 1-sigma error at each point.
fwhm: The fwhm of the smoothing kernel, in degrees. Defaults to
1/100 of the map size
truncate: The radius at which to truncate the smoothing. This can
considerably speed up execution time. See PROCEDURE for
details. Defaults to 2 * fwhm.
emap: Set to a named variable to hold the estimated variance
On output, the map.map array will be populated with the smoothed
map. Pixels which are further than (truncate) from the nearest
source are set to NAN.
This function uses weighted mean smoothing with a Gaussian
V(x,y) = sum(w_i * val_i) / sum(w_i)
w_i = 1/dval_i^2 * exp[-((x-x_i)^2 + (y-y_i)^2) / 2 sigma^2]
w_i truncated to 0 at (x - x_i)^2 + (y - y_i)^2 > truncate^2
The variance map is given by
dV(x,y) = sum(w_i * dval_i^2) / sum(w_i)^2
In other words, it is the interpolated value of the variance, divided
by the "effective" number of points which contribute to each pixel.
The procedure only calculates weights out to the truncation radius
to speed things up. The procedure takes O(stamp_size * n_data)
operations, instead of O(map_size * n_data).
March 2010: Written by Chris Beaumont. Meant to replace
Dec 2010: Added /SAMPLE to REBIN calls. cnb.