Functions to open FIT files, make histograms and fit a Gaussian. Code is untested.
from astropy.io import fits
import numpy as _np
import matplotlib.pyplot as _plt
from scipy.optimize import curve_fit
import glob
import ntpath
import math as _m
def OpenImage(FileName):
"""Opens a FIT image and puts pixel ADU values into an array"""
ImageFile = fits.open(FileName)
ImageFile0 = ImageFile[0]
Data = ImageFile0.data
A = Data.flatten() * 2.39
return A
def OpenDirectory(Directory):
filenames = glob.glob(Directory+"*.fit")
return filenames
def Histogram(Array, Bins):
"""Plots histogram from flattened array of data, automatic range for gaussian distribution in array."""
h_test = _plt.hist(Array, Bins, (min(Array), max(Array))) #Histogram #1 over entire array
_plt.clf()
test_y = h_test[0]
test_x = h_test[1]
PeakPosition = _np.argmax(test_y) #Find peak of histogram
Peak = test_x[PeakPosition] #Find position of peak along x axis
BottomofRange = Peak-300 #Tailor new range around peak, very vague estimate, could probably use FWHM instead? Or a percentage of peak
TopofRange = Peak+300
EstRange = (BottomofRange, TopofRange)#Estimate 1 of range
Test2_h = _plt.hist(Array, Bins, EstRange) #Histogram #2 has slightly better range, still not great
_plt.clf()
Test2_y = Test2_h[0]
Test2_x = Test2_h[1]
Counter = 0 #loop over the y values (Counts in histogram) to find the first non-zero value.
for i in Test2_y:
if i == 0:
Counter+=1
if i >0:
break
Test2RangeLower = Test2_x[Counter] #The next test range is between all values that are non-zero,
#this should be changed to a number dependent on the values in the histogram,
#because they might not all be non-zero.
Counter = 0 #loop over reversed list, to cut off any non-zero values.
for i in reversed(Test2_y):
if i == 0:
Counter+=1
if i >0:
break
Test2RangeUpper = Test2_x[len(Test2_y) - Counter]
Test2Range = (Test2RangeLower, Test2RangeUpper)
Test3_h = _plt.hist(Array, Bins, Test2Range)
_plt.clf()
Test3_x = Test3_h[1]
Test3_y = Test3_h[0]
Counter = 0
for i in Test3_y:
if i<2000: #Instead of value here could have a percentage of peak count?
Counter +=1
if i>=2000:
break
NewRangeLower = int(round(Test3_x[Counter]))
Counter = 0
for i in reversed(Test3_y):
if i<2000:
Counter +=1
if i>=2000:
break
NewRangeHigher = int(round(Test3_x[len(Test3_y) - Counter]))
FinalRange = (NewRangeLower, NewRangeHigher)
Difference = NewRangeHigher-NewRangeLower
bins = (Difference/2)-7
h = _plt.hist(Array, bins, FinalRange)
xedge = h[1]
y = h[0]
xcentre = (xedge + _np.roll(xedge, -1))/2.0
x = xcentre[0:len(y)]
return x, y, NewRangeLower, NewRangeHigher, bins
def GaussFunc(x,a,x0,sigma):
return a*_np.exp(-(x-x0)**2/(2*sigma**2))
def FitAGaussian(Array): #Not Working
"""Fits a Gaussian Distribution over histogram"""
h = Histogram(Array, 100)
x = h[0]
y = h[1]
Range = (h[2], h[3])
bins = h[4]
n = len(x)
mean = sum(x)/n
error = []
for i in y:
error.append(_m.sqrt(i))
popt, pcov = curve_fit(GaussFunc, x, y, (1, mean, 1), absolute_sigma=True)
_plt.clf()
_plt.plot(x, GaussFunc(x, *popt), color = 'red')
_plt.hist(Array,bins, Range)
return popt, pcov
def FindExposureTime(path):
head, tail = ntpath.split(path)
B = tail
#B = ntpath.basename(head)
ExposureTime = float(B[:3])
return ExposureTime
def mean_fit(t, mean0, mu):
return mean0 + mu*t
def MakeMeanPlot(Directory):
PictureList = OpenDirectory(Directory)
Mean = []
ExposureTime = []
for i in PictureList:
Exp = FindExposureTime(i)
Array = OpenImage(i)
popt, pcov = FitAGaussian(Array)
Mean.append(popt[1])
ExposureTime.append(Exp)
_plt.clf(),
_plt.scatter(ExposureTime,Mean, label = 'Data', color = 'red')
popt, m_pcov = curve_fit(mean_fit, ExposureTime, Mean, absolute_sigma= 'True')
MF=[]
for i in ExposureTime:
MF.append(mean_fit(i,popt[0], popt[1]))
print MF
_plt.plot(ExposureTime, MF, label = r'Fit $\nu = ' + str(_np.round(popt[0], 2))+ ' + ' + str(_np.round(popt[1], 2)) + ' t $')
_plt.xlabel('Exposure Time')
_plt.ylabel('Mean')
_plt.title('Linear fit of mean vs exposure time of dark frames')
_plt.legend()
_plt.show()
def Sigma_fit(t, mu, si_0):
return _np.sqrt((mu * t) + (si_0**2))
def MakeSigmaPlot(Directory):
PictureList = OpenDirectory(Directory)
sigma = []
ExposureTime = []
for i in PictureList:
Exp = FindExposureTime(i)
Array = OpenImage(i)
popt, pcov = FitAGaussian(Array)
sigma.append(abs(popt[2]))
ExposureTime.append(Exp)
_plt.clf(),
_plt.scatter(ExposureTime,sigma , label = 'Data', color = 'red')
popt, m_pcov = curve_fit(Sigma_fit, ExposureTime, sigma, absolute_sigma= 'True')
SF=[]
for i in ExposureTime:
SF.append(Sigma_fit(i,popt[0], popt[1]))
print SF
_plt.plot(ExposureTime, SF, label = r'Fit $\nu = ' + str(_np.round(popt[0], 2))+ ' + ' + str(_np.round(popt[1], 2)) + ' t $')
_plt.xlabel('Exposure Time')
_plt.ylabel('Sigma')
_plt.title('Linear fit of Sigma vs exposure time of dark frames')
_plt.legend()
_plt.show()
MakeSigmaPlot("C:\Users\Blake\Documents\Msci_project/")
--
BekiChafer - 13 Oct 2017</verbatim>