Pro boxanalysesoarimage,date=date,SingleImage=SingleImage,$
	       threshold=threshold,Bordersize=Bordersize,$
	       SearchRegion=Searchregion,VR=VR
If (not keyword_set(SingleImage)) then SingleImage = 1
If (not keyword_set(date)) then date = '11'
If (not keyword_set(Threshold)) then Threshold = 1500
If (not keyword_set(SubBack)) then SubBack = 0
If (not keyword_set(SearchRegion)) then SearchRegion = 5
If (not keyword_set(VR)) then VR = 7
If (not keyword_set(BorderSize)) then BorderSize = VR+1
If (not keyword_set(Pitch)) then Pitch = 15
If (not keyword_set(Thickness)) then Thickness = 40
;PRIMITIVE SUBTRACTION AND FLAT FIELD
CalibrationDir='/nfs/slac/g/ki/ki08/lsst/CPcampaign/SOAR/2005-05-09/'+$
	      'Calibration/Images/'
BiasFrameFile=CalibrationDir+'Zero09May2005.fits'
fits_read,BiasFrameFile,Bias,BiasHeader
FlatFrameFile=CalibrationDir+'SFlat_09May2005R.fits'
fits_read,FlatFrameFile,Flat,FlatHeader
FlatNormal=float(Flat)/max(Flat)
ZeroFrameFile='/nfs/slac/g/ki/ki08/lsst/CPcampaign/SOAR/2005-05-09/'+$
              'Raw/zero.1024.fits'
fits_read,ZeroFrameFile,ZeroFrame,ZeroHeader
;SOAR IMAGES IN RESPECTIVE 2005-05-{9,10,11,12} DIRECTORIES
SOARDir='/nfs/slac/g/ki/ki08/lsst/CPcampaign/SOAR/2005-05-'+Date+$
	'/Raw/'
SOARFiles=file_search(SOARDir+'*.fits',count=NumSOARFiles)
Naxis1=SXPAR(ZeroHeader,'NAXIS1')
Naxis2=SXPAR(ZeroHeader,'NAXIS2')
Im=fltarr(Naxis1,Naxis2)

for i=0, 143 do begin; NumSOARFiles-2 do begin
  fits_read,SOARFiles(i),Im1,header1,/header_only
  ExpTime=SXPAR(header1,'EXPTIME')
  WindDir=SXPAR(header1,'WINDDIR')
  TimeObs=SXPAR(header1,'TIME-OBS')
  TelFocus=SXPAR(header1,'TELFOCUS')

  ;Get rid of nasty pixels
  If ExpTime eq 30.0 then begin
     fits_read,SOARFiles(i),Im1,header1
     fits_read,SOARFiles(i+1),Im2,header2
     Diff=float(Im2)-float(Im1)
     DiffHist=histogram(Diff,Nbins=40,locations=locDiffHist)
     ImHist=histogram(Im1,Nbins=40,locations=locImHist)

     SatPix=where(Im1 gt 50000)
     Im1(SatPix)=0                  ;BAD PIXELS
     Im1(0:1536)=0                  ;TWO BAD ROWS

     Hits=where(Im1 gt threshold)
     Im(Hits) = Im1(Hits)
  endif
endfor

Im1=Im
Im1MB=Im1(BorderSize-1:Naxis1-1-Bordersize,BorderSize-1:Naxis2-1-BorderSize)

Pass = 0
for i=-SearchRegion,SearchRegion do begin
  for j=-SearchRegion,SearchRegion do begin
    if not ((i eq 0) and (j eq 0)) then begin
      FirstCut=where(Im1MB gt Im1(Bordersize-1-i:Naxis1-1-BorderSize-i,$
                            BorderSize-1-j:Naxis2-1-BorderSize-j))
      if Pass eq 0 then begin
        Highest=FirstCut
      endif else begin
        Highest=[Highest,Firstcut]
        Highest=Highest(sort(Highest))
        Dups=where(Highest(0:N_elements(Highest)-2) eq Highest(1:N_elements(Highest)-1))
        Highest=Highest(Dups)
      endelse
      Pass=Pass+1
    endif
  endfor
endfor

;/////////////////////////////////////////////////////////////////////
NumHits=N_elements(Highest)
XCosmics=(Highest mod (Naxis1+1-2*BorderSize))+BorderSize-1
YCosmics=(Highest/(Naxis1+1-2*BorderSize))+BorderSize-1
XSep=XCosmics(1:NumHits-1)-XCosmics(0:NumHits-2)
YSep=YCosmics(1:NumHits-1)-YCosmics(0:NumHits-2)
PossibleDups=where(abs(Xsep) lt (VR+1) and abs(Ysep) lt (VR+1))
;/////////////////////////////////////////////////////////////////////
;FILL ARRAYS WITH DATA FROM COSMIC RAYS
Lambda1Arr=fltarr(NumHits)        ;The long second moment
Lambda2Arr=fltarr(NumHits)        ;The short second moment
EArr=fltarr(NumHits)              ;Scalar Ellipticity
FluxArr=fltarr(NumHits)           ;Total Counts in a region
TrackLength2DArr=fltarr(NumHits)  ;Length of the track in Pixels 2D
TrackLength3DArr=fltarr(NumHits)  ;Lenght of the track in physical units
PerpCountsArr=fltarr(NumHits)     ;Counts scaled by length
HitSizeArr=fltarr(NumHits)        ;Total number of pixels affected
HitWidthArr=fltarr(NumHits)       ;HitWidth as judged by Length
XExtentArr=Uintarr(NumHits)       ;Extent of Hit in X pixels
YExtentArr=Uintarr(NumHits)       ;Extent of Hit in Y pixels
FillFactorArr=fltarr(NumHits)     ;Amount of space filled in hit
   
