# default values may be fine, but this is how to change them
set gridlo1 = 0.0
set gridhi1 = 2.5
set maxlev = 4
# to illustrate deeper recursion
set imglev = 6
# for convenience with gnuplot
set plotmode = 2
# specify the mass model
setlens 1 1
alpha 1 0 0 0 0 0 0 0 0 1
0 0 0 0 0 0 0 0 0 0
# find some sample images
findimg 0.2 0.0 sis1.img
findimg 0.2 0.3 sis2.img
# write grid and critical curve(s) fot plotting
plotgrid sis.grid
plotcrit sis.crit
#
quit
Here are the image (left) and source (right) planes. The grids are shown in blue, while the critical curve and caustic are shown in red . (Note that gravlens/lensmodel does not actually make plots; it writes data files that you can plot with your favorite plotting software. I have used gnuplot here.)
To better understand this situation, let's turn it on its side and add a third dimension to show how the lens mapping takes the image plane grid and folds its over on itself to obtain thes source plane grid.
Let's zoom in near the critical curve to see the recursive subgridding. The subgridding near (1.2, 0) was done to find one of the images (see below).
Now let's find and plot some image positions. The source plane (right) shows two source positions, and the image plane (left) shows the corresponding sets of images. Although we specified a "singular" isothermal sphere, the code imposes a small core radius for numerical stability, so each source produces three images. The two central images are at basically the same place, and the blue point covers the green point.
set gridlo1 = 0.0
set gridhi1 = 2.5
set maxlev = 4
set imglev = 6
set plotmode = 2
#
setlens 1 1
alpha 1 0 0 0.5 0 0 0 0 0 1
0 0 0 0 0 0 0 0 0 0
#
findimg 0.0 0.4 sie1.img
findimg 0.1 0.2 sie2.img
findimg 0.18 0.27 sie3.img
#
plotgrid sie.grid
plotcrit sie.crit
#
quit
Here are the grids and critical curve and caustic for this example. This example illustrates one minor limitation of the gridding algorithm: the bottom of the critical curve is cut off slightly. The code find the critical curve by looking for tiles where the four vertices have different signs. If the critical curve juts into the edge of a tile without overlapping a vertex, the code will miss it. However, this should not be a problem for most analyses. If it is a concern, you can increase the resolution of the grid using the ngrid2 variables.
And here are some sample image configurations. As before, the code finds a central image for each multiply-imaged source, but the three central images are all basically on top of one another.
So far we have plotted the grid, critical curve and caustic, and image positions. There are many other things you can plot as well. For example, you can plot the time delay surface corresponding to source #2 (the blue points) above by adding the following line to the input file:
plottdel 0.1 0.2 -2.0 2.0 200 -2.0 2.0 200 sie2tdel.fits 3
In this example the code writes the data to a FITS file, which I have plotted using ds9. In this color scheme magenta and blue are low values while green, yellow, and red are high values. You can see that the images found above correspond to minima, saddle points, and maxima of the time delay surface, as they should.
Finally, you can also work with an extended source. Here is how to specify an extended source, in this case a Sersic profile with n=0.5 (a Gaussian), with a half-light radius of 0.03. We use a circular source for simplicity, but you can give it an arbitrary ellipticity and position angle. The position is chosen to correspond to source #2 above. The last entry is the keyword macro to indicate that we want to use the extended source for macrolensing. (By contrast, a source labeled micro would be treated as a point source for macrolensing and an extended source only for microlensing.)
setsource 1
sersic 1.0 0.1 0.2 0.0 0.0 0.03 0 0.5 macro
0 0 0 0 0 0 0 0
To plot the images, use the commend SBmap2 for "Surface Brightness map". This command takes the full surface brightness distribution in the source plane (including all sources given in the setsource command) and produces the corresponding surface brightness distribution in the image plane.
SBmap2 -2.0 2.0 1000 -2.0 2.0 1000 1 SBmap.fits 3
Here is the result: