;return_soarellipticity
;return_soarellipticity
;
;PURPOSE:
;       This script takes a 25x25 array of floating point numbers that represent
;       the tilts of a wavefront across the ShackHartmann wavefront sensor on 
;       the SOAR telescope and returns the ellipticity of a star that would 
;       be formed by the wavefront in the case of geometric optics.
;
;INPUTS:
;	XSlopes-A 25x25 array with the xslope values
;	YSlopes-A 25x25 array with the yslope values

function return_covariance,XSlopes,YSlopes
  Num=25
  Cxl=fltarr(num)
  Cxt=fltarr(num)
  Cyl=fltarr(num)
  Cyt=fltarr(num)
  Covariances=fltarr(4,num)
  Index=findgen(625)
  XSlopesVer=XSlopes((index mod 25)*25+index/25)
  YSlopesVer=YSlopes((index mod 25)*25+index/25)
  X_off=0
  Y_off=0
  for spacing=0, num-1 do begin

    CxlArr=fltarr(num*num)
    CxtArr=fltarr(num*num)
    CylArr=fltarr(num*num)
    CytArr=fltarr(num*num)

    ;SCAN FIRST ALONG THE ROWS

    for column_index=0, num-1 do begin
       for row_index=0, num-1 do begin
          spot_index=column_index*num+row_index

          ;If we have reached the edge, start on the next row
          If (row_index+spacing) ge num then break
            ;Check to see if both have valid values
            If (XSlopes(spot_index) ne -10) and $
               (XSlopes(spot_index+spacing) ne -10) $
               then begin
               ;add offset in to tweek the plots
                    CxlArr(spot_index)=$
                    (XSlopes(spot_index)+x_off)*$
                    (XSlopes(spot_index+spacing)+x_off)
                    CxtArr(spot_index)=$
                    (YSlopes(spot_index)+y_off)*$
                    (Yslopes(spot_index+spacing)+y_off)
             endif
        endfor
     endfor

     Correlated=where(CxlArr ne 0)
     If Correlated(0) ne -1 then begin
        Cxl(spacing)=mean(CxlArr(Correlated))
        Cxt(spacing)=mean(CxtArr(Correlated))
        endif else begin
          print,'no correlations',spacing
          Cxl(spacing)=-10
          Cxt(spacing)=-10
        endelse
     for column_index=0, num-1 do begin
        for row_index=0, num -1 do begin
           spot_index=column_index*num+row_index

     ;If we have reached the edge, start on the next row
     If (row_index+spacing) ge num then break
     ;Check to see if both have valid values
         If (YSlopesVer(spot_index) ne -10) and $
            (YSlopesVer(spot_index+spacing) ne -10) $
            then begin
              ;add offset in to tweek the plots
              CylArr(spot_index)=$
                    (YSlopesVer(spot_index)+x_off)*$ 
                    (YSlopesVer(spot_index+spacing)+x_off)
              CytArr(spot_index)=$
                    (XSlopesVer(spot_index)+y_off)*$
                    (XslopesVer(spot_index+spacing)+y_off)
         endif
       endfor
     endfor

     Correlated=where(CylArr ne 0)
     If Correlated(0) ne -1 then begin
           Cyl(spacing)=mean(CylArr(Correlated))
           Cyt(spacing)=mean(CytArr(Correlated))
     endif else begin
           print,'no correlations',spacing
           Cyl(spacing)=-10
           Cyt(spacing)=-10
     endelse

endfor
Covariances(0,*)=Cxl
Covariances(1,*)=Cxt
Covariances(2,*)=Cyl
Covariances(3,*)=Cyt

return,Covariances

end