This is the second component of the galaxy modeling unit, which involves simulating a spectroscopic observation (the intensity of light as a function of position and wavelength).
Remember that the overall goal here is to get comfortable hacking at code. How many cells you actually work through is not necessarily a great example of this -- there's a lot of subtleties in code, so sometimes working slowly but carefully can teach you a lot more than trying to get through all the prepared material.
In [ ]:
# only necessary if you're running Python 2.7 or lower
from __future__ import print_function
from __builtin__ import range
In [ ]:
import numpy as np
In [ ]:
# import plotting utility and define our naming alias
from matplotlib import pyplot as plt
# plot figures within the notebook rather than externally
%matplotlib inline
The model for the spatial intensity of the galaxy we observe (i.e. the distribution of brightness on the sky) has two basic components:
Let's initialize each of these in turn using Python's built-in ability to store/define functions.
Our galaxy intensity will be a radially-symmetric exponential profile defined as
$$ I_{\textrm{gal}} = I_0 \exp \left(- a \frac{r}{r_e} \right) $$where $I_0$ is a normalization factor (which we will take to be 1 here), $a = 1.68$ is a constant that governs how quickly light "falls off" towards the outskirts, and $r_e$ is the effective radius of the galaxy.
In [ ]:
# Galaxy intensity model: Exponential
def prof_expo(r, re):
a = 1.68
return np.exp(-a * r / re)
In addition to the "standard" way of defining functions shown above, Python has an additional method that can be used to define simple functions using lambda. These are convenient if you ever want to define functions "on the fly" for any reason, as discussed here. The same function defined in lambda notation is illustrated below.
In [ ]:
# defining galaxy intensity using exponential
prof_expo2 = lambda r, re: np.exp(-1.68 * r / re)
Check that these give the same result using a grid of radii. Feel free to use np.linspace, np.arange, or any other tool to generate the set of radii for testing.
In [ ]:
re = 1. # effective radius
radius = ... # radii
gal1 = ... # function 1 (def)
gal2 = ... # function 2 (lambda)
# numerical checks
# plot results
Experiment with different values for $a$ and $r_e$ and different gridding for radius to see how they change the resulting distribution. Also feel free to play around with the above plot until you have something you're happy with. (You can never spend too much time making good plots!)
Let's turn this radial model into a 2-dimensional image in $\mathbf{x}$ and $\mathbf{y}$. First, let's create a grid of points in x and y to plot. We'll set our galaxy to have an effective radius of 0.5 arcseconds. Check to make sure you've set $a$ back to the default value $a=1.68$ again before you continue.
In [ ]:
# galaxy effective radius
re = 0.5 # in arcsec
# define 2-D grid
Nx, Ny = 1000 + 1, 1050 + 1 # number of grid points in x and y (1 padded for the edge)
x_grid = np.linspace(-5. * re, 5. * re, Nx) # grid in x direction
y_grid = np.linspace(-5. * re, 5. * re, Ny) # grid in y direction
The next step is to turn our two 1-D arrays into a 2-D grid in $\mathbf{x}$ and $\mathbf{y}$. More explicitly, we want to compute the intensity of our galaxy at (x_grid[i], y_grid[j]) for all i and j elements in x_grid and y_grid, respectively. This means we need to create a new array of values x and y of length $N_x \times N_y$.
Spend some time thinking about how you could create these arrays, which contain all possible combinations of x_grid and y_grid, so that if you plotted them using plt.plot you would get a grid. Feel free to try your hand at coding something up below before moving onto the "Pythonic" way of solving this problem.
In [ ]:
# space for experimenting with computing a 2-D grid from 2 1-D grids
As with many problems that come up often but are tedious to implement by hand, Python has a function to do this! The solution is shown below.
In [ ]:
# mesh (x_grid, y_grid) into a new set of 2-D (x, y) arrays
x, y = np.meshgrid(x_grid, y_grid) # x,y for our 2-D grid
print(x, x.shape)
print(y, y.shape)
We see that our new x and y are the same shape, and contain $N_x \times N_y$ elements. Unfortunately, the number of points on our axes is flipped: we wanted 1001 points in the $x$ direction and 1051 in the $y$ direction. This is because the default option in np.meshgrid uses 'xy' (Cartesian) indexing instead of 'ij' (matrix) indexing, which ends up flipping the order of the axes. Since we're using matrices, we actually need to specify the latter.
This result shows how important it is to understand exactly what the "default" options for a particular package are. If you're not exactly sure how something is looking, it's always in your best interest to sanity-check the outputs!
In [ ]:
# *properly* mesh (x_grid, y_grid) into a new set of 2-D (x, y) arrays
x, y = np.meshgrid(x_grid, y_grid, indexing='ij') # x,y for our 2-D grid
# print array and array shapes
print(x, x.shape)
print(y, y.shape)
Let's now visually check the results.
In [ ]:
plt.plot(x, y, '.');
Woah -- what is going on here (assuming your notebook didn't crash making the plot)? Why do we havea bunch of different colored lines?
There's actually a lot to unpack here, so let's take our time getting into exactly what plt is doing here.
First, when plt sees a 2-D array, rather than "flattening" the results and treating it like a 1-D array (i.e. just a string of numbers), it instead plots each column of the array separately. In other words, plt.plot(x, y) actually calls plt.plot(x[:, 0], y[:, 0]), plt.plot(x[:, 1], y[:, 1]), etc., plotting your data $N_y$ times. Since we're using the default setting, each new plot is also automatically assigned a new color based on the default color scheme in matplotlib, which is what leads to the color vomit shown on screen.
Let's check this out by plotting using just a few columns.
In [ ]:
# select 10 columns of the array
x_temp, y_temp = x[:, 15:20], y[:, 15:20] # example of array slicing
plt.figure()
plt.plot(x_temp, y_temp, '.');
# print array and array shape
print(x_temp, x_temp.shape)
print(y_temp, y_temp.shape)
The second thing you might notice is that although we're using points specified by '.' in plt, what we get looks like a thick line. This is also plt working as intended: our points are just so dense that they overlap with each other in that particular region. We can better see that actual results by "thinning" our arrays. This works using the format x[start:stop:step].
In [ ]:
# select 10 columns of the array
x_temp, y_temp = x[::20, 15:20:1], y[::20, 15:20:1] # example of array slicing/thinning
plt.figure()
plt.plot(x_temp, y_temp, '.');
# print array shape
print('x:', x_temp.shape)
print('y:', y_temp.shape)
So with that in hand, let's now see what our grid looks like (thinned by 20). To help out with plotting, we're also going to "flatten" our arrays from 2-D to 1-D before plotting them.
In [ ]:
# thin grid by a factor of 20
x_temp, y_temp = x[::20, ::20].flatten(), y[::20, ::20].flatten() # slicing/thinning/flattening
plt.figure()
plt.plot(x_temp, y_temp, '.', markersize=2)
# print array shape (2-D vs flattened)
print(x[::20, ::20].shape, x_temp.shape)
With our 2-D grid of points in (x, y), we are now ready to compute our 2-D galaxy image. First, we can compute the radius using
after which we can plug our compute radii $\mathbf{r}$ into our function prof_expo for the exponential profile.
In [ ]:
r = ... # 2-D grid of radii
model_gal = prof_expo(r, re) # 2-D grid of galaxy intensity
Rather than trying to struggle with getting images to look good in plt.plot (or, alternately, plt.scatter), we're instead going to use plt.imshow (which is designed for this). Some examples are shown here.
In [ ]:
# plotting our galaxy profile
plt.figure()
# default plot
plt.imshow(model_gal)
# more detailed plot
#plt.imshow(model_gal.T, # take the transpose to flip x and y in plot
# origin='lower', # specify the origin to be at the bottom not the top
# extent=[x_grid[0], x_grid[-1], y_grid[0], y_grid[-1]], # specify [left, right, bottom, top] positions
# cmap='magma', interpolation='none') # additional options
#plt.xlabel('x [arcsec]')
#plt.ylabel('y [arcsec]')
#plt.title('Intrinsic Galaxy Profile')
#plt.colorbar(label='Intensity') # add a colorbar
Notice that there's some weird stuff going on with the default imshow plot:
x has 1051 elements instead of y), andThese should be fixed in the detailed version, which specifies much of these directly. Switch the plot to instead use the detailed version. Read through the documentation and take a look at the different options used in the more detailed version to get a sense of what's changed. Feel free to play around with different options until you have something you like.
Extra Challenge: Try to rotate the label on the colorbar 180 degrees so it now reads vertically in the other direction.
Extra Extra Challenge: See if you can get the color scheme to function logarithmically rather than linearly using plt.imshow's norm argument.
Although the image above looks great, we rarely ever see it in practice from the ground because of the atmospheric point-spread function (PSF). If you imagine our image above as an infinite (rather than finite) collection of points, what the PSF does is turn every one of those points into a small blob. Our final observed image is then, in math terms, a convolution of our galaxy model and our PSF model.
PSFs can be very complicated, but in most cases can be approximated by a Moffat profile quite well, which is like a Normal (i.e. Gaussian, "bell curve") distribution but with heavier ("fatter") tails. This has the form:
$$ I(r \,|\, \alpha, \beta) = \frac{\beta - 1}{\pi \alpha^2} \left[ 1 + \left(\frac{r}{\alpha}\right)^2 \right]^{-\beta} $$where $I(r \,|\, \alpha, \beta)$ stands for the intensity as a function of radius given $\alpha$ and $\beta$, where $\alpha$ and $\beta$ are constants that describe the overall shape of the profile. In general, $\alpha$ describes the "size" while $\beta$ describes how quickly the PSF "falls off" towards the edges (similar to the effective radius $r_e$ we defined above).
When astronomers talk about typical "seeing" conditions, they often are describing the full width at half maximum (FWHM), which measures the width ("diameter") of the distribution at half its maximum value. This relates to $\alpha$ as
$$ \textrm{FWHM}(\alpha, \beta) = 2 \times \alpha \times \sqrt{2^{1/\beta} - 1} \quad . $$Let's now define our PSF profile in terms of the FWHM, fixing $\beta = 4.765$.
In [ ]:
# PSF Model: Moffat
def prof_moffat(r, fwhm):
# compute constants
beta = 4.765 # beta is fixed
bnorm = np.sqrt(2.**(1. / beta) - 1) # constant, computed from beta
alpha = fwhm / 2. / bnorm # alpha, computed from beta and FWHM
# compute PSF
norm = (beta - 1.) / (np.pi * alpha**2)
psf = norm * (1 + (r / alpha)**2)**-beta
return psf
Let's assume a typical seeing of 0.8 arcsec. Compute and plot the PSF profile below. Take the code we used to compute and plot our galaxy profile earlier as a guide. Again, feel free to play around with the parameters to see how they change the PSF and plot; just make sure that $\beta$ is set back its default value of $\beta=4.765$ before moving on.
In [ ]:
# define our typical seeing
psf_fwhm = 0.8 # FWHM [arcsec]
# compute our psf
model_psf = ...
# plot our psf
...
As mentioned above, our final observed image is a convolution of our galaxy model and our PSF model. Convolutions can be tricky to deal with. Naively, we could imagine doing a convolution as follows:
r(x, y) grid into a PSF, making sure we keep the same overall amplitude $I(r(x,y))$ for each PSF.This type of thinking represents the type of thinking that is useful when coding: breaking down each part of a larger problem into a bite-size chunk that is easier to code up before combining things together at the end. (Figuring out where loops are needed versus where array operations can be done instead is also incredibly useful.)
In addition to this discrete, "computing"-oriented way of thinking about the problem, convolutions also lend themselves naturally to a more theoretical/math-y approach (i.e. we're smearing out an infinite collection of points from one continuous function using another continuous function). From that perspective, it turns out there are some really neat tricks we can do to actually compute more "exact" convolutions even with a small set of data points.
As with most applications that are relatively common, Python has a bunch of functions/packages useful for convolutions! We'll use a specific version fftconvolve (which uses fast Fourier transforms) from scipy.signal (a set of functions in scipy useful for signal processing) to do our convolution.
Import fftconvolve from scipy.signal, use it to compute the observed galaxy model, and plot the results.
Note: make sure to compute the convolution using mode='same', otherwise you won't be able to plot your results! Check the documentation for additional details.
In [ ]:
# compute convolution of galaxy and PSF
model_obs = fftconvolve()
model_obs /= np.max(model_obs) # normalize result to 1.
# plot our result
...
Compare our "convolved" image with the original galaxy model and the original PSF. This "smeared out" image now represents the galaxy we'd actually observe through the atmosphere from our ground-based telescope!
Let's imagine we're interested in taking a spectrum of this galaxy. Typically, this works by putting a slitmask on the telescope. These are typically made of metal, and are designed to blocks out everything except for the small amount of light passing through slits that have been drilled in the mask. The reason this has to be done is that the light from each slit is going to be spread out on the physical detector as a function of wavelength, so that blue light (i.e. shorter wavelengths) on one end and red light (i.e. longer wavelengths) are on the other. A general picture of this is shown below.
Using the slit dimensions specified below, see if you can compute (and plot) the galaxy image that would be seen by our telescope through the slitmask. I've included a solution below, but see if you can come up with your own version before moving on.
Extra Challenge: Try and compute the "slit loss" (the fraction of light that is not blocked by the slit). This tells us the overall "efficiency" of our observation so far.
In [ ]:
# slit dimensions [arcsec]
slit_width = 0.8
slit_height = 7.
# post-slit galaxy model
model_spec =
# plotting the results
...
A solution is given below.
In [ ]:
# one possible solution using boolean algebra (0=False, 1=True)
# slit dimensions [arcsec]
slit_width = 0.8
slit_height = 7.
# slit model
model_slit = (abs(x) <= slit_width / 2.) # set all x's within slit_width / 2. of 0. to 1; otherwise set to 0
model_slit *= (abs(y) <= slit_height / 2.) # *also* set all y's within slit_height / 2. to 1; otherwise set to 0
# post-slit observed galaxy model
model_spec = model_slit * model_obs
# plot results (basic)
plt.figure(figsize=(14, 6)) # create figure object
plt.subplot(1, 2, 1) # split the figure into a grid with 1 row and 2 columns; pick subplot 1
plt.imshow(model_slit.T) # plot slitmask model
plt.colorbar(label='Transmission')
plt.subplot(1, 2, 2) # pick subplot 2
plt.imshow(model_spec.T) # plot combined galaxy+slit model
plt.colorbar(label='Intensity')
plt.tight_layout()
We're almost done simulating our galaxy. The last step is to account for the pixel scale of our observation. Telescopes don't have infinite resolving power, but have a minimum resolution determined by the size of the pixels on the charge-coupled device (CCD). This usually is expressed in units of arcseconds per pixel. Since our simulated observation above likely has a much higher resolution than the typical resolving scale, we need to bin it to lower resolution.
numpy has a few functions to accomplish this. The most relevant for our purposes are histogram as histogram2d, which compute the histogram of input data. There are also some nifty shortcuts for plotting these using the plt.hist and plt.hist2d functions.
As a warm up, let's first bin a set of normally distributed random numbers just to get familiar with the function.
In [ ]:
rand = np.random.normal(loc=0., scale=1., size=100000) # normally distributed random numbers
plt.hist(rand); # plot a histogram
Take a look at some of the additional arguments that can be passed to plt.hist. See if you can:
Extra Challenge: Ensure that each histogram is properly normalized (i.e. integrates to 1.), cumulative (each bin is the sum of all previous bins), has 100 evenly spaced bins from -5 to 5, and ignores any values outside of those boundaries.
Extra Extra Challenge: Plot 10 realizations of random numbers such that they "stack" on top of each other, each with a different color based on a pre-defined colormap (see the previous notebook for an example of how to define a colormap).
Now let's look at weighted histograms. Our original exponential profile has a series of radii (stored in radius) corrsponding to a set of intensities (gal1 and gal2). We want to bin this to lower resolution using plt.hist. To do this, we just bin up the radius to lower resolution, where we weight each point by the intensity.
In [ ]:
plt.plot(radius, gal1) # plot original relation
plt.hist(radius, weights=gal1, normed=True); # plot a histogram, normalized based on our weights
Let's bin this using np.histogram to see what a histogram output looks like.
In [ ]:
counts, bin_edges = np.histogram(radius, weights=gal1, normed=True)
print(counts)
print(bin_edges)
With all that done, let's now define the bins set by the pixel scale resolution of the Multiple Mirror Telescope (MMT) and Magellan Infrared Spectrograph (MMIRS) instrument. This is a spectrograph operated by a joint venture between the Smithsonian and the University of Arizona that you could one day use if you end up going to either institution!
Note that the code below uses a number of methods that you might not be familiar with. Take some time to see if you understand everything that's going on well enough that you could explain what the code below is doing to a friend.
In [ ]:
# the pixel scale of the MMIRS instrument
pix_scale = 0.2012 # [arcsec/pix]
# define our bins
x_bin = np.arange(x_grid[0], x_grid[-1] + pix_scale, pix_scale) # bin edges in x with width=pix_scale
y_bin = np.arange(y_grid[0], y_grid[-1] + pix_scale, pix_scale) # bin edges in y with width=pix_scale
# define the centers of each bin
x_cent, y_cent = 0.5 * (x_bin[1:] + x_bin[:-1]), 0.5 * (y_bin[1:] + y_bin[:-1])
# the total number of bins
Nx_pix, Ny_pix = len(x_cent), len(y_cent)
print(Nx_pix, Ny_pix)
Using the bins/bin centers computed above:
model_spec observation to the MMIRS pixel scale using np.histogram2d.plt.imshow.Note that the code below should not work as is: not only are the variables placeholders (i.e. they're undefined), histogram2d also requires that all the input data is 1-D, which means we have to flatten the input arrays (if needed).
In [ ]:
# convert from arcsec to pixels
model_pix, x_edges, y_edges = np.histogram2d(Xs, Ys, weights=Model, bins=[x_bin, y_bin]) # bin over pixel scale
In [ ]:
# plot results
plt.imshow(model_pix.T)
Let's take a step back to survey what we've accomplished:
This together gives us a reasonable idea of the relative distribution of photons we'd expect to see on our detector if we observed our galaxy at a specific wavelength. Cool, huh? :)