--
JeanDuffy - 25 Oct 2013
Extraction of useful information from a fit file
import os
import pyfits
import matplotlib.pyplot as plt
import numpy as np
import math
# Plot a fits file
os.chdir("C:/Users/Jean/Desktop/Python/Bsc/Fits") # Changes the directory
f=pyfits.open("Atlas_pleione_1sec.fit") # Open the fit file "Atlas_pleione_1sec.fit" which opens it as an HDU (Header Data Unit) list and set to f.
i=f[0] # Set i to be an HDU object.
d=i.data # Set d to be an array of data.
plt.figure('Atlas & Pleione'),plt.imshow(d), plt.title('Atlas & Pleione') # Plots the array and creates a specific figure for the plot
print pyfits.getheader("Atlas_pleione_1sec.fit") # Gets and shows all the Header information from the fit file
print
print 'Mean of d =',d.mean() # Prints the mean of d
print
print 'Standard Deviation of d =',d.std() # Prints the standard deviation of d
print
print 'Total Flux =',d.sum() # Prints the sum of d, whihc is the total flux
print
print 'Sum of the ADU in the x-axis',d.sum(0) # Prints the sum of values of the pixels in direction 0
print 'Sum of the ADU in the y-axis',d.sum(1) # Prints the sum of values of the pixels in direction 1
print
d2=d[720:760,870:915] # Create a sub area, which is the subarea around Atlas (Brighter sar)
d3=d[220:265,835:885] # Create a sub area around Pleione
# Plot the sub area of Atlas and Pleione on a new figure
plt.figure('Sub Areas'),
plt.subplot(211), plt.imshow(d3), plt.title('Pleione'), plt.text(155,20,d3.sum()), plt.text(110,20,r'Total Flux='),
plt.subplot(212), plt.imshow(d2), plt.title('Atlas'), plt.text(155,20,d2.sum()), plt.text(110,20,r'Total Flux='),
y1=d.sum(1) # Sum the values of all the pixels in direction 1
y1_2=y1[0:300] # Subset of y1 covering the smaller peak
x1=np.arange(0,len(y1),1) # Create an array of length y1, starting from 0 in increments of 1
xmean=(x1[600:]*y1[600:]).sum()/y1[600:].sum() # The mean of the values in a specified region
print 'Weighted mean =',xmean
print
# Plot x1 against y1 and calculate the maximums and their y position on the plot for Figure 3
plt.figure(3),plt.plot(x1,y1), plt.title("ADU against y1"),
plt.xlabel('y1'), plt.ylabel('ADU')
max_y1=max(y1) # Maximum y value
max_y1_2=max(y1_2) # Maximum y value in subset y1_2
max_x1=x1[y1.argmax()] # x value corresponding to max y value
max_x1_2=x1[y1_2.argmax()] # x value corresponding to max y value in subset y1_2
print 'max_x1 =',max_x1 # Print the maximum y value
print 'max_x1_2 =',max_x1_2 # Print the maximum y value in subset y1_2
print 'max_y1 =',max_y1 # Print the x value corresponding to max y value
print 'max_y1_2 =',max_y1_2 # Print the x value corresponding to max y value in subset y1_2
# Plot x0 against y0 and calculate the maximums and their y position on the plot for Figure 4
y0=d.sum(0) # Sum the values of all the pixels in direction 0
y0_2=y0[0:880] # Subset of y0 covering the smaller peak
x0=np.arange(0,len(y0),1) # Array of legnth y0 strarting at 0 in increments of 1
plt.figure(4),plt.plot(x0,y0), plt.title("Other one"),
print
max_y0=max(y0) # Maximum y value
max_y0_2=max(y0_2) # Maximum y value in subset y0_2
max_x0=x0[y0.argmax()] # x value corresponding to max y value
max_x0_2=x0[y0_2.argmax()] # x value corresponding to max y value in subset y0_2
print 'max_x0 =',max_x0 # Print the maximum y value
print 'max_x0_2 =',max_x0_2 # Print the maximum y value in subset y0_2
print 'max_y0 =',max_y0 # Print the x value corresponding to max y value
print 'max_y0_2 =',max_y0_2 # Print the x value corresponding to max y value in subset y0_2
print
# Calculating the angular seperation
d0=max_x0-max_x0_2 # Distance in pixels of the stars in direction 0
d1=max_x1-max_x1_2 # Distance in pixels of the stars in direction 1
d=math.sqrt((d0**2)+(d1**2)) # Total distance between the stars
print 'Distance between the peaks on the CCD in pixels =',d
print
P=0.61 # Plate scale in arcseconds per pixel
theta=P*d # Theta - Angular seperation, P - Plate Scale, d - Seperation on CCD in pixels
print 'Angular seperation in arcseconds =',theta
plt.show()
# Calculate the absolute magnitude where m is the apparent magnitude and dis is the distance to the star in light years
m=3.62
dis=435
print 'Absolute magnitude =',m-5*(math.log10(dis/32.6))