could later make into an add_substructure() function for GravitationalLens.

Let's just straight up hard code an SIE, and plan on refactoring it later to be more general. Let's use the equations in Evans and Witt 2001, http://arxiv.org/pdf/astro-ph/0108374v1.pdf These only apply to "scale free" ie singular isothermal density profiles, but they are fast to compute and will enable interestingly evil features like boxiness later down the line. You should read the Saas Fee lectures too - if Yashar doesn;t have a copy of this that he can lend you, you can download them. Here's an online library of useful works in gravitational lensing: https://groups.diigo.com/group/gravitational-lensing - you should be able to find the Saas Fee lectures there.

Note that the SIS model is fully analytic, and so is a good test of any SIE code! You should think about this a bit: how will you prove to me that your code is calculating deflection angles accurately?

FYI and for future reference, Barkana 1998 worked out how to compute a singular power law elliptical mass distribution (SPEMD): http://arxiv.org/pdf/astro-ph/9802002v1.pdf

In order to quantify bias in measured subhalo mass from numerical integration rather than analytic.

Because we don't need it any more! Read this: http://git-scm.com/book/en/Git-Branching-Basic-Branching-and-Merging

hey warren,
you're first mock data was great. a very good start.
now let's generate another mock, with an elliptical lens and some orientation angle.
(not to get a perfect einstein ring).
Also, did you figure out the scaling of the lens mass (are you giving it in terms of sigma? like km/s?, if so, then it's fine, you should be using a sigma close t0 200 km/s).

Y

generate two sims, one without phase errors, another with variable phase errors (using a realistic phase screen). Run the subhalo finder on both, and compare subhalo delta_chi2 maps.

Find fits image of galaxy on internet somewhere

Run through EvilLens to get lensed fits image.

Input to CASA and use simobserve to get simulated ALMA data.

In order to avoid kappa(0,0) = infinity.

Learning inheritance! :-)

make 2 simulations (including noise, and one or more subhalos). One with decorrelation, another without.
Run the subhalo finder on both and compare the delta_chi2 maps of subhalos.

Since it performs ~10000-40000 (and possibly more) integrals that take ~0.1-1 seconds each, but do not depend on each other at all.

This would yield substantial improvements in speed, which is critical to create ALMA resolution images using numerically calculated deflection angles.

We may also be able to vectorize the integration, which would give additional speed ups, but perhaps it is best to do one thing at a time.

Deflection angles are reasonably good, so more complex mass distributions can now be read in and used. Maybe next we should make our first predicted images, ready to feed to the ALMA simulator?

Read this: http://astropy.readthedocs.org/en/latest/io/fits/#creating-a-new-image-file
Then write def write(self,map="kappa",to="kappa.fits") so that you can call it with SIElens.write(map="kappa",to="analytic_sie_kappa.fits") or similar. Notice how the use of kwargs allow more literate calls!

Specifically, the pixscale in the loaded file has changed from .1 arcsec to 1.0 arcsec. I am thinking that it it has to do with the way the header of the files that we write are interpreted in the set_pixscale function. Is there maybe an operation we can add to the write_kappa_to(fits) function to change the headers of the fits file so that the set_pixscale function understands them correctly?

python Lens_prototype.py

now runs a quick test of the sigmacrit and distance calculations. Note the units of SigmaCrit! These should be Msun / Mpc^2 !

This way we can just call lens.add_subhalo to add a subhalo. Important that we use subhalo mass/position in the main halo to determine tidal radius. Should add a nice evil touch of realism.

CASA has some fantastic tutorials. Start with this one and go through the tutorial to learn how to simulate ALMA observations from your fits images:
http://casaguides.nrao.edu/index.php?title=Simulation_Guide_for_New_Users_(CASA_4.2)

generate a few mock observations for identical lens/source/noise/etc with only different antenna configurations (i.e., uv coverage). Calculate subhalo delta_chi2 maps for them. Which antenna configuration results in best subhalo constraints? (i.e., higher delta_chi2).

Repeat this for a few different source structures. Presumably for sources with smaller clumps longer baselines should help.

Hi Warren! Phil's here visiting today, and we are chatting about the ALMA data, and potential problems with it.

There are three effects I think we need to worry about, amplitude errors, phase gain errors and decoherence. The first two are antenna-based effects, and the phase errors have previously been dealt with in the Cycle 0 lens fitting. The decoherence is a baseline-based corruption and will suppress the visibility amplitude on the longest baselines. This is effectively reducing the high-spatial frequency structure in the maps, precisely the signature of DM substructure. We should probably simulate this in the data and see what it does.

It seems likely that we can add this effect to uncorrupted simulations of ALMA data, because CASA probably doesn't include these effects.

Can you assign this issue to me, and I will get back to you with a suitable reference for the atmospheric structure function, so that you can simulate it in a realistic manner?

Dan

Yours sims were great. The signal to noise is way too off though.

in CASA's simobserve task there's a variable called "inbright"
the description says: "scale surface brightness of brightest pixel e.g. "1.2Jy/pixel"
So, for example, if our entire image is just zero, with only one non-zero pixel at its center, if we set inbright=0.1 Jy then the whole source will be assumed to have a flux of 0.1 Jy.

