;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

data=mrdfits('../CMOS/dark_160K_33_min_1_8.0_40.fits',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(1)-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=gauss2dfit(gauss_area,gauss_params,x,y)
				print,gauss_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

END