;Generate_PhaseScreen
;
;PURPOSE:
;	This code will return the wavefront slopes from a Phase Screen.
;	It uses either a finite centered, backward, or forward difference
;	scheme.  It also accepts an input of how many pixels to average over
;	in order to obtain the slope.

;INPUTS:
;	PhaseScreen:
;		An NxN array of floating point numbers
;	ApDim:
;		The dimension of the array that the slopes will be 
;		returned in.
;	DifferenceMeth:
;            1 = Finite Centered Difference
;            2 = Finite Forward Difference
;            3 = Finite Backward Difference
;	     4 = Compact Pade Implicit Scheme
;	ApNum:
;		The number of apertures to average over to get the 
;		slopes.
;
;KEYWORDS:
;	ShackHartman: 1 - We calculate an average slope for each ApNum x ApNum
;			  array and return an average for those
;		      0 - We return every slope calculated
; 
;OUTPUTS
;	SlopesFinal = ApDimxApDim array
;	SlopesFinal(*,*,0) = XSlopes
;	SlopesFinal(*,*,1) = YSlopes

function return_slopes, PhaseScreen, ApDim, DifferenceMeth, ApNum, $
	 ShackHartman = ShackHartman

  If not(keyword_set(ShackHartman)) then ShackHartman = 0

  PhaseSize=Size(PhaseScreen)		;NxN array
  N=PhaseSize(1)			;N is in Size element 1

  XSlopesF=fltarr(N-2,N-2)			;Full Array of Slopes
  YSlopesF=fltarr(N-2,N-2)			;**********
  XSlopesS=Fltarr(ApNum*ApDim,ApNum*ApDim)	;Sampled Array of Slopes
  YSlopesS=Fltarr(ApNum*ApDim,ApNum*ApDim)	;***********

  If ShackHartman then begin 			;Return Binned Slopes
    SlopesFinal=Fltarr(ApDim,ApDim,2)		
  Endif Else begin				;Return All Slopes
    SlopesFinal=Fltarr(ApNum*ApDim, ApNum*ApDim, 2) 
  EndElse

  if DifferenceMeth eq 1 then begin
    XSlopesF(*,*)=(PhaseScreen(2:N-1,1:N-2)-PhaseScreen(0:N-3,1:N-2))/2
    YSlopesF(*,*)=(PhaseScreen(1:N-2,2:N-1)-PhaseScreen(1:N-2,0:N-3))/2
  endif else if DifferenceMeth eq 2 then begin
    XSlopesF(*,*)=(PhaseScreen(1:N-1,1:N-2)-PhaseScreen(0:N-2,1:N-2))/2
    YSlopesF(*,*)=(PhaseScreen(1:N-2,1:N-1)-PhaseScreen(1:N-2,0:N-2))/2
  endif else if DifferenceMeth eq 3 then begin
    XSlopesF(*,*)=(PhaseScreen(0:N-2,1:N-2)-PhaseScreen(1:N-1,1:N-2))/2
    YSlopesF(*,*)=(PhaseScreen(1:N-2,0:N-2)-PhaseScreen(1:N-2,1:N-1))/2
  endif else if DifferenceMeth eq 4 then begin
    SlopesF=return_padecompact_deriv(PhaseScreen)
    XSlopesF(*,*)=SlopesF(*,*,0)
    YSlopesF(*,*)=SlopesF(*,*,1)
  endif
  XSlopesS(*,*)=XSlopesF(0:ApNum*ApDim-1,0:ApNum*ApDim-1)
  YSlopesS(*,*)=YSlopesF(0:ApNum*ApDim-1,0:ApNum*ApDim-1)

  If ShackHartman then begin
    SlopesFinal(*,*,0)=rebin(XSlopesS(*,*),ApDim,ApDim)
    SlopesFinal(*,*,1)=rebin(YSlopesS(*,*),ApDim,ApDim)
  EndIf else begin
    SlopesFinal(*,*,0)=XSlopesS(*,*)
    SlopesFinal(*,*,1)=YSlopesS(*,*)
  EndElse

  return,SlopesFinal

end