<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 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)
Step1:  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
Step2:  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.
Step3:  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
Step4:  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)
Step5:  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.
Step6:  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.
Step7:  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?
Step8:  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.
Step9:  Problem 3) Sersic Profile
The Sersic profile for galaxies can be written as
Step10:  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.
Step11:  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?
<ASSISTANT_TASK:>
Python Code:
import numpy as np
import matplotlib.pyplot as plt
from astropy.io import fits
%matplotlib notebook
Explanation: 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)
End of explanation
gal_im = fits.getdata( # complete
# complete
# complete
plt.colorbar()
plt.tight_layout()
Explanation: 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?
End of explanation
plt.hist( # complete
plt.yscale("log")
Explanation: 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.
End of explanation
unusual_pix = # complete
gal_median_filt = # complete
# complete
# complete
plt.tight_layout()
Explanation: 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.
End of explanation
def y_ellipse(x, h, k, a, b, theta):
    
    # complete
    # complete
    # complete
    # complete
    # complete
    # complete
    return # complete
Explanation: 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.
End of explanation
x_grid = np.linspace(-3.605, 3.605, 1000)
y_vals = np.vstack((np.empty_like(x_grid), np.empty_like(x_grid)))
for x_num, x in enumerate(x_grid):
    roots = y_ellipse( # complete
    y_vals[:,x_num] = np.sort(roots)
    
plt.plot(np.append(x_grid, x_grid[::-1]), np.append(y_vals[0], y_vals[1][::-1]))
plt.xlim(-4, 4)
plt.ylim(-4, 4)
Explanation: 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.
End of explanation
def mean_radial_flux(im_data, h, k, a, b, theta):
    
    # complete
    # complete
    # complete
    # complete
    # complete
    # complete
    return # complete
Explanation: 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.
End of explanation
a5_flux = # complete
a20_flux = # complete
print("The mean flux at 5 pixels is {:.6f}".format( # complete
print("The mean flux at 20 pixels is {:.6f}".format( # complete
Explanation: 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?
End of explanation
r_grid = # complete
mean_rad_prof = # complete
for # complete
plt.plot( # complete
Explanation: 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.
End of explanation
def squared_error( # complete
    
    # complete
    # complete
    # complete
Explanation: 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.
End of explanation
from scipy.optimize import minimize
res = minimize( # complete
                # complete
                # complete
print("The Sersic index is {:.3f}".format( # complete
Explanation: 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.
End of explanation
# complete
Explanation: 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?
End of explanation
<END_TASK>