;This Routine is intended to search through a 33 read fits file and find any cosmic rays, whereby it should attempt to fit a 2-d Gaussian to the pixel distribution and return the parameters of that Gaussian


;READ, "Enter Name of File:",File_name,format=string
;READ, "Enter Row:",row,format=long
;READ, "Enter Read:",Read,format=long

openw,1,"PSF.txt",/append
;header=string(how are you!C)

printf,1,FORMAT='(%"Background\tAmplitude\tX Pixel\t\tY Pixel\t\tSigmax\t\tSigmay\t\tSigmar\t\tFile_No")'
printf,1 ,"-----------------------------------------------------------------------------------------------------------------------------"

Voltage=string(2,format='(I1)')
print,Voltage

for file_index=1,40 do begin

	if file_index lt 10 then name='_0'+string(file_index,format='(I1)')
	if file_index ge 10 then name='_'+string(file_index,format='(I2)')

	filename='./dark_160K_33_min_1_'+Voltage+'.0'+name+'.fits'
	print,filename


data=mrdfits(filename,0,My_header)
fldata=float(data)
poss_cos=0
poss_cos_oned=0
cum_sum=fltarr(1024,1024)
difference=fltarr(1024,1024)
my_mean=fltarr(1024,1024)
pixel_check=fltarr(33)

for i=5,32 do begin

;Add up all frames to get cumulative sum
cum_sum=cum_sum+fldata(*,*,i)
	if i ne 32 then begin
		difference=fldata(*,*,(i+1))-fldata(*,*,i)
		poss_cos_oned=where(difference gt 10000 and fldata(*,*,i) gt 0)
		;negatives=where(difference lt 0)
		;If the difference between one read and the last is greater than 1000		     ;0 fit a Gaussian to that pixel's surroundings

		if poss_cos_oned(0) ne -1 then begin
			poss_cos=array_indices(difference,poss_cos_oned)
			my_size=size(poss_cos)
			if my_size(0) eq 1 then begin
				length=0
				endif else begin 
				length=my_size(2)-1 
			endelse
			for j=0,length do begin
				
				;Define area to fit Gaussian to
				xstart=poss_cos(0,j)-5
				ystart=poss_cos(1,j)-5
				xfinish=xstart+9
				yfinish=ystart+9
				x=findgen(10)+xstart	
				y=findgen(10)+ystart
				gauss_area=fldata(xstart:xfinish,ystart:yfinish,32)
				
				my_gauss_elip=mpfit2dpeak(gauss_area,gauss_elip_params,x,y)
				my_gauss_circ=mpfit2dpeak(gauss_area,gauss_circ_params,x,y,/circular)
				;printf,1,gauss_elip_params(0),gauss_elip_params(1),gauss_elip_params(2),gauss_elip_params(3),gauss_circ_params(2),name, FORMAT='(3F0)'
				printf,1,FORMAT='(%"%d\t\t%d\t\t%4.0f\t\t%4.0f\t\t%6.4f\t\t%6.4f\t\t%6.4f\t\t%s")',gauss_elip_params(0),gauss_elip_params(1),gauss_elip_params(4),gauss_elip_params(5),gauss_elip_params(2),gauss_elip_params(3),gauss_circ_params(2),name

				print,gauss_elip_params
				print,gauss_circ_params
				endfor
						

		endif

	endif 

;Keep track of one pixel over 33 reads to see what's happening
;pixel_check(i)=fldata(910,286,i)

endfor

;my_mean=cum_sum/33
;diff=fldata(*,*,32)-my_mean(*,*)


;mwrfits,mean,'../CMOS/mean.fits',My_header
;mwrfits,diff,'../CMOS/diff.fits',My_header

;print,My_header
;print,pixel_check
;print,poss_cos

;b=data(row,*,read)

;plot,b


endfor
close,1

END