

In [2]:

    
import os
from astropy.io import fits
from astropy.wcs import WCS
from astropy.modeling import models, fitting
import numpy as np 
from scipy import optimize
from matplotlib import pyplot as plt
%matplotlib inline
#%matplotlib notebook
import aplpy

Open data file



In [3]:

    
home = os.path.expanduser("~")
fitsdir = home + "/Dropbox/GalacticCenter/fits/"

#ROOTfile = "./release_galactic_skymap_fixed_head.fits"
ROOTfile = fitsdir+"SgrA_galactic_uncorrelated_excess_map_fixedhead.fits"
HDUlist = fits.open(ROOTfile)
HDUlist.info()
gal_wcs = WCS(HDUlist[0].header)
print(gal_wcs)
HDUlist.close()









    



Filename: /Users/mbuchove/Dropbox/GalacticCenter/fits/SgrA_galactic_uncorrelated_excess_map_fixedhead.fits
No.    Name         Type      Cards   Dimensions   Format
0    UncorrelatedExcess  PrimaryHDU      16   (240, 240)   float32   
WCS Keywords

Number of WCS axes: 2
CTYPE : 'GLON-TAN'  'GLAT-TAN'  
CRVAL : 359.94420870804902  -0.0460128714594389  
CRPIX : 120.0  120.0  
PC1_1 PC1_2  : 1.0  0.0  
PC2_1 PC2_2  : 0.0  1.0  
CDELT : -0.025000000000000001  0.025000000000000001  
NAXIS    : 240 240



In [4]:

    
fig_apl = aplpy.FITSFigure(ROOTfile)
#fig_apl.show_colorscale( smooth=3, vmin=-5., vmax=5.) #cmap='hot',

#col_bar = fig_apl.show_colorbar()
#fig_apl.recenter(0., 0., width=3., height=2.)

#col_bar



In [4]:

    
excess_data, skymap_header = fits.getdata(ROOTfile, header = True)

excess_data_N = np.nan_to_num(excess_data) # rm nan's for fitting



In [5]:

    
#print(skymap_header)
# G0.9+0.1
#print(excess_data_N[121:123,83:85])
#numpy.where(excess_data_N > 100)
# looks like the x/y axes are flipped, and there is a shift of 1 vs ds9

Create coordinate arrays



In [6]:

    
x_bins = excess_data.shape[0]
y_bins = excess_data.shape[1]

# Create x and y indices
x = np.linspace(0, x_bins-1, x_bins)
y = np.linspace(0, y_bins-1, y_bins)
x, y = np.meshgrid(x, y)
coords = x, y

x_world, y_world = fig_apl.pixel2world(x,y)
coords_world = x_world, y_world

print(fig_apl.world2pixel(359.944,-0.046)) #SgrA
print(fig_apl.world2pixel(0.050,-0.150)) #J1746
print(fig_apl.world2pixel(.87,.08)) #G09









    



(120.00834831926974, 120.00051485836507)
(115.76835104537822, 115.84050006370752)
(82.965071714460734, 125.04094076750306)

Fitting Central Gaussian



In [7]:

    
SgrA_bounds = {'amplitude':(10,75), 'x_mean':(119,121), 'y_mean':(119,121), 'x_stddev':(0.25,3.), 'y_stddev':(0.25,3.), 'theta':(-1.,1.)}
SgrA_gauss = models.Gaussian2D(amplitude=30, x_mean=120, y_mean=120, x_stddev=1, y_stddev=1, theta=0, bounds=SgrA_bounds)
SgrA_bounds_2 = {'amplitude':(0.,25.), 'x_mean':(119,121), 'y_mean':(119,121), 'x_stddev':(3,10), 'y_stddev':(3,10), 'theta':(-1.,1.)}
SgrA_gauss_2 = models.Gaussian2D(amplitude=10, x_mean=120, y_mean=120, x_stddev=3, y_stddev=3, theta=0, bounds=SgrA_bounds_2)
J1746_bounds = {'amplitude':(5,15), 'x_mean':(115,117), 'y_mean':(115,117), 'x_stddev':(3,10), 'y_stddev':(3,10), 'theta':(-1.,1.)}
J1746_gauss = models.Gaussian2D(amplitude=10, x_mean=116, y_mean=116, x_stddev=3, y_stddev=3, theta=0, bounds=J1746_bounds)
sum_gauss = SgrA_gauss + SgrA_gauss_2 + J1746_gauss



In [8]:

    
fit_p = fitting.LevMarLSQFitter()
fit_gauss = fit_p(sum_gauss, x, y, excess_data_N, maxiter=5000)
# try to add in sigma, or do this with curve_fit 

for i in range(0,len(fit_gauss.parameters)):
    paramline = fit_gauss.param_names[i] + ": " + str(fit_gauss.parameters[i])
    print(paramline)









    



