from astropy.io import fits import matplotlib.pyplot as plt import numpy as np from scipy.optimize import curve_fit from scipy.special import gamma import os #============================================================================================================================ # ============================================================================= # plt.rc('text', usetex=True) # plt.rc('font', family='serif') # plt.rcParams.update({'font.size':14}) # ============================================================================= #============================================================================================================================ def OSCheck(): print(os.getcwd()) def OpenFits(file): """Open fit/fits file. Example OpenFits('Atlas.fits') d is an array of the size of the CCD eg. 1020x1530""" f = fits.open(file) d = f[0].data vminVal = 0.95*np.median(d) vmaxVal = 1.1*np.median(d) cmapVal = 'Greys' plt.imshow(d, vmin=vminVal, vmax=vmaxVal, cmap=cmapVal, origin='lowerleft',) plt.show() return d def PSF(file,x,y): """PSF (Point Spread Function) requires an array, d, and line, l an int, along the x-axis""" f = fits.open(file) d = f[0].data plt.plot(d[y,x-100:x+100]) PSFTitle = "PSF Of Position (" + str(x)+','+str(y) + ') in file: ' + str(file) plt.title(PSFTitle) plt.show() def Histogram(file,x,y): """"x,y point about which the histogram will be made""" f = fits.open(file) d = f[0].data array = [] for i in range(x-50,x+50): for j in range(y-50,y+50): array.append(d[i,j]) print('array length', len(array)) HISTTitle = "Histogram of area (" + str(x) +','+str(y)+') in file: ' + str(file) plt.title(HISTTitle) plt.hist(array, bins=2500) plt.show() # ============================================================================= # def CurveFit(d,r,c,size,EGAIN): # # """Function to produce a fit against supplied data. The image (d(array)) multiplied by the EGAIN(a float), is then sliced into a box of 'size(int)'x'size' centre at row (r(int)) and column (c(int)) """ # #define fit function # def func(x,*theta): # theta0, theta1,theta2,theta3,theta4 = theta # s = theta0 #signal # b = theta1 #background # nu = theta2 #degrees of freedom # mu = theta3 #mean # sigma = theta4 #standard deviation # return s*gamma((nu+1.)/2.)/(np.sqrt(nu*np.pi)*sigma*gamma(nu/2.))*(1+(x-mu)**2/(nu*sigma**2.))**(-(nu+1)/2.)+b*1/(x[-1]-x[0]) # # def FWHM(nu,sigma): # return 2*sigma*np.sqrt(nu*(2.**(2./(nu+1.))-1.)) # def FWHMg(sigma): # return 2*np.sqrt(2*np.log(2))*sigma # # #set of data # y = d[r,c-int(size/2.):c+int(size/2.)]*EGAIN # x = np.linspace(1,len(y),len(y)) # sig = np.sqrt(y) # #set default parameter values # p0= np.array([1.2e5,2222.,1.,50.,4.]) # # p0= np.array([1.2e5,2200.,54.,4.]) # numPoints = len(x) # numPar = len(p0) # ndof = numPoints-numPar # thetaHat, cov = curve_fit(func,x,y,p0,sig,absolute_sigma=True) # # #minimized chi-squared, dof etc. # numPoints = len(x) # numPar = len(p0) # ndof = numPoints-numPar # chisq = sum(((y-func(x,*thetaHat))/sig)**2) # print "chisq = ", chisq, ", ndof = ",ndof # # #print fit parameters # print "\n", "Fitted parameters and standard deviations:" # sigThetaHat = np.sqrt(np.diag(cov)) # for i in range(len(thetaHat)): # print "thetaHat[", i, "] = ", thetaHat[i], " +- ", sigThetaHat[i] # # # # print "FWHM(StudentsT) = ", FWHM(thetaHat[2],thetaHat[4]) # print "FWHM(Gaussian) = ", FWHMg(thetaHat[4]) # # #Set up plot # plt.clf() # fig, ax = plt.subplots(1,1) # plt.gcf().subplots_adjust(bottom=0.15) # plt.gcf().subplots_adjust(left=0.15) # plt.errorbar(x, y, yerr=sig, xerr=0, color='black', fmt='o',ms=2, label='data') # plt.xlabel(r'$x^{th}$ pixel') # plt.ylabel(r'Num of photoelectrons') # plt.title(r'Atlas') #Title # xMin=0 # xMax=len(y) # yMin=0 # yMax=np.max(d[r,c-int(size/2.):c+int(size/2.)]*EGAIN)+5000. # plt.xlim(xMin, xMax) # plt.ylim(yMin, yMax) # xPlot = np.linspace( xMin, xMax, 500) # # # plt.bar(x,y,alpha=0.5,color='b') # # fit = func(xPlot, *thetaHat) # plt.plot(xPlot, fit, 'red', linewidth=2, label='fit result') # ticks = [0,20,40,60,80,100] # plt.xticks(ticks,np.linspace(c-int(size/2.),c+int(size/2.)-1,6)) # # handles, labels = ax.get_legend_handles_labels() # handles = [handles[1], handles[0]] # labels = [labels[1], labels[0]] # handles = [handles[0][0], handles[1]] # plt.legend(handles, labels, loc='upper right', fontsize=14, frameon=False) # # plt.show() # =============================================================================