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
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
test_x = h_test

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
Test2_x = Test2_h

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
Test3_y = Test3_h

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
y = h
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
y = h

Range = (h, h)
bins = h
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):
B = tail
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)
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, popt))
print MF

_plt.plot(ExposureTime, MF, label = r'Fit $\nu = ' + str(_np.round(popt, 2))+ ' + ' + str(_np.round(popt, 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))
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, popt))
print SF

_plt.plot(ExposureTime, SF, label = r'Fit $\nu = ' + str(_np.round(popt, 2))+ ' + ' + str(_np.round(popt, 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>

Physics WebpagesRHUL WebpagesCampus Connect • Royal Holloway, University of London, Egham, Surrey TW20 0EX; Tel/Fax +44 (0)1784 434455/437520

Topic revision: r4 - 19 Oct 2017 - LunaMarieBorella1 Home   Public Web  P P P P P View Edit        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