amplitude_0: 34.588952332
x_mean_0: 120.15161081
y_mean_0: 119.250931146
x_stddev_0: 1.97346044478
y_stddev_0: 2.07791146282
theta_0: -1.0
amplitude_1: 4.38247771279
x_mean_1: 121.0
y_mean_1: 119.0
x_stddev_1: 10.0
y_stddev_1: 3.0
theta_1: -0.139425938354
amplitude_2: 5.0
x_mean_2: 115.0
y_mean_2: 117.0
x_stddev_2: 4.1051940876
y_stddev_2: 5.41531102346
theta_2: -1.0



In [9]:

    
SgrA_double = SgrA_gauss + SgrA_gauss_2
fit_SgrA = fitting.LevMarLSQFitter()
SgrA_fit = fit_SgrA(SgrA_gauss, x, y, excess_data_N, maxiter=2500)

print(SgrA_fit.parameters)









    



[  36.51477917  119.85366859  119.07277224    3.            2.38902626
    0.16745351]



In [10]:

    
background_model = models.Box2D



In [11]:

    
print(fit_p.fit_info['message'])
#print(type(fit_gauss.x_mean_0))
print(fit_gauss.x_stddev_0.value)

print(fig_apl.pixel2world(fit_gauss.x_mean_0.value,fit_gauss.y_mean_0.value))
print(fig_apl.pixel2world(fit_gauss.x_mean_1.value,fit_gauss.y_mean_1.value))

zero_pixel = fig_apl.world2pixel(0,0)
print(zero_pixel)

print(fig_apl.pixel2world(zero_pixel[0]-fit_gauss.x_stddev_1.value,zero_pixel[1]+fit_gauss.y_stddev_1.value))

#fig_ps = pyplot.figure()
#pyplot.imshow(pointsource_excess)









    



Both actual and predicted relative reductions in the sum of squares
  are at most 0.000000
1.97346044478
(359.94041843558801, -0.064739592003644389)
(359.9192086928137, -0.071012863112965202)
(117.76834689699827, 121.84051525404625)
(0.24999671220973596, 0.07499876182130695)



In [12]:

    
pointsource_filename = fitsdir + "/central_point_excesses.fits"
pointsource_excess = fit_gauss(x,y)

if not os.path.exists(pointsource_filename):
    fits.writeto(pointsource_filename, pointsource_excess, skymap_header, clobber=True)



In [13]:

    
fig_apl = aplpy.FITSFigure(pointsource_filename)
fig_apl.show_colorscale(vmin=-1,vmax=50)
#fig_apl

Excess after subtraction of central sources



In [14]:

    
central_subtraction_filename = fitsdir + "/central_subtraction_excess.fits"
residual_skymap = excess_data - pointsource_excess 

fits.writeto(central_subtraction_filename, residual_skymap, skymap_header, clobber=True)



In [15]:

    
fig_apl = aplpy.FITSFigure(central_subtraction_filename)
fig_apl.show_colorscale( smooth=3, cmap='hot' ) #cmap='hot',
col_bar = fig_apl.show_colorbar()
fig_apl.recenter(0., 0., width=3., height=2.)









    



INFO:astropy:Auto-setting vmin to -7.502e+00
INFO:astropy:Auto-setting vmax to  9.329e+00
/Users/mbuchove/Programs/Anaconda3/anaconda/lib/python3.5/site-packages/aplpy/normalize.py:115: RuntimeWarning: invalid value encountered in less
  negative = result < 0.
/Users/mbuchove/Programs/Anaconda3/anaconda/lib/python3.5/site-packages/matplotlib/collections.py:590: FutureWarning: elementwise comparison failed; returning scalar instead, but in the future will perform elementwise comparison
  if self._edgecolors == str('face'):






    



INFO: Auto-setting vmin to -7.502e+00 [aplpy.core]
INFO: Auto-setting vmax to  9.329e+00 [aplpy.core]

Fitting G0.9+0.1



In [16]:

    
G09_bounds = {'x_mean':(81,85), 'y_mean':(123,127)}
G09_gauss = models.Gaussian2D(amplitude=20,x_mean=83,y_mean=125,x_stddev=2.,y_stddev=2.,bounds=G09_bounds)

fitter_G09 = fitting.LevMarLSQFitter()
fit_G09 = fit_p(G09_gauss, x, y, excess_data_N, maxiter=1000)

print(fit_G09.parameters)









    



[   3.48033553   85.          123.            5.18659608    8.77453814
   -1.74155461]



In [17]:

    
G09_filename = fitsdir + "/G09_excesses.fits"
G09_excess = fit_gauss(x,y)

fig_G09 = plt.figure()
plt.imshow(G09_excess)

#fits.writeto(pointsource_filename, pointsource_excess, skymap_header, clobber=False)









    Out[17]:





