<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 Introduction to Photometry
Dora Föhring, University of Hawaii Institute for Astronomy
Aim
Step1:  0.1 Directory Set up
Step2:  0.2 Display images
Step3:  1. Photometry set up
Select part of the image for ease of display.
Step4:  Define starting values. Fill in values here
Step5:  Aperture photometry set up. Play around with adjusting the aperture radii sizes and see the resulting image under 'Tests'
Step8:  1.1 Centroiding
Step9:  Use centroiding algorithm to find the actual centers of the targe and comparison.
Step10:  Inspect PSF to see whether shift makes sense
Step12:  1.2 Aperture Photometry
Science Aperture
Step14:  Sky annulus
Step15:  Extract values from regions
Create mask
Step16:  Define new regions where the target and comparison are centered.
Step17:  Place mask on region
Step18:  Place mask on sky annulus slice.
Step19:  1.3 Tests
a. Display image with target and comparison centers before and after centroiding
Step20:  b. Disply image with aperture mask and sky annulus
Step21:  2. Photometry
2.1 Calculate SNR
Calculate Signal-to-Noise Ratio. CCD noise = sqrt(signal + background + dark current + read noise). Ignore dark current and read noise here.
Step22:  Sum all flux inside target and comparison apertures and divide by number of pixels to get average count per pixel.
Step23:  2.2 Optimize photometry aperture
Step24:  2.3 Calculate the target's magnitude and uncertainty
Given the comparison is of known magnitude of 19.4
<ASSISTANT_TASK:>
Python Code:
import numpy as np
import matplotlib.pyplot as plt
import astropy.io.fits as fits
## make matplotlib appear in the notebook rather than in a new window
%matplotlib inline
Explanation: Introduction to Photometry
Dora Föhring, University of Hawaii Institute for Astronomy
Aim: Demonstrate photometry on a series of bias and flat field corrected images of a Near Earth Asteroid. 
0. Prerequisites
End of explanation
datadir = ''
objname  = '2016HO3'
Explanation: 0.1 Directory Set up
End of explanation
def plotfits(imno):
    img = fits.open(datadir+objname+'_{0:02d}.fits'.format(numb))[0].data
    f = plt.figure(figsize=(10,12))
    #im = plt.imshow(img, cmap='hot')
    im = plt.imshow(img[480:580, 460:600], cmap='hot')
    plt.clim(1800, 2800)
    plt.colorbar(im, fraction=0.034, pad=0.04)
    plt.savefig("figure{0}.png".format(imno))
    plt.show()
