from astropy import constants as const
print(const.c)

print(const.c.to('km/s'))

print(const.c.to('pc/yr'))

from astropy import units as u
dist_sun = 8.5 * u.kpc
dist_sun

dist_sun.to(u.m)

(u.count / u.m ** 2 / u.s / u.deg ** 2 * (0.1 * u.deg) ** 2  / u.pixel).to(u.count / u.cm ** 2 / u.s / u.pixel)

from astropy.table import Table
fermi_2FGL = Table.read('Data/gll_psc_v08.fit')

print(fermi_2FGL)

from astropy.coordinates import ICRS, Galactic

# Get crab postion from online database 
crab_position = ICRS.from_name('Crab')

# Convert to galactic coordinates
crab_position.transform_to(Galactic)

# Print as string
print(crab_position.to_string())

from astropy.wcs import wcs
from astropy.io import fits

# Read WCS information from fits header
header = fits.getheader('Data/galactic_center_fermi.fits')

# Set up WCS transformation
wcs_transformation = wcs.WCS(header)

# Transform pixel to world coordinates
wcs_transformation.wcs_pix2world(100, 100, 0)


# Transform world to pixel coordinates
wcs_transformation.wcs_world2pix(359.95, 0.05, 0)

import numpy as np
from astropy.modeling import models, fitting

# Generate fake data
y, x = np.indices((101, 101))
model = models.Gaussian2D(10, 50, 50, 10, 10)
data = model(x, y) + 0.5 * np.random.randn(101, 101)

# Fit model
fitter = fitting.NonLinearLSQFitter()
fitted_model = fitter(model, x, y, data)

# Plot data and model
figsize(16, 5)
subplot(1, 2, 1)
imshow(data, origin='lower')
title('Faked data')

subplot(1, 2, 2)
imshow(data - fitted_model(x, y), origin='lower')
title('Residuals')
show()

from astropy.io import fits

# Read fits file and print info
hdulist = fits.open('Data/galactic_center_fermi.fits')
hdulist.info()

# Get and show image data 
data = hdulist['PRIMARY'].data
figsize(16, 5)
imshow(data, origin='lower')
colorbar()

from astropy.convolution import Tophat2DKernel, convolve
from astropy.io import fits

# Read data
data = fits.getdata('Data/galactic_center_fermi.fits')

# Define tophat kernel
tophat_kernel = Tophat2DKernel(10)

# Filter data
filtered_data_top = convolve(data, tophat_kernel)

# Plot image and kernel
figsize(16, 5)
subplot(1, 2, 1)
imshow(filtered_data_top, origin='lower')
title("Filterd data")

subplot(1, 2, 2)
imshow(tophat_kernel, origin='lower', interpolation='none')
title("Kernel Image")

from astropy.cosmology import WMAP9, Planck13 

# Redshift
z = 4

# Calculate distance using WMAP9
print "WMAP: {0}".format(WMAP9.comoving_distance(z))

# Calculate distance using Planck13
print "Planck: {0}".format(Planck13.comoving_distance(z))

import aplpy
from astropy.io import fits

# Set up FITS figure, draw fermi diffuse model in the first energy bin
ff = aplpy.FITSFigure('Data/gll_iem_v02_P6_V11_DIFFUSE.fit', slices=[0])
ff.show_colorscale(vmin=0)
ff.show_arrows(-10, 30, 10, -30, color='white')
ff.add_label(0, 40, "Galactic Center", size=22, color="white")
ff.show_colorbar()

from astroquery import fermi

# Request fermi crab data 
results = fermi.FermiLAT.query_object('M1', energyrange_MeV='1000, 100000', obsdates='2013-01-01 00:00:00, 2013-01-02 00:00:00')
print(result)

# Open fits file
from astropy.io import fits
crab_event_list = fits.open(results[1])

from gammapy.stats import significance_on_off

# Compute significance
significance_on_off(n_on=18, n_off=97, alpha=0.1, method='lima')

from gammapy.spectrum import crab

for ref in crab.refs.keys():
    emin, emax = 1, 10 # TeV
    int_flux = crab.int_flux(emin, emax, ref)
    print("{0:15s}: {1:.3g}".format(ref, float(int_flux)))

from astropy.convolution import Gaussian2DKernel, Tophat2DKernel, convolve
from astropy.io import fits

# Read data
data = fits.getdata('Data/galactic_center_fermi.fits')

# Define filter kernels
gauss_kernel = Gaussian2DKernel(5)
tophat_kernel = Tophat2DKernel(5)

# Filter data
filtered_data_gauss = convolve(data, gauss_kernel, boundary='extend')
filtered_data_top = convolve(data, tophat_kernel, boundary='extend')