<matplotlib.image.AxesImage at 0x10bba2b70>



In [18]:

    
residual_G09 = residual_skymap - G09_excess



In [19]:

    
fig_residual_G09 = plt.figure(figsize=(10,10))
fig_residual_G09.suptitle("Sgr A* Subtracted Excess Map")
ax = fig_residual_G09.add_subplot(2,1,1)
ax.set_xlabel("Galactic Longitude")
ax.set_ylabel("Galactic Latitude")
plt.axis([50,190,80,150])
plt.imshow(residual_G09)
cbar = plt.colorbar()
#cbar.solids.set_edgecolors('face')

ax = fig_residual_G09.add_subplot(2,1,2)
ax.set_xlabel("Galactic Longitude")
ax.set_ylabel("Galactic Latitude")
plt.axis([50,190,80,150])
plt.imshow(residual_G09)
cbar = plt.colorbar()
#cbar.solids.set_edgecolors('face')









    



/Users/mbuchove/Programs/Anaconda3/anaconda/lib/python3.5/site-packages/matplotlib/collections.py:590: FutureWarning: elementwise comparison failed; returning scalar instead, but in the future will perform elementwise comparison
  if self._edgecolors == str('face'):

FITS figure of partially subtracted map



In [20]:

    
outFilename = fitsdir + "/subtracted_skymap_G09.fits"

fits.writeto(outFilename, residual_G09, skymap_header, clobber=True)

Other stuff..



In [21]:

    
fig_apl = aplpy.FITSFigure(outFilename)
fig_apl.show_colorscale(cmap='hot', smooth=5) # vmin=0, vmax=8
col_bar = fig_apl.show_colorbar()
#col_bar









    



INFO:astropy:Auto-setting vmin to -1.048e+01
INFO:astropy:Auto-setting vmax to  9.341e+00
/Users/mbuchove/Programs/Anaconda3/anaconda/lib/python3.5/site-packages/aplpy/normalize.py:115: RuntimeWarning: invalid value encountered in less
  negative = result < 0.
/Users/mbuchove/Programs/Anaconda3/anaconda/lib/python3.5/site-packages/matplotlib/collections.py:590: FutureWarning: elementwise comparison failed; returning scalar instead, but in the future will perform elementwise comparison
  if self._edgecolors == str('face'):






    



INFO: Auto-setting vmin to -1.048e+01 [aplpy.core]
INFO: Auto-setting vmax to  9.341e+00 [aplpy.core]



In [22]:

    
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 [23]:

    
init_guess_G0901 = (250, 0.08, 0.87, 0.5, 0.5, 50)
p_opt_G0901, p_cov_G0901 = optimize.curve_fit(TwoD_Gaussian, coords_world, residual_skymap.ravel(), p0=init_guess_G0901)

print(p_opt_G0901)









    



---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
<ipython-input-23-d4d9fb5315fa> in <module>()
      1 init_guess_G0901 = (250, 0.08, 0.87, 0.5, 0.5, 50)
----> 2 p_opt_G0901, p_cov_G0901 = optimize.curve_fit(TwoD_Gaussian, coords_world, residual_skymap.ravel(), p0=init_guess_G0901)
      3 
      4 print(p_opt_G0901)

/Users/mbuchove/Programs/Anaconda3/anaconda/lib/python3.5/site-packages/scipy/optimize/minpack.py in curve_fit(f, xdata, ydata, p0, sigma, absolute_sigma, check_finite, **kw)
    559     # NaNs can not be handled
    560     if check_finite:
--> 561         ydata = np.asarray_chkfinite(ydata)
    562     else:
    563         ydata = np.asarray(ydata)

/Users/mbuchove/Programs/Anaconda3/anaconda/lib/python3.5/site-packages/numpy/lib/function_base.py in asarray_chkfinite(a, dtype, order)
    666     if a.dtype.char in typecodes['AllFloat'] and not np.isfinite(a).all():
    667         raise ValueError(
--> 668             "array must not contain infs or NaNs")
    669     return a
    670 

ValueError: array must not contain infs or NaNs



In [ ]:

    
#fig_residual = pyplot.figure(figsize=(8,8))
fig_residual = pyplot.figure()
plt.axis([50,190,80,150])
plt.imshow(residual_skymap)
plt.Circle((120,120), 50)

cbar = plt.colorbar()
cbar.solids.set_edgecolors('face')

#fig_residual.set_axes



In [ ]:

    
# print(x)
# print(x_world)
# print(x_world[2,6])
# print(y_world[239,6])

print(fig_apl.pixel2world(100,180))
print(fig_apl.pixel2world(101,180))

#excess_data_N[81:85,123:127]
#xworld = fig_apl.pixel2world()

models.Gaussian2D?

fermi tools maximum likelihood Wilke's theorem