numb = 1 
plotfits(numb)
numb = 2
plotfits(numb)
Explanation: 0.2 Display images
End of explanation
partimg = fits.open(datadir+objname+'_01.fits'.format(numb))[0].data[480:580, 460:600]
Explanation: 1. Photometry set up
Select part of the image for ease of display.
End of explanation
targcen = np.array([##,##])  ## target center
compcen = np.array([##,##])  ## comparison center
Explanation: Define starting values. Fill in values here:
End of explanation
searchr = 6  ## search box size
ap_r    = 2   ## aperture radius
sky_inner = 3
sky_outer = 5
Explanation: Aperture photometry set up. Play around with adjusting the aperture radii sizes and see the resulting image under 'Tests'
End of explanation
def cent_weight(n):
    
    Assigns centroid weights
    
    wghts=np.zeros((n),np.float)
    for i in range(n):
        wghts[i]=float(i-n/2)+0.5
    return wghts
def calc_CoM(psf, weights):
    
    Finds Center of Mass of image
    
    cent=np.zeros((2),np.float)
    ### Write Equations for finding Center of Mass here ###
    return cent
Explanation: 1.1 Centroiding: Center of Mass
Calculate Center of Mass (CoM) defined as: $\bar{x} = \frac{\sum A_i x_i}{\sum A_i }$, $\bar{y} = \frac{\sum A_i  y_i}{\sum A_i }$.
End of explanation
## Cut a box between search limits, centered around targcen
targbox = partimg[targcen[0]-searchr : targcen[0]+searchr, targcen[1]-searchr : targcen[1]+searchr]
weights = cent_weight(len(targbox))
tcenoffset = calc_CoM(targbox, weights)
print(tcenoffset)
tcenter = targcen + tcenoffset
Explanation: Use centroiding algorithm to find the actual centers of the targe and comparison.
End of explanation
plt.plot(sum(targbox))
plt.show()
compbox = partimg[compcen[0]-searchr : compcen[0]+searchr, compcen[1]-searchr : compcen[1]+searchr]
compw = cent_weight(len(compbox))
ccenoffset = calc_CoM(compbox,compw)
ccenter = compcen + ccenoffset
print(tcenter)
compw
Explanation: Inspect PSF to see whether shift makes sense
End of explanation
def circle(npix, r1):
    
    Builds a circle
    
    pup=np.zeros((npix,npix),np.int)
    for i in range(npix):
        for j in range(npix):
            r=np.sqrt((float(i-npix/2)+0.5)**2+(float(j-npix/2)+0.5)**2)
            if r<=r1:
                pup[i,j]=1
    return pup
Explanation: 1.2 Aperture Photometry
Science Aperture
End of explanation
def annulus(npix, r_inner,r_outer=-1.):
    
    Builds an annulus
    
    pup=np.zeros((npix,npix),np.int)
    for i in range(npix):
        for j in range(npix):
            #### Fill in annulus form here ####
            if ((r<=r_outer)&(r>=r_inner)):
                pup[i,j]=1
    return pup
Explanation: Sky annulus
End of explanation
circmask = circle(searchr*2, ap_r)
annmask = annulus(searchr*2, sky_inner, sky_outer)
Explanation: Extract values from regions
Create mask
End of explanation
newtarg = partimg[int(round(tcenter[0]))-searchr : int(round(tcenter[0]))+searchr, int(round(tcenter[1]))-searchr : int(round(tcenter[1]))+searchr]
newcomp = partimg[int(round(ccenter[0]))-searchr : int(round(ccenter[0]))+searchr, int(round(ccenter[1]))-searchr : int(round(ccenter[1]))+searchr]
Explanation: Define new regions where the target and comparison are centered.
End of explanation
targaper = newtarg * circmask
compaper = newcomp * circmask
Explanation: Place mask on region
End of explanation
targann = newtarg * annmask
compann = newcomp * annmask
Explanation: Place mask on sky annulus slice.
End of explanation
im = plt.imshow(partimg, cmap='hot')
plt.clim(1800, 2800)
plt.scatter(targcen[1], targcen[0], c='g', marker='x')
plt.scatter(compcen[1], compcen[0], c='g', marker='x')
plt.scatter(tcenter[1], tcenter[0], c='b', marker='x')
plt.scatter(ccenter[1], ccenter[0], c='b', marker='x')
plt.show()
Explanation: 1.3 Tests
a. Display image with target and comparison centers before and after centroiding
End of explanation
im = plt.imshow(targaper, cmap='hot')
plt.clim(1800, 2800)
plt.show()
im = plt.imshow(targann, cmap='hot')
plt.clim(1800, 2800)
plt.show()
Explanation: b. Disply image with aperture mask and sky annulus
End of explanation
def calcsnr(target, bg):
    signal = target - bg
    noise = np.sqrt(signal + bg)
    snr = signal / noise
    return snr, noise   
Explanation: 2. Photometry
2.1 Calculate SNR
Calculate Signal-to-Noise Ratio. CCD noise = sqrt(signal + background + dark current + read noise). Ignore dark current and read noise here.
End of explanation
targc = np.sum(targaper) / np.sum(circmask)
targbg= np.sum(targann) /  np.sum(annmask)
compc = np.sum(compaper) /  np.sum(circmask)
compbg= np.sum(compann) /  np.sum(annmask)
snr, noise = calcsnr(targc, targbg)
print(snr)
snr, noise = calcsnr(compc, compbg)
print(snr)
Explanation: Sum all flux inside target and comparison apertures and divide by number of pixels to get average count per pixel.
End of explanation
## Write code here that tries a range of photometry apertures and finds the best SNR ##
print(bestaper)
print(snr)
Explanation: 2.2 Optimize photometry aperture
End of explanation
targc = circle(searchr*2, ap_r)*newtarg
targskyc = annulus(searchr*2, sky_inner, sky_outer)*newtarg
compc = circle(searchr*2, ap_r)*newcomp
compskyc = annulus(searchr*2, sky_inner, sky_outer)*newcomp
ratio = np.sum(compc)/np.sum(targc)
### complete here ###
### complete here ###
### complete here ###
refmag = 19.4
### complete here ###
print("Measured Magnitude = {:0.3f} ± {:0.3f}".format(mag, sigmamag))
Explanation: 2.3 Calculate the target's magnitude and uncertainty
Given the comparison is of known magnitude of 19.4
End of explanation
<END_TASK>