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>

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

 
This site is powered by the TWiki collaboration platformCopyright © 2008-2017 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