-- 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

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))

This topic: Students/UnderGraduates > UserList > JeanDuffy > JeanDuffyBScProject > JeanDuffyPythonCode
Topic revision: r3 - 31 Oct 2013 - JeanDuffy

Copyright © 2008-2021 by the contributing authors. All material on this collaboration platform is the property of the contributing authors.
Ideas, requests, problems regarding RHUL Physics Department TWiki? Send feedback