In [ ]:

    

%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
from astropy.io import fits


    

In [ ]:

    

fname = "F090W_flat_1024_1_webbpsf.fits"
hdulist = fits.open(fname)
#print(hdulist[0].header)
pixel_scale = hdulist[0].header['PIXELSCL']
print pixel_scale


    

In [ ]:

    

data = hdulist[0].data
plt.imshow(np.log10(data))


    

In [ ]:

    

pixel_range = 512.*np.sqrt(2)
angle_range = pixel_range*pixel_scale

npix  = 1000
pixel = 10**( np.log10(pixel_range)*np.arange(npix)/float(npix-1) )
angle = pixel * pixel_scale

idx_max = np.where(data[:,:] == data[:,:].max())
print idx_max[0], idx_max[1]
print data.max()
print data[idx_max[0],idx_max[1]]


image_sum = np.sum(data)
print image_sum


    

In [ ]:

    

nx = hdulist[0].header['NAXIS1']
ny = hdulist[0].header['NAXIS2']

i = np.arange(nx)
j = np.arange(ny)

ix, jy = np.meshgrid(i,j)

print ix[:,1]
print jy[:,1]


    

In [ ]:

    

flux = np.zeros(npix)
for ii in range(npix):
    fi = np.where( np.sqrt( (ix-idx_max[0])**2 + (jy-idx_max[1])**2 )<pixel[ii])
    flux[ii] = np.sum(data[fi])


    

In [ ]:

    

plt.plot(pixel,-2.5*np.log10(flux))
plt.xscale('log')
plt.yscale('log')
plt.ylim([1.0e-2,10])


    

In [ ]:

    

plt.plot(angle,-2.5*np.log10(flux))
plt.xscale('log')
plt.yscale('log')
plt.ylim([1.0e-2,10])
plt.xlim([0.03,10])
plt.xlabel(r"$Aperture\,[\mathrm{arcsec}]$")
plt.xticks(fontname="Times New Roman")
plt.yticks(fontname="Times New Roman")
plt.ylabel(r"$Correction\,[\mathrm{mag}]$")
plt.axes().set_aspect('0.75')
plt.text(3,3,'F090W',fontname="Times New Roman",fontsize=12)
plt.savefig("F090W_aperture_correction.png",bbox_inches="tight")


    

In [ ]: