import numpy as np 
import matplotlib.pyplot as plt
import pyfits
from scipy.ndimage import median_filter

def waverage(x,y):
    return (sum(x))/sum(y)

def standev(x,xbar,N):
    return np.sqrt(np.sum((x-xbar)**2)/(N))

f= pyfits.open("Calib_cd_6.50.fit")
d = f[0].data # 2d array data
h = f[0].header # fits header data
px = d.sum(0) #converts to 1D array as well as sum the vlaues
py = d.sum(1) # same a px but does the sum at differet axis
Gain = h.get("EGAIN") # gets te gain value
fgx=px*Gain # gives the number of photoelectrons

#fgx= median_filter(gx, size=4)#  uses a median filter to effectively remove hot pixels

secdiffy = np.diff(fgx,2)
xint = np.arange(0, 763, 1)

bgsubmaxnorm = (fgx[0:-2] - fgx[0:-2].mean())/ fgx[0:-2].max()
secdiffynorm = secdiffy / secdiffy.max()

xmer=xint[secdiffynorm < -0.1]

list_of_lists=[] #creates an empty array
sub_list=[]
sub_list.append(xmer[0]) #include the first elememt of the array

for i in range(0,len(xmer)-1): # loops over from 0 to the lenght of xmer -1, main function of this loop is to form a list of lists which groups of peak values
    if xmer[i]+1==xmer[i+1]: # if the current element +1 is equal to the next element goes into this statement
        sub_list.append(xmer[i+1]) # adds the value of the next element to a list
    else: # however if this element's value +1 does not equal the next element then goes through this statement
        if len(sub_list)>0:
            list_of_lists.append(sub_list)# adds all the values of the sublist into list_of_lists
        sub_list = [] # starts a new sub_list
        sub_list.append(xmer[i+1])

list_of_lists.append(sub_list)#adds all of the sub_list into list_of_lists 
print list_of_lists

xy=[]
y=[]

for i in list_of_lists:# goes through groups of peak values and finds the weighted average for each group.
    #print 'i------>', i
    #print 'px------>', px[i]
    multiple= bgsubmaxnorm[i]*i
    xy.append(multiple)
    y.append(bgsubmaxnorm[i])

num=[]
denom=[]
for i in list_of_lists:#This is a loop for trying to find the rms of each of the group to give the width/std of the curve.
    multiple= bgsubmaxnorm[i]*(np.square(i))
    num.append(multiple)
    denom.append(bgsubmaxnorm[i]) 

for i in range(0,len(num)):
    print waverage(num[i],denom[i])# gives a value too large to be the width so different equation used. 

x=[]
meana=[]
totN=[]
for i in list_of_lists:# this obtains the variables to use for the standard deviation equation for the width.
    x.append(i)
    meana.append(np.mean(i))
    #totN.append(np.sum(bgsubmaxnorm[i]))
    totN.append(len(i)-1)

print 'x--->', list_of_lists    
for i in range(0,len(x)):# goes through the list of values obtained in the pervious for loop and calculated the standard deviation of each.
    if totN[i]!=0:
        if standev(x[i],meana[i],totN[i])!=0: 
            print 'std---->', 2*2.355*standev(x[i],meana[i],totN[i])
    
    
#print 'x---->', x
#print 'm---->', meana
#print 't---->', totN
plt.figure(1)
plt.title('Camera angle for 6.5mm(without median filter, normalised values)')
plt.plot(bgsubmaxnorm)
plt.plot(xint, secdiffynorm)
plt.xlabel('Pixel values')
plt.ylabel('Normalised values')
plt.figure(2)
#plt.subplot(211)
plt.title('Camera angle for 6.5mm with peaks located (without median filter)')
plt.xlabel('Pixel values')
plt.ylabel('Photoelectrons')
plt.plot(fgx, linewidth='2.0')
for i in range(0,len(xy)):
    #if len(xy[i])>1: # this statement is here to get rid of any peaks that could be due to hot pixels as their "width" is too small.
        #print 'waverage', waverage(xy[i],y[i])
        plt.axvline(waverage(xy[i],y[i]), linewidth='2.0', color='red')
plt.show()