In [1]:
from astropy.io import fits
from astropy import wcs
import wcsaxes
import numpy as np
import aplpy
import pandas
import matplotlib.pyplot as plt
from astropy.modeling import models, fitting
from scipy import optimize
%matplotlib inline
from os.path import expanduser
In [2]:
home = expanduser("~")
gc_dir = home + "/Dropbox/GalacticCenter/"
filename_skymap = gc_dir + "fits/release_galactic_skymap.fits"
HDUlist = fits.open(filename_skymap)
HDUlist.info()
# get excess data
skymap_data, skymap_header = fits.getdata(filename_skymap, header = True)
In [3]:
# world coordinate system
wcs_SgrA = wcs.WCS(filename_skymap)
wcs_GC = wcs.WCS()
w = wcs.WCS(HDUlist[0].header)
#fig = plt.figure(figsize=(8,8))
#ax = fig.add_subplot(1, 1, 1, projection = wcs_SgrA)
#wx, wy = w.wcs_pix2world(0,0)
In [4]:
skymap_data_n = np.nan_to_num(skymap_data)
img = plt.imshow(skymap_data_n, origin='lower')
ax = img.axes
ax.grid(color='white', alpha=1, ls='solid')
ax.set_xlabel("Galactic Longitude")
ax.set_ylabel("Galactic Latitude")
plt.colorbar()
Out[4]:
In [5]:
def TwoD_Gaussian(coord, amp, x0, y0, sigma_x, sigma_y, base):
"""simple model of two-dimensional gaussian"""
x, y = coord
g = amp*np.exp(-((x-x0)**2/(2*sigma_x**2)+(y-y0)**2/(2*sigma_y**2))) + base
return g.ravel()
# TwoD_Gaussian
In [6]:
nx, ny = skymap_data.shape
# Create x and y indices
x = np.linspace(0, nx-1, nx)
y = np.linspace(0, ny-1, ny)
x, y = np.meshgrid(x, y)
coords = x, y
Run fit optimization
In [7]:
guess_SgrA_world = (1000, 0, 0, 0.2, 0.2, 100)
guess_SgrA_pix = (100, 120, 120, 5, 5, 1)
guess_G09_world = (250, 0.08, 0.87, 0.5, 0.5, 50)
#guess_G09_pix = (250,)
#guess_J1745_world = (300)
#guess_J1745_pix = (300)
# Fitting for uncorrelated map, by pixel
p_opt, p_cov = optimize.curve_fit(TwoD_Gaussian, coords, skymap_data_n.ravel(), p0=guess_SgrA_pix, maxfev=2500)
#p_opt, p_cov = optimize.curve_fit(TwoD_Gaussian, coords_world, excess_data_N.ravel(), p0=guess_world_G09,maxfev=2500)
print(p_opt)
print(p_cov)
#print(y[0:3,8:12])
#init_guess = (1000,120,120,1,1,100)
#p_opt_G09, p_cov_G09 = optimize.curve_fit(TwoD_Gaussian, coords, residual_skymap.ravel(), p0=init_guess_G09)
#p_opt, p_cov = scipy.optimize.curve_fit(TwoD_Gaussian, (x,y), dataN.ravel(), p0=init_guess)
In [8]:
pointSource_excess = TwoD_Gaussian(coords, *p_opt)
#print(pointSource_excess.reshape(240, 240)[115:125, 115:125])
residual_skymap = skymap_data_n - pointSource_excess.reshape(nx, ny)
fig_residual = plt.figure()
plt.axis([50, 190, 80, 150])
plt.imshow(residual_skymap)
plt.colorbar()
#print(residual_skymap[115:125, 115:125])
fig_PS = plt.figure()
plt.imshow(pointSource_excess.reshape(nx, ny))
Out[8]:
In [9]:
with open(gc_dir+'/spectralPoints/SgrA_4tels_noPoor_dNdE_TeV.txt') as infile:
filetext = infile.read()
linelist = filetext.split(sep='\n')
for index, line in enumerate(linelist):
paramlist = line.split()
if len(paramlist) < 2:
del linelist[index]
for line in linelist:
print(line)
infile.close()
lines = list(open(gc_dir+"/spectralPoints/SgrA_4tels_noPoor_dNdE_TeV.txt", "r"))
linesSanitized = map(lambda each:each.strip("\n"), lines)
print(linesSanitized)
Celsius = [39.2, 36.5, 37.3, 37.8]
Fahrenheit = map(lambda x: (float(9)/5)*x + 32, Celsius)
print(*Fahrenheit)