function chi_Resid_Corr, image_deconv, gradient, LAG_AUTO_CORR=Lag_ac, $
						 RATIO=ratio, SPECTRUM=spectrum
; Used by procedure Max_Resid_Like,
; Based on:
;	"Incorporation of Spatial Information in Bayesian Image Reconstruction:
;			the Maximum Residual Likelihood (MRL) Criterion",
; by Robert Pina and Richard Puetter UCSD 1992.
; coded by Frank Varosi NASA/GSFC 1992.
; F.V. 1992, added new positivity constraint function as suggested by R.P.,
;	so that final deconv image result = positivity( image_deconv ).

  common Deconv_Data, imaged, psf, psf_FT
  common Max_Resid_Like, close, phase2
  common chi_Resid_Corr, sigma, Lagac, mode
  common chi_Resid_Corr1, epsilon
  common chi_Resid_Corr2, resid, rac

	resid = imaged - convolve( positivity( image_deconv, EPS=epsilon ), $
							psf, FT_PSF = psf_FT )
	if N_elements( close ) NE 1 then close=0
	if N_elements( sigma ) NE 1 then sigma=0

	if N_elements( Lag_ac ) GT 0 then begin
		Lag = Lag_ac(0)
		if N_elements( Lag_ac ) GT 1 then Lagmin=Lag_ac(1) else Lagmin=1
	  endif else begin
		Lag = Lagac(0)
		if N_elements( Lagac ) GT 1 then Lagmin=Lagac(1) else Lagmin=1
	   endelse

	npix = N_elements( resid )
	sim = size( resid )
	Lag = (Lag > 0) < ( min( sim(1:2) )/2 -1 )
	Lagmin = ( Lagmin < Lag ) > 0

	if (Lag LE 3) then begin

	   rac = fltarr( 2*Lag+1, Lag+1 )

	   for L = Lagmin, Lag do begin
	    for i=1-L,L-1 do rac(Lag+i,L) = total( resid * shift( resid, i,L ) )
		for j=0,L do rac(Lag+L,j) = total( resid * shift( resid, L,j ) )
		for j=1,L do rac(Lag-L,j) = total( resid * shift( resid,-L,j ) )
	   endfor

	  endif else begin

		rac = convolve( resid, /AUTO_CORREL )

		s = sim/2 - Lag
		s(2) = s(2) + Lag
		L = sim/2 + Lag		;get the non-redundant elements only.

		rac = rac(s(1):L(1),s(2):L(2))
		rac(0:Lag-1,0) = 0

		if (Lagmin GT 0) then begin
			L1 = Lagmin-1
			rac( Lag-L1:Lag+L1, 0:L1 ) = 0
		   endif
	   endelse

	if (sigma GT 0) then begin
		sig2 = npix * sigma^2
		sig4 = npix * sigma^4
	  endif else begin
		sig = stdev( resid )
		sig2 = npix * sig^2
		sig4 = npix * sig^4
	   endelse

	if (Lagmin EQ 0) then rac(Lag,0) = (rac(Lag,0) - sig2)/sqrt(2)

	if keyword_set( spectrum ) then begin

		Lags = indgen( Lag ) + 1
		modes = (2*Lags+1)*(Lags+1) - Lags
		L = Lagmin-1
		modes = modes - ( (2*L+1)*(L+1) - L ) - 2
		chis = fltarr( Lag )
	    	for j=1,Lag do chis(j-1) = total( rac( Lag-j:Lag+j, 0:j )^2 )
		chis = chis/sig4
		if keyword_set( ratio ) then chis = chis/modes
		return, [ [modes], [chis] ]

	  endif else begin

		mode = (2*Lag+1)*(Lag+1) - Lag
		if (Lagmin GT 0) then begin
			L = Lagmin-1
			mode = mode - ( (2*L+1)*(L+1) - L )
		   endif
		mode = (mode-2) > 1
		chisq = total( rac * rac )/sig4
		if keyword_set( ratio ) then return, chisq/mode
	   endelse

	if N_params() GE 2 then begin

		resid = convolve( resid, FT_PSF=psf_FT, /CORREL )
		rac(Lag,0) = rac(Lag,0) / sqrt(2)
		gradient = fltarr( sim(1), sim(2) )

		for L = Lagmin, Lag do begin

			for i=1-L,L-1 do gradient = gradient + rac(Lag+i,L) * $
				( shift( resid, -i,-L ) + shift( resid, i, L ) )

			for j=0,L do gradient = gradient + rac(Lag+L,j) * $
				( shift( resid, -L,-j ) + shift( resid, L, j ) )

			for j=1,L do gradient = gradient + rac(Lag-L,j) * $
				( shift( resid, -L, j ) + shift( resid, L,-j ) )
		  endfor

		rac(Lag,0) = rac(Lag,0) * sqrt(2)
		if (close) then  fac=4*(chisq - mode)/sig4  else  fac=2/sig4
		gradient = -fac * positivity( image_deconv,/DERIV ) * gradient
	   endif

if (close) then  return,(chisq - mode)^2  else  return,chisq
end