# ---
# jupyter:
#   jupytext:
#     text_representation:
#       extension: .py
#       format_name: light
#       format_version: '1.5'
#       jupytext_version: 1.14.4
#   kernelspec:
#     display_name: Python 3
#     language: python
#     name: python3
# ---

# # Surface Brightness Profiles
# Measuring the radius of a galaxy can be tricky.  Where do you define the 'edge'?  One way astronomers often define a galaxy radius is by modelling how the light from the galaxy falls off with radius, and using the properties of this fit to find the radius that encompasses half the light (i.e. the half-light radius).  Today, we are going to compare how the surface brightness of an elliptical and spiral galaxy differ as a function of radius.

# ![Screen%20Shot%202020-07-20%20at%203.15.17%20PM.png](attachment:Screen%20Shot%202020-07-20%20at%203.15.17%20PM.png)

# We are going to do make a couple plots like this one (in a semi-simplified way) using aperature photometry on two nearby galaxies, M89 and NGC0628.  Let's begin by loading the code and data!

from astropy.io import fits
import matplotlib.pyplot as plt
import numpy as np
from astropy import units as u
from astropy.coordinates import SkyCoord
from astropy.wcs.wcs import WCS
from matplotlib import colors


def ap_phot(x,y,rad,data):
    sum = 0.0
    r = int(rad)+1
    i = x - 2*rad -1.
    j = y - 2*rad
    while (i <= x+2*r):
        j = y - 2*rad
        i = i + 1
        while (j <= y+2*r):
            if(np.sqrt((x-i)**2+(y-j)**2) <= r):
                sum = sum + data[int(i)][int(j)]
                j = j + 1.0
            else:
                sum = sum
                j = j + 1.0
    return sum


spiral_file = 'ngc628.fits'
elliptical_file = 'M89_elliptical.fits'
hdu_s = fits.open(spiral_file)
hdu_M89 = fits.open(elliptical_file)
s_data = hdu_s[0].data
M89_data = hdu_M89[0].data

# Now, let's take a look at the two galaxies I've chosen.  How would you classify them?

plt.imshow(s_data, cmap='gray', norm=colors.LogNorm())
plt.colorbar()

plt.imshow(M89_data, cmap='gray', norm=colors.LogNorm())
plt.colorbar()

# Great!  Now we can get to the aperature photometry.  We're going to do a quick-and-dirty version of this for now, so we will simply use circular aperatures centered on each galaxy.  The central RA and Dec for each are:
#
# NGC0628
# $$\rm{RA} = 24.174050, \rm{Dec} = 15.783461 $$
#
# M89
# $$\rm{RA} = 188.9159, \rm{Dec} = 12.5563 $$
#

# I've gone ahead and made one aperature with a five arcsec radius for M89.  We'll walk through M89 together, then you all will do a similar process for NGC0628.

M89_ra = 188.9159
M89_dec = 12.5563
M89_nuc = SkyCoord(M89_ra, M89_dec, unit='deg', frame='icrs')
wcs = WCS(hdu_M89[0].header)
coords = wcs.wcs_world2pix(M89_ra, M89_dec, 0)
e_phot = ap_phot(coords[0], coords[1], 5.0/1.7, M89_data)
print(e_phot)

# This would be the flux for the center of this galaxy. We want to measure surface brightness as a function of radius, so we will need multiple aperatures each growing in size.  I've gone ahead and made a set for M89.  To do this I used a for loop.  For loops repeat a task a specified number of times. So, in this for loop, I have told python to plot circles with radii for each value in our radius array.  I then tell python to extract the flux in that circle.  Everything stays the same for each step except for the size of the radius.  In the last step, I add the new flux value to an array with the command 'append'.  Append takes the existing flux array and sticks our new value onto the end.  

# +
fig, ax = plt.subplots()
coords = wcs.wcs_world2pix(M89_ra, M89_dec, 0)
radius = [5.0, 10.0, 15.0, 20., 30., 40., 50., 60., 70., 80., 90., 100., 150.]
cdelt = 1.7 #arcsec/pixel conversion
rad_pix = [x/cdelt for x in radius]
flux = []
for x in rad_pix:
    circle = plt.Circle((coords[0], coords[1]), x, color='cyan', fill=False)
    ax.add_artist(circle)
    e_phot = ap_phot(coords[0], coords[1], x, M89_data)
    flux = np.append(flux, e_phot)

ax.imshow(e_data, cmap='gray_r', norm=colors.LogNorm(), origin='lower')
print(flux)
# -

# Now we have all the data we need to plot surface brightness as a function of radius! To get a surface brightness, we need to divide by the area (so it's a brightness per area).  We extracted our fluxes over a circular aperature, so the area is $\pi r^2$.  I've also added some axis labels with the matplotlib xlabel and ylabel commands.

area = [np.pi * x**2 for x in radius]
plt.scatter(radius, flux/area)
plt.xlabel('Radius (arcsec)')
plt.ylabel('Surface Brightness (arbitrary units)')
plt.show()

# I decided not to worry about converting the counts to a more reasonable flux unit, since we are mostly interested in the shape of the curve.  What do you notice?
#
#
# Once you are done making observations, repeat the above process for NGC0628!  How is it different than M89?







# Since this is an in-class activity, I skipped a few steps for expediancy.  What do you think these steps are?  What might you do differently if you were designing a project to examine how surface brighness changes with radius in different galaxy types?



StarcoderdataJupyterScriptsDedupFiltered