NoSpots=0                         ;Number of "Spots" 
NoMuons=0                         ;Number of "Muons"
NoWorms=0                         ;Number of "Worms"
NoPossStars=0                     ;Number of "Possible Star" Candidates
NoDeltas=0                        ;Number of "Delta Rays"
;/////////////////////////////////////////////////////////////////////
for Cosmics=0, NumHits-1 do begin
  ;Return Coordinate in overall Image
  XCosmic=(Highest(Cosmics) mod (Naxis1+1-2*BorderSize))+BorderSize-1
  YCosmic=(Highest(Cosmics)/(Naxis1+1-2*BorderSize))+BorderSize-1
  ;Create Analysis Region
  XRegion=rebin(findgen(2*VR+1)-VR,2*VR+1,2*VR+1)
  YRegion=rebin(transpose(findgen(2*VR+1)-VR),2*VR+1,2*VR+1)
  IRegion=Im1(XCosmic-VR:XCosmic+VR,YCosmic-VR:YCosmic+VR)
  IRegionNorm=float(IRegion)/total(IRegion)
  ;Get Average X and Y Positions.  Redefine Analysis Region around them.
  XAvgReg=fix(total(XRegion*IRegionNorm))
  YAvgReg=fix(total(YRegion*IRegionNorm))
  IRegion=Im1(XCosmic-VR+XAvgReg:XCosmic+VR+XAvgReg,$
              YCosmic-VR+YAvgReg:YCosmic+VR+YAvgReg)
  FluxArr(Cosmics)=total(IRegion)
  IRegionNorm=float(IRegion)/FluxArr(Cosmics)
  ;Calculate Second Moments and find them in rotated frame
  I11=total(XRegion^2*IRegionNorm)
  I22=total(YRegion^2*IregionNorm)
  I12=total(XRegion*YRegion*IRegionNorm)
  E1=(I11-I22)/(I11+I22)
  E2=(2*I12)/(I11+I22)
  XMat=[[I11,I12],[I12,I22]]
  Evals=Hqr(Elmhes(XMat),/Double)
  Lambda1Arr(Cosmics)=Real_part(max(Evals))
  Lambda2Arr(Cosmics)=Real_part(min(Evals))
  EArr(Cosmics)=E1^2+E2^2
  ;Get Track Length and Number of Perpendicular Counts
  XsCosmic=where(IRegionNorm ne 0) mod (2*VR+1)
  YsCosmic=where(IRegionNorm ne 0)/ (2*VR+1)
  XExtentArr(Cosmics)=Max(XsCosmic)-Min(XsCosmic)
  YExtentArr(Cosmics)=Max(YsCosmic)-Min(YsCosmic)
  TrackLength3DArr(Cosmics)=sqrt(Pitch*XExtentArr(Cosmics)^2+$
                                 Pitch*YExtentArr(Cosmics)^2+$
                                Thickness^2)
  TrackLength2DArr(Cosmics)=sqrt(XExtentArr(Cosmics)^2+$
                                 YExtentArr(Cosmics)^2)+1
  PerpCountsArr(Cosmics)=FluxArr(Cosmics)*Thickness/TrackLength3DArr(Cosmics)
  HitPixels=where(IRegionNorm gt 0)
  HitSizeArr(Cosmics)=N_elements(HitPixels)
  HitWidthArr(Cosmics)=HitSizeArr(Cosmics)/TrackLength2DArr(Cosmics)
  FillFactorArr(Cosmics)=HitSizeArr(Cosmics)/$
                         float((1+XExtentArr(Cosmics))*(1+YExtentArr(Cosmics)))
  ;Try to find gaps in the pixel arrays HS=Horizontally Sorted VS=Vertically
  PossibleBreaks=-1
  PossibleDeltas=-1
  DefiniteBreaks=-1
  if N_elements(HitPixels) gt 2 then begin
   HitPixels=HitPixels(sort(HitPixels))
   CosmicsV=(XsCosmic-1)*(2*VR+1)+YsCosmic
   XsCosmicHS=XsCosmic(sort(HitPixels))
   YsCosmicHS=YsCosmic(sort(HitPixels))
   XsCosmicVS=XsCosmic(sort(CosmicsV))
   YsCosmicVS=YsCosmic(sort(CosmicsV))
   XsLines=XsCosmic(sort(XsCosmic))
   YsLines=YsCosmic(sort(YsCosmic))
   XsLineBreaks=XsLines(1:HitSizeArr(Cosmics)-1)-$
		XsLines(0:HitSizeArr(Cosmics)-2)
   YsLineBreaks=YsLines(1:HitSizeArr(Cosmics)-2)-$
		YsLines(0:HitSizeArr(Cosmics)-1)
   XsBreakH=XsCosmicHS(1:HitSizeArr(Cosmics)-1)-$
           XsCosmicHS(0:HitSizeArr(Cosmics)-2)
   YsBreakH=YsCosmicHS(1:HitSizeArr(Cosmics)-1)-$
           YsCosmicHS(0:HitSizeArr(Cosmics)-2)
   XsBreakV=XsCosmicVS(1:HitSizeArr(Cosmics)-1)-$
           XsCosmicVS(0:HitSizeArr(Cosmics)-2)
   YsBreakV=YsCosmicVS(1:HitSizeArr(Cosmics)-1)-$
           YsCosmicVS(0:HitSizeArr(Cosmics)-2)
   DefiniteBreaks=where((XsLineBreaks gt 1) or (YsLineBreaks gt 1))
   PossibleBreaksH=where((abs(XsBreakH) gt 1) and (abs(YsBreakH) gt 1))
   PossibleBreaksV=where((abs(XsBreakV) gt 1) and (abs(YsBreakV) gt 1))
    
   PossibleDeltas=where((((abs(XsBreakH) gt 1) and (abs(YsBreakH) eq 0)) or $
                         ((abs(XsBreakH) eq 0) and (abs(YsBreakH) gt 1))) or $
                        (((abs(XsBreakV) gt 1) and (abs(YsBreakV) eq 0)) or $
                         ((abs(XsBreakV) eq 0) and (abs(YsBreakV) gt 1))))

  endif 
  ;Diagnose Type of Hit: Spot, Worm, or Delta
  if DefiniteBreaks(0) ne -1 then begin
     print,"Double Hit in Viewing Region"
     window,xsize=512,ysize=512
     tvscl,congrid(alog(IRegion),512,512,/center)
     print,XExtentArr(Cosmics),YExtentArr(Cosmics),FillFactorArr(Cosmics)
     print,Earr(Cosmics)
     PossibleDeltas=-1
     stop
  
  endif else if PossibleDeltas(0) ne -1 then begin
     NoDeltas=NoDeltas+1
     window,xsize=512,ysize=512 
     tvscl,congrid(alog(IRegion),512,512,/center)
     print,XExtentArr(Cosmics),YExtentArr(Cosmics),FillFactorArr(Cosmics)
     print,Earr(Cosmics)
     PossibleDeltas=-1
     stop
  endif else if TrackLength2DArr(Cosmics) eq 1 then begin
     NoSpots=NoSpots+1
  endif else if ((FillFactorArr(Cosmics) gt .5) and $
                 (EArr(Cosmics) lt .15)) then begin
     NoWorms=NoWorms+1
  endif 