Now, if we had two non-zero pixels in the image, where one was half the intensity of the other, by setting inbright=0.1 Jy, we're assuming that the source has a total flux of 0.15 (0.1 for the brightest pixel + 0.05 from the second pixel at half intensity).

For out simulations, right now, we want our entire image to have a flux of, say 0.08 Jy. This is the typical observed flux of bright strong lenses in continuum (continuum = dust emission).

so you have to work out what value to give inbright so that your source has a flux of 0.08 Jy.

Update the saboteur such that the phase screen can be made with custom cell sizes. Important for dealing with small distance observational effects such as decoherence.

Important that cell sizes do not change normalization of phase screen. Should write out the math to show that this is the case.

I started a directory called "evillens" and added an init.py file. This is to enable us to do:

import evillens as evil
lens = evil.GravitationalLens(0.4,1.5)

Can you mv Lens_prototype.py into evillens/lens.py and check that the above commands work from teh top level directory, please?

We can now read_image_to(fitsfile), making an output of our image (or data cube) for input to CASA, but I'm wondering if CASA will want more information in the header (such as reference frequency, bandwidth, RA, DEC) than our code currently produces.

It might be that you specify this in CASA (using inwidth ,indirection, etc.), but I wanted to be sure.

In performing the double integral numerically, a way to avoid losing a significant amount of time (which is otherwise lost to a deep series of for loops) is to just use array arithmetic (this is what I would assume that current codes do). The operation I am trying to perform is:

sum( kappa[!=i,!=j] * (x[i,j]-x[!=i,!=j)/((x[i,j]-x[!=i,!=j])**2+(y[i,j]-y[!=i,!=j])**2))

where x[i,j] and y[i,j] are the x and y coordinates at which we are trying to find the deflection angle, and x[!=i,!=j], y[!=i,!=j], kappa[!=i,!=j] are the arrays of x, y, and kappa values, but excluding the elements x[i,j], y[i,j], kappa[i,j](since a point cannot deflect itself). If we can do this, then each element of the deflection angle can be found in a single line of code, and we can get the full arrays cheaply. Python doesn't seem to have a trivial way to do this though.

Do you know of any good way to perform the above sum?

1. check mcmc chains, look converged?
2. plot best-fit model from mcmc chains.
3. make dirty images for model and data and plot residuals

Create simulations with identical source/noise/uv coverage, etc. Ensure that the data span a realistic frequency range (e.g., observe in band 4 with 2 spws in each of the two sidebands, so approx 12GHz apart at center). Calculate delta_chi2 for various choices of frequency averaging: average all spws together, average two in each SB together, divide spws into sub-windows so that effective uv size of long baseline averaging remains close to 12m (dnu/nu*baseline ~ 12m). Determine effect of this averaging on the data.

Just give a quick look through and find everything bad that you can, then heckle me about it.

If in a notebook, render on screen; if not, write a PNG image.

image grid has pixscale = kappa_pixscale / n where n is the oversampling factor

grids are offset by a suitable diagonal vector.

you've made a lot of progress!
You can now generate lensed images. While you're working on your code to improve different aspects of it, you can start using CASA to simulate ALMA observations of your sky models.
You can download CASA here:
http://casa.nrao.edu/casa_obtaining.shtml
it's pretty easy to install.

Allows for lensing of multiwavelength images.

@yasharhezaveh @drphilmarshall , is it safe to do the bilinear interpolation on each plane of the data cube separately? It seems to work on rgb jpeg images at least (conserves color), but I want to be sure before I commit anything.

Add documentation explaining the plots etc.

Analytic SIE lens axis ratio parameter currently extends along semi-major axis without becoming correspondingly smaller on semi-minor axis. This is leading to a mass underestimate from using analytic models to estimate the mass (i.e. the SIELens.get_mass_inside(r) function), which is starting to become a nuisance.

And keep documentation cells up to date! :-) This is going to be fun.

I see tens of percent error in the central regions... Is this good enough, @yasharhezaveh? I'm guessing no, but I might be wrong. What's the acid test, @wmorning ?

Currently, the centroid and axis ratio are being set automatically in the initialization of the SIElens. It would be better if we could specify them explicitly in the init function (as kwargs, that way they don't need to be specified if one doesn't want to do so). This would allow us to customize our lens more appropriately. However, our last attempt to include this feature failed, possibly due to problems involving the super function not approving of how we were trying to do so.

@drphilmarshall you might already know what to do about this.

Current units of ThetaE are not scaled properly to cancel with SigmaCrit (end up with units of Mpc). This implies that it is missing a length dimension in the denominator.

Most likely reason is that we have R_E instead of ThetaE. I think that we need to add a factor of 1/Dd, and maybe some trig also (so as not to induce any error through the small angle approximation)

@drphilmarshall what are your thoughts?

Then:

1. Issue each task as you work through the plan.
2. Crunch through them, closing each issue when its desired output is visible in notebook! :-)