import numpy as np import matplotlib.pyplot as plt import pyfits from scipy.ndimage import median_filter from scipy import signal def weigave(x,y): return sum(x)/sum(y) def standev(x,xbar, n): return np.sqrt(np.sum((x-xbar)**2)/n) def PeakFinder(filename, cutoff, value): f= pyfits.open(filename) 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 Gain = h.get("EGAIN") # gets te gain value fgx=px*Gain # gives the number of photoelectrons 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() #ffgx= secdiffynorm ffgx= median_filter(secdiffynorm, size=4) xmer=xint[ffgx < cutoff] print xmer #----------------------------------------------------------------------------------- # Finding range for where peaks could be 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 #-----------------------------------------------------------------------------------# # Finding weighted average to give position of peaks 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]) #print weigave(multiple,bgsubmaxnorm[i]) #-----------------------------------------------------------------------------------# # Finding the sigma value of the curve using "RMS" equation 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= fgx[i]*(np.square(i-np.mean(i))) print 'multiple-->',multiple, '\n' print 'i->',i,'mean-->',np.mean(i) num.append(multiple) denom.append(fgx[i]) print 'figx->',fgx[i] for i in range(0,len(num)): if len(xy[i])>value: print sum(num[i]),'bot',sum(denom[i]) print 'sigma--->',weigave(num[i],denom[i]) #-----------------------------------------------------------------------------------# # Finding the sigma value of the curve using "standard deviation" equation 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) 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---->', standev(x[i],meana[i],totN[i]) #-----------------------------------------------------------------------------------# # Plotting the data plt.figure(1) plt.title('Normalised plot of spectrum of Vega at 06.50mm and its second derivative') plt.plot(bgsubmaxnorm) plt.plot(xint, ffgx) plt.xlabel('Pixel values') plt.ylabel('Normalised values') plt.figure(2) #plt.subplot(211) plt.title('Plot of spectrum of Vega at 06.50mm with located peaks') plt.xlabel('Pixel values') plt.ylabel('Photoelectrons') plt.plot(fgx, linewidth='2.0') for i in range(0,len(xy)): if len(xy[i])>value: # this statement is here to get rid of any peaks that could be due to hot pixels as their "width" is too small. plt.axvline(weigave(xy[i],y[i]), linewidth='2.0', color='red') plt.show()