# Plot data
figsize(16, 12)
subplot(1, 2, 1)
imshow(filtered_data_gauss, origin='lower')
title('Gauss smooting')

subplot(1, 2, 2)
imshow(filtered_data_top, origin='lower')
title('Tophat smoothing')
show()

from IPython.display import Math

# The King or beta function 
Math("$K(x, \sigma, \gamma) = \\frac{1}{2 \pi \sigma^2} \\left(1 - \\frac{1}{\gamma}\\right)"
     "\\cdot \\left[1 + \\frac{x^2}{2\gamma\sigma^2}\\right]^{-\gamma}$")

# Fermi PSF can be described by core and tail component
Math("PSF(x) = f_{core}K(x, \sigma_{core}, \gamma_{core}) + (1 - f_{core})K(x, \sigma_{tail}, \gamma_{tail})")

# Lets try to model the Fermi PSF with astropy
from numpy import pi
from astropy.modeling.models import Beta2D
from astropy.convolution import Model2DKernel

# Define core and tail parameters
# Values are taken from: 
# http://fermi.gsfc.nasa.gov/ssc/data/analysis/documentation/Cicerone/Cicerone_LAT_IRFs/IRF_PSF.html
f_core = 0.05
sigma_core, sigma_tail = 0.5399, 1.063
gamma_core, gamma_tail = 2.651, 2.932

A_core = 1 / (2 * pi * sigma_core ** 2) * (1 - 1 / gamma_core)
A_tail = 1 / (2 * pi * sigma_tail ** 2) * (1 - 1 / gamma_tail)
psf_core = Beta2D(A_core, 0, 0, 2 * gamma_core * sigma_core ** 2, gamma_core)
psf_tail = Beta2D(A_tail, 0, 0, 2 * gamma_tail * sigma_tail ** 2, gamma_tail)

# Define PSF kernel
psf_kernel = f_core * Model2DKernel(psf_core, x_size=91) + (1 - f_core) * Model2DKernel(psf_tail, x_size=91)

# Plot PSF image and profile
figsize(16, 5)
subplot(1, 2, 1)
imshow(log(psf_kernel))
colorbar()
title('PSF image')

subplot(1, 2, 2)
plot(psf_kernel.array[psf_kernel.center[0], psf_kernel.center[0]:-1])
loglog()
ylim(1e-6, 1)
xlim(1, 5e1)
title('PSF radial profile')
show()

import numpy as np
from astropy.io import fits
from astropy.wcs import wcs
from astropy.convolution import Gaussian2DKernel, convolve

# Read significance, counts and background map
counts = fits.getdata('Data/galactic_center_fermi.fits')
header = fits.getheader('Data/galactic_center_fermi.fits')

# Set up WCS transformation
wcs_transformation = wcs.WCS(header)

# Smooth data with Gaussian kernel of width 0.1°
kernel = Gaussian2DKernel(header['CDELT2'] / 0.1)
smoothed_counts = convolve(counts, kernel)

# Find peak significance
y_max, x_max = np.unravel_index(np.argmax(smoothed_counts), smoothed_counts.shape)

# Transform pixel to Galactic coordinates
lon_max, lat_max = wcs_transformation.wcs_pix2world(x_max, y_max, 0)
print("GLON: {0}, GLAT: {1}".format(lon_max, lat_max))

# Sum up excess in a radius of 0.3°
y, x = np.indices(smoothed_counts.shape)
r2 = (x - x_max) ** 2 + (y - y_max) ** 2
mask = r2 < (0.3 / header['CDELT2']) ** 2
total_counts = counts[mask].sum()

print("Total counts in 0.3 deg: {0}".format(total_counts))

import numpy as np
from astropy.io import fits
from astropy.table import Table
from astropy.wcs import WCS

# Read Fermi LAT catalog source positions
fermi_catalog = Table.read('Data/gll_psc_v08.fit')
lon = fermi_catalog['GLON']
lat = fermi_catalog['GLAT']

# Read ROSAT map
image, header = fits.getdata('Data/ROSAT.fits', header=True)

# Instantiate WCS object
wcs = WCS(header)

# Find pixel positions of LAT sources. Note we use ``0`` here for the last
# argument, since we want zero based indices (for Numpy), not the FITS
# pixel positions.
x, y = wcs.wcs_world2pix(lon, lat, 0)

# Find the nearest integer pixel
x = np.round(x).astype(int)
y = np.round(y).astype(int)

# Find the ROSAT values (note the reversed index order)
values = image[y, x]

# Print out the values
print(values)