This page describes how to use the matlab routine invert to invert "observations" of gravitational lenses. This routine accompanies the article Inverting Gravitational Lenses by Peter R. Newbury and Raymond J. Spiteri.

The Figure below is a 128-by-128 pixel simulation showing two gravitationally lensed arcs of light. The arcs are highly distorted images of a uniformly bright, elliptical disks of light hidden behind the lens. The colour of the arcs is a measure of the magnification, with high magnification at the red end of the spectrum and low magnification at the blue end of the spectrum.

Downloading the software	Inverting the lens	Changing the routine

Downloading the software

1. The main program is the matlab routine invert.m. The latest version is "friendlier" than the previous versions (the changes are described in Inverting the lens below.)

Be sure to save this file as "invert.m" in a directory where matlab can find it.

2. The inversion routine uses the matlab optimization routine fminsearch which evaluates a user-supplied misfit function penalty.m.

Be sure to save this file as "penalty.m".

3. The inversion scheme also requires a set of "observations" to invert, like the observations shown above.

As described below, these observations must first be created with the forward modeling routine lens and a corresponding parameter file, SIS2.

Inverting the lens

0. First you need a set of observations to invert.

These are created by first running the forward-modeling routine lens on a parameter file. The observations above were created from the parameters stored in the file SIS2.

matlab>> lens('SIS2')

Use the 'w -- write image to file' command to save the observations. The observations will be stored in the data file "SIS2.mat" (in general, "parameterfile.mat") using matlab's internal data storage format.

1. Run the inversion routine (be sure to enclose the parameter file in single quotes.)

matlab>> invert('SIS2')

The routine 'invert' reads the observation parameters (the center of the field-of-view, the size of the observation, and the size of the pixels) from the parameter file, SIS2. No, invert does not "cheat" and read the unknown model parameters...

2. The simulation above should appear momentarily on the desktop.

If the simulation does not appear and the matlab console indicates an error has occured, check that the parameter file is valid, following the instructions on altering the parameter file. The position of the figure window has been coded into the routine, and you may wish to adapt it to better fit your screen. Look for the line

set(gcf,'Position',[300 300 800 400]);

which generates a 800-by-400 pixel figure window starting 300 pixels over and 300 pixels up from the bottom left of the screen. The Figure window should have a 2:1 width:height ratio to best display the graphics.

3. The console prompts you for further action:

matlab>> invert('SIS2')
 
  Inverting lens
==================
Initial lens mass is sigma=600 km/s.
Left-click on the observations for the initial
location of the lens centre.
 
  Move the centre of the SIS with the n, e, s, w keys.
  Use +/- keys to increase/decrease lens mass.
  Press q to quit.

Left-click on the observations plot to locate the center of the SIS. If you've played with the forward modeling routine lens, you'll recognize the circular symmetry of these lenses. A good initial guess is some point on the diagonal between the two lensed arcs, closer to the small arc. For example, left-click near the point (2.0,2.0).

4. After a computation, a second plot will pop up which shows the source plane.

The speed of this compuation depends on the platform you are using. On a PIII at 700 Mhz with 128 Mbytes RAM, the computation takes about one second.

Each pixel in the observations is traced back to the source plane (via the lens equation) and marked with a dot. Then the ellipse with the smallest area which encloses all the points is found. This ellipse is drawn in magenta in both Figures. As well, the Einstein Ring of the SIS lens is drawn in red over the observations.

5. The matlab console gives you instructions and information about the current "solution":

 
Current misfit = ###.######


6. With the pointer/cursor over the Figure window, use single keystrokes to proceed.

n,e,s,w   Move the center of the SIS north, east, south, or west (as a reminder, these directions appear in the observation window) one pixel in the desired direction.

Note: previous versions of invert use the up, right, down, left arrow keys on the number pad. This was changed when we tried, painfully, to run invert on a laptop with an embedded numpad.

+ / -   Increase or decrease the mass parameter of the SIS, in increments of 50 km/s. (This can be changed in the code for finer adjustment.)

After each keystroke, a computation is performed and the two Figures will be updated and the new value for the misfit will be displayed in the matlab console.

7. Adjust the model parameters.

Continue altering the model parameters to decrease the area of the surrounding ellipse, as indicated by the value of the misfit.

For the example described in the article, first press the '+' to increase the mass from 600 km/s up to 750 km/s (3 keystrokes). You'll notice a big drop in the misfit. Then move the mass center towards (2.5,1.4), which will "fine-tune" the solution.

8. Quit

When all possible keystrokes increase the misfit, you are within one increment the minimum of the misfit. The keystroke 'q' quits the interactive loop, and the results are displayed on the matlab console. The results in the article are described by:

...
Current misfit = 3.722799
Solution:
Information about deflector:
 
Lens centre ( 2.563, 1.373)
SIS mass, sigma=750.000
 
Information about the source:
 
Source centre ( 0.997, 0.301)
Source semi-major axis a= 1.978
Source semi-minor axis b= 0.599
Source ellipticity e=1-b/a=0.697
Source orientation phi=119.7


Your results may be slightly different than these because of differences in the initial guess in Step 3.

Changing the inversion routine

• The size of the increments, particular the mass increment of 50 km/s, can be changed. Either hard-code the changes in the invert.m routine, or add more actions, like 'h' half the increment size and 'H' double the increment size.

• As we mention in Section 5 of the article, an ambitious upgrade is to get matlab to find the optimal set of model parameters. Essentially, you can use this whole invert routine (without the graphics) as the penalty function which is called by fminsearch from a higher loop. After find the BEST best-ellipse, you can use the graphics again.

• A more ambitious upgrade is to change the implementation of the optimization routine. We use a modification of Courant's classical penalty method, but as we mention in Section 5 of the article, there are other possibilities for finding the optimal set of model parameters.