Galaxy photometry

Version 0.1

Today our goal is to compute the 1D radial profile of a simulated galaxy. While many software packages exist to do this, we will build this same functionality from scratch. Following this, we will be able to fit a radial profile to a galaxy image and measure its size. In practice, this problem is more complicated, but this example is, nevertheless, instructive.

---

By M Alpaslan (NYU) & AA Miller (CIERA/Northwestern & Adler)

Problem 1) Visualize the Galaxy Image Data

The simulated galaxy for today's problem comes from the Illustris Simulation, which is a large cosmological simulation of galaxy formation.

Download the fits image for our simulated galaxy: https://northwestern.box.com/s/a8qp2afv6ed9r36v0y0vyc3vf82s1bw0

(Note - The Illustris team has made simulated images of many different galaxies available.)

Problem 1a

Using your favorite tool, display an image of the galaxy using a logarithmic stretch. Do you notice anything unusual about the image?

Problem 1b

Plot a histogram of the $\log$ of the pixel values in the simulated image.

Hint - you will likely want more than the matplotlib default number of bins.

Problem 1c

Informed by your histogram, replace any highly discrepant pixels within the image with the median pixel value for the image.

Hint - plot the results to confirm this operation occurred as expected.

Note - this procedure is strongly discouraged with real data: if weird things are happening in an image there is a reason, and it is best to understand the origins of the peculiarities. For our current purposes, this step will simplify our analysis.

Problem 2) Measuring Radial Intensities

Before we can measure the radial profile of this galaxy, we need to compute the flux intensity at a given radius. By convention, this is done by computing the mean flux along an ellipse centered on the galaxy (where the "radius" in this case is the semi-major axis).

An ellipse with arbitrary centroid and rotation angle can be described as

$$ \frac{((x-h)\cos{\theta} - (y-k)\sin{\theta})^2}{a^2} + \frac{((x-h)\sin{\theta} + (y-k)\cos{\theta})^2}{b^2} = 1,$$

where $h, k, a, b, \theta$ are the x centroid, y centroid, semi-major axis, semi-minor axes, and the rotation angle of the elipse, respectively.

This can be re-written as (it's worth checking this on the board or scratch paper):

$$\left(\frac{\sin^2{\theta}}{a^2} + \frac{\cos^2{\theta}}{b^2}\right)y^2 + \left(-2\frac{((x-h)\cos{\theta} + k\sin{\theta})\sin{\theta}}{a^2} + 2\frac{((x-h)\sin{\theta} - k\cos{\theta})\cos{\theta}}{b^2}\right)y + \left(\frac{((x-h)\cos{\theta} + k\sin{\theta})^2}{a^2} + \frac{((x-h)\sin{\theta} - k\cos{\theta})^2}{b^2} - 1\right) = 0$$

Problem 2a

Write a function to solve for $y$ given $x$ and the ellipse parameters.

Hint - the function np.roots may be helpful.

Problem 2b

Plot an elipse using your equation with $h = 0$, $k = 0$, $a = 4$, $b = 2$, and $\theta = \pi/6$.

Hint - most of this code has been provided for you.

Problem 2c

Write a function to measure the average flux of every pixel in the image that is intersected by an ellipse with parameters $k, h, a, b, \theta$.

Hint - np.floor() will be useful for selecting the correct pixel.

Hint 2 - brute force is not elegant, but it is sufficient as long as you throw away imaginary solutions.

For this and the remaining problems, assume that the galaxy is centered at (x, y) = (128, 138), that the ellipticity (1 - b/a) is fixed at 0.3, and that $\theta = 0.8377$. [In principle we could ask you to determine each of these values, but that is outside the scope of the current problem.]

Problem 2d

Measure the mean "radial" flux of the galaxy at $a = 5$ and $a = 20$. Do these answers make sense?

Problem 2e

Plot the "mean radial profile" of the galaxy.

Hint - when choosing a grid of radii to measure, recall that the function does not make sub-pixel estimates of the flux.

Problem 3) Sersic Profile

The Sersic profile for galaxies can be written as:

$$\ln I(R) = \ln I_0 - kR^{1/n},$$

where $I$ is the intensity (or flux), $I_0$ is the intensity at $R = 0$, $k$ is a normalization constant, and $n$ is the "Sersic index", which describes the curvature of the light profile.

We will now use a non-linear optimization routine from scipy.optimize to determine the value of the Sersic index.

Problem 3a

Create a function squared_error() that calculates the sum of the squared difference between $\ln I(R)$ and $\ln I_0 - kR^{1/n}$.

Hint - the 3 unknown parameters ($I_0$, $k$, and $n$) should be passed to the function as a single tuple called params.

Problem 3b

Using the L-BFGS-B method of scipy.optimize.minimize(), determine the value of the Sersic index.

Hint - you will want to provide a reasonable first guess along with bounds for the minimization function.

Challenge Problem

The image that we have been working with reflects the full sensitivity of the Illustris simulation. Ground based images are observed through the atmosphere. Convolve the image with a PSF (assume a circular gaussian with FWHM = 12 pixels) and re-measure the Sersic index. How does your answer compare?