using_pyfits.ipynb

Kernel: Unknown Kernel

Use of PyFits

#Importing libraries
import numpy as np
%pylab inline
import matplotlib.pylab as plt
import scipy.integrate as integ
import scipy.interpolate as interp

#PyFITS library
import pyfits

Populating the interactive namespace from numpy and matplotlib

#Loading data with pyfits
spectras, header = pyfits.getdata("NGC5947.V500.rscube.fits",0, header=True)
#==================================================================
#Parameters
#==================================================================
#Minim wavelenght
LMIN = 3749
#Maxim wavelenght
LMAX = 7501
#Size of spectras
nspec = len(spectras)

#Wavelenghts of all spectras
wavelenght = np.linspace( LMIN, LMAX, nspec )
#Detecting number of pixels
nx = len(spectras[0])
ny = len(spectras[0,0])

#==================================================================
#Importing data from FITS format to ascii data
#==================================================================
for i in xrange(nx):
    for j in xrange(ny):        
        np.savetxt( "spectra_%d_%d.dat"%(i,j), np.transpose([wavelenght, spectras[:,i,j]]), fmt="%d\t%1.4f" )

def Integrated_Galaxy( wlmin=LMIN, wlmax=LMAX ):
    galaxy = np.zeros( (nx,ny) )
    #Calculating index associated to wlmin and wlmax
    i_wlmin = int(float(wlmin-LMIN)/(LMAX-LMIN)*nspec)
    i_wlmax = int(float(wlmax-LMIN)/(LMAX-LMIN)*nspec)

    #Calculating brightness
    for i in xrange(nx):
        for j in xrange(ny):
            #erg s^-1 cm^-2 \AA^-1
            #Simpson method
            galaxy[i,j] = integ.simps( spectras[i_wlmin:i_wlmax,i,j], wavelenght[i_wlmin:i_wlmax] )
            
    #Returning galaxy with integrated bright
    return galaxy

Full integrated galaxy

#Galaxy Picture
galaxy = Integrated_Galaxy(  )

plt.figure( figsize=(8,8) )
plt.imshow(galaxy[::-1], cmap='jet')

<matplotlib.image.AxesImage at 0x7f50c07a02d0>

RGB contributions

#Color Ranges
#Red
red = 5800, LMAX
#Green
green = 4900, 5800
#Blue
blue = LMIN, 4900

plt.figure( figsize=(20,10) )

#Total brightness
galaxy = Integrated_Galaxy(  )

#Red
plt.subplot(1,3,1)
red_galaxy = Integrated_Galaxy( red[0], red[1] )
red_galaxy /= np.max(red_galaxy)
plt.imshow(red_galaxy[::-1], cmap='jet')
plt.title("RED")

#Green
plt.subplot(1,3,2)
green_galaxy = Integrated_Galaxy( green[0], green[1] )
green_galaxy /= np.max(green_galaxy)
plt.imshow(green_galaxy[::-1], cmap='jet')
plt.title("GREEN")

#Blue
plt.subplot(1,3,3)
blue_galaxy = Integrated_Galaxy( blue[0], blue[1] )
blue_galaxy /= np.max(blue_galaxy)
plt.imshow(blue_galaxy[::-1], cmap='jet')
plt.title("BLUE")

<matplotlib.text.Text at 0x7f50c0a9b750>

False colour image

#RGB galaxy
rgb_galaxy = np.zeros( (nx,ny,3) )

for i in xrange(nx):
    for j in xrange(ny):
        
        rgb_galaxy[i,j] = [red_galaxy[i,j],green_galaxy[i,j],blue_galaxy[i,j]]

plt.figure(figsize=(8,8))
plt.imshow(rgb_galaxy[::-1,:,:])

<matplotlib.image.AxesImage at 0x7f50c0c1c150>

Spectrum of the galaxy

spectrum = np.zeros( nspec )
for i in xrange(nspec):
    spectrum[i] = np.mean(spectras[i,:,:])
    

plt.figure( figsize=(12,8) )
plt.plot( wavelenght, spectrum )
plt.xlim( (LMIN,LMAX) )
plt.grid()
plt.xlabel( "Wavelenght [$\AA$]" )
plt.ylabel( "Flux [$1$x$10^{-16}$ erg s$^{-1}$ cm$^{-2}$ $\AA^{-1}$]" )

Halpha = wavelenght[np.argsort(spectrum)[-1]]
print "Halpha line in:", Halpha, "angstroms"
print "Real Halpha line in", 6562.8, "angstroms"
print "Redshift of", abs(Halpha-6562.8)/6562.8

Halpha line in: 6693.0 angstroms
Real Halpha line in 6562.8 angstroms
Redshift of 0.01983909307

H-Alpha emission of the galaxy

#Halpha emission
Halpha = wavelenght[np.argsort(spectrum)[-1]]
Ha_range = Halpha-5, Halpha+5

#Galaxy Picture in Halpha
Ha_galaxy = Integrated_Galaxy( Ha_range[0], Ha_range[1] )

plt.figure( figsize=(8,8) )
plt.imshow(Ha_galaxy[::-1], cmap='gray')

<matplotlib.image.AxesImage at 0x7f50bee8f690>

Dominant emission lines

emission_line = np.zeros( (nx, ny) )
for i in xrange(nx):
    for j in xrange(ny):
        emission_line[i,j] = wavelenght[np.argsort(spectras[:,i,j])[-1]]

#Dominant emission lines
plt.figure( figsize=(8,8) )
plt.imshow(emission_line[::-1], cmap='jet')
plt.colorbar()

<matplotlib.colorbar.Colorbar instance at 0x7f50beb22050>