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()
# =============================================================================