endfor

;////////////////////////////////////////////////////////////////////////
;HISTOGRAMS AND SUCH

!p.multi=[0,2,2]
HitSizeHist=histogram(HitSizeArr,Nbins=25,max=25,locations=locHitSize)
plot,locHitSize,HitSizeHist,psym=10,title='Num pixels hit'
TrackLength3DHist=histogram(TrackLength3DArr,Nbins=10,$
                          locations=locTrack3DLength)
plot,locTrack3DLength,TrackLength3DHist,psym=10,title='TrackLength'
TrackLength2DHist=histogram(TrackLength2DArr,Nbins=15,$
                          locations=locTrack2DLength)
plot,locTrack2DLength,TrackLength2DHist,psym=10,title='TrackLength 2D'
HitWidthHist=histogram(HitWidthArr,Nbins=15,$
                          locations=locHitWidth)
plot,locHitWidth,HitWidthHist,psym=10,title='Track Width'
stop
lambda1Hist=histogram(Lambda1Arr,Nbins=10,$
                      locations=locLambda1)
plot,locLambda1,Lambda1Hist,psym=10,title='Lambda 1'
lambda2Hist=histogram(Lambda2Arr,Nbins=10,$
                      locations=locLambda2)
plot,locLambda2,Lambda2Hist,psym=10,title='Lambda 2'
lambda2Hist=histogram(Lambda2Arr,Nbins=10,$
                      locations=locLambda2)
plot,Lambda1Arr,Lambda2Arr,psym=2,xtitle='Lambda 1',ytitle='Lambda 2'
plot,lambda2arr,perpcountsarr*10^(-4.),psym=2,yrange=[0,2],xrange=[0,2],$
                      xtitle='Lambda 2',ytitle='Perp Counts'